A Targeted Next-Generation Sequencing Panel to Genotype Gliomas

Gliomas account for the majority of primary brain tumors. Glioblastoma is the most common and malignant type. Based on their extreme molecular heterogeneity, molecular markers can be used to classify gliomas and stratify patients into diagnostic, prognostic, and therapeutic clusters. In this work, we developed and validated a targeted next-generation sequencing (NGS) approach to analyze variants or chromosomal aberrations correlated with tumorigenesis and response to treatment in gliomas. Our targeted NGS analysis covered 13 glioma-related genes (ACVR1, ATRX, BRAF, CDKN2A, EGFR, H3F3A, HIST1H3B, HIST1H3C, IDH1, IDH2, P53, PDGFRA, PTEN), a 125 bp region of the TERT promoter, and 54 single nucleotide polymorphisms (SNPs) along chromosomes 1 and 19 for reliable assessment of their copy number alterations (CNAs). Our targeted NGS approach provided a portrait of gliomas’ molecular heterogeneity with high accuracy, specificity, and sensitivity in a single workflow, enabling the detection of variants associated with unfavorable outcomes, disease progression, and drug resistance. These preliminary results support its use in routine diagnostic neuropathology.


Introduction
Human gliomas are primary tumors of the central nervous system (CNS) that include different subtypes for their histopathology features and degree of malignancy. Glioblastoma is the most aggressive, invasive, and undifferentiated glioma with a median survival of 15 months [1]. In 90% of cases, it is a primary glioblastoma that arises without previous disease, while secondary glioblastoma is rare and develops from a lower grade astrocytoma [2,3]. Hardly histologically distinguishable, primary and secondary glioblastomas present different clinical features and affect patient groups of different ages and sexes. For example, secondary glioblastomas manifest predominantly in younger patients (median age of 45 years), with a better prognosis compared to primary glioblastomas (median age  [27]. In this article, we describe the development and validation of a custom targeted NGS approach for the analysis of 13 glioma-related genes (ACVR1, ATRX, BRAF, CDKN2A, EGFR, H3F3A, HIST1H3B, HIST1H3C, IDH1, IDH2, P53, PDGFRA, PTEN), a 125 bp region of the TERT promoter (TERTp), and 54 single nucleotide polymorphisms (SNPs) located along chromosomes 1 and 19 for the reliable assessment of their copy number alterations (CNAs).

Sample Collection
Tissue samples from patients who underwent surgical resection of newly diagnosed gliomas (n = 34) were collected at the Neurosurgery Unit of Fondazione IRCCS Ca' Granda  While enormous advances in the field of both molecular biology and genetics have elucidated the complex etiopathogenesis of gliomas and identified new potential treatment options that target molecular receptors and signaling pathways, only limited progress has been achieved in clinical practice to help the transition towards molecular medicine and personalized approaches. Several molecular signatures can influence tumor biology and impact clinical outcomes, offering the possibility to stratify patients into diagnostic, prognostic, and therapeutic clusters [4]. To this end, the application of different highthroughput technologies, such as NGS, may aid to enable the comprehensive detection of relevant genetic alterations for gliomas according to the WHO classification, as well as to identify those associated with drug resistance [25,26]. However, due to its complexity and challenging bioinformatics pipeline, NGS has not yet entered into routine diagnostic practice, fully replacing established low-resolution techniques [27].

Sample Collection
Tissue samples from patients who underwent surgical resection of newly diagnosed gliomas (n = 34) were collected at the Neurosurgery Unit of Fondazione IRCCS Ca' Granda Ospedale Maggiore Policlinico and processed at the Laboratory of Experimental Neurosurgery and Cellular Therapy of the same Institution. Histological subtypes were classified on the basis of morphologic and immunohistochemical characteristics according to the WHO 2016 classification criteria [28]. The clinical and molecular characteristics of the patients (described below in Section 3) were extracted from medical records and included sex,  R172T, c.515G > T p.R172M, c.516G > A p.R172R, c.516G > C p.R172S, c.516G > T p.R172S), and TERTp mutations (G228A and G250A) was performed using the MassARRAY Analyzer 4 system (Agena Bioscience, San Diego, CA, USA), based on MALDI-TOF mass spectrometry. Tumor proliferation was assessed by Ki-67 (MIB-1) staining using immunohistochemistry, in a range between 15% and 80%. Immunohistochemistry was performed using the BenchMark Ultra automatic system (Ventana Medical Systems) and was reviewed by two independent pathologists.
The Institutional Review Board approved the protocol (IRB#1670/2015), and all patients provided informed consent. All methods were performed in accordance with ethical regulations.
Notably, all patients followed Stupp therapeutic protocol consisting of maximal safe surgical resection, followed by radiotherapy (60 Gy) with concomitant daily TMZ chemotherapy, and then maintenance treatment with TMZ for 6 to 12 months. The approved conventional schedule was a daily dose of 120 to 200 mg per square meter of body-surface area for 5 days of every 28-day cycle [29].

DNA Extraction and Dosage
DNA from FFPE samples was extracted with a DNeasy Blood and Tissue Kits Isolation kit (Qiagen, Germantown, MD, USA), quantified by a NanoDrop ND-3000 spectrophotometer, and assessed for quality by microcapillary electrophoresis on a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The DNA yield was quantified using a TaqMan RNase P Detection Reagent Kit (Thermo Fisher Scientific, Waltham, MA, USA) on Real Time PCR AriaDx (Agilent Technologies, Santa Clara, CA, USA).

NGS Panel Design
A custom "on-demand" NGS panel was designed using custom IonAmpliSeq Designer software (https://ampliseq.com, accessed on 3 May 2020, Thermo Fisher Scientific, Waltham, MA, USA). The panel included 441 amplicons with a size range of 125-175 bp distributed across two primer pools (Pool_1: 221 amplicons, Pool_2: 220 amplicons) covering 47.22 kb and targeting the whole coding region of 13 glioma-related genes (ACVR1, ATRX, BRAF, CDKN2A, EGFR, H3F3A, HIST1H3B, HIST1H3C, IDH1, IDH2, P53, PDGFRA, PTEN), a 125 bp region of the TERT promoter, and 54 SNPs located along chromosomes 1 and 19 ( Table 1). The in silico coverage was 99.11%. The complete design of the panel is available in Supplementary Material S1 (Tables S1 and S2). Targeted genes and SNPs are detailed below and were selected according to the WHO diagnostic requirements for molecular testing of gliomas [12], a previous gene-set used for targeted strategies [30], and their high mutation frequency in gliomas.

Library Preparation
Automated library preparation was carried out with 0.67 ng of input DNA using an Ion AmpliSeq Kit for Chef DL8 (DNA to Library, 8 samples/run) on the Ion Chef System (Thermo Fisher Scientific, Waltham, MA, USA). Library quality and molarity were assessed by using an Ion Library TaqMan Quantitation Kit (Thermo Fisher Scientific, Waltham, MA, USA) on the Aria Dx Real-Time PCR System (Agilent Technologies, Santa Clara, CA, USA). Serial dilutions of the E. coli DH10B control library (Thermo Fisher Scientific, Waltham, MA, USA) were prepared and run in triplicate to generate a standard curve. Molar concentration of libraries was determined by using the Delta R baseline-corrected raw fluorescence calculated with Aria DX Real-Time PCR Software (Agilent Technologies, Santa Clara, CA, USA). Barcoded libraries were super-pooled in equimolar concentrations using the strategies suggested for combining libraries prepared with different panels for equal coverage to obtain final molarity of 60 pM each.

Chip Loading and Sequencing
Loading of the Ion 540 chips was carried out using an Ion 540 Kit-Chef (Thermo Fisher Scientific, Waltham, MA, USA) following manufacturer instructions. High throughput sequencing runs were carried out on an Ion Gene Studio S5 system (Thermo Fisher Scientific, Waltham, MA, USA). Run planned in the S5 Torrent Suite (v. 5.12.2) had the following parameters: analysis parameters, default; reference library, hg19; target regions, IAD199863.bed file; read length, 30X; flows, 550; base calibration mode, default. Plugins used were Coverage Analysis, Ion Reporter Uploader, and Variant Caller (so-matic_low_stringency_540).

Analytical Specificity and Sensitivity
To evaluate the performance of variant calls in detecting the true genotype for each sample, we estimated diagnostic accuracy, sensitivity, and specificity of the test by comparing our results with the previous data obtained using different technology (MassARRAY) [31]. True positives (TPs) were defined as correctly genotyped samples whose alternative variants were detected by both our filtering pipeline as well as from MassARRAY collected data. True negatives (TNs) were samples known to be wild-type (wt) from MassARRAY data. False positives (FPs) and false negatives (FNs) were considered as samples carrying variants detected by our pipeline but not from expected data and vice versa, respectively. Accuracy was calculated as follows: (TP + TN)/(TP + FP + TN + FN); sensitivity was calculated as follows: TP/(TP + FN); specificity was calculated as follows: TN/(TN + FP). The accuracy of NGS-based analysis for somatic genomic variants is also based on the correct detection of mutations at low variant allele frequencies (VAF) [32]. In the NGS method, this parameter is assessed by measuring the limit of detection (LOD). Different factors such as nucleic acid input, processing methods, or sequencing depth can influence LOD [33]. Moreover, in cancer screening, tissue heterogeneity caused by the co-existence of multiple mutations in tumor cells or the occurrence of multiple clones in the same tumor is another factor that can affect this metric [34]. To assess the accuracy of our method in the detection of somatic mutations at any VAF in a tumor sample, we measured the LOD using diluted DNA samples. Pre-validated genomic DNA samples were used as reference material and contained different types of mutations. The mix of known reference DNA dilutions with wild-type DNA was used as template, and a theoretical LOD ranging between 1 and 10% was assessed.

Variant Calling and Prioritization
Data obtained after sequencing runs were analyzed using the IonTorrent Suite v.5.12.3 workflow that allows the assessment of several parameters such as Key Signal, Test Fragment, Read Length, ISP Density, and other features related to sequencing run performance. Sequences were aligned to the hg19 reference genome, and variant calling was performed using Ion Reporter version 5.18 (Thermo Fisher Scientific, Waltham, MA, USA). Variants were filtered using the following filter chain: (a) filtered coverage ≥ 100; (b) 50 ≤ allele read-count ≤ 100,000; (c) location in exonic, intronic, upstream, splicesite_5, splicesite_3; (d) minimum variant allele frequency was set to 0.05 to call somatic mutations variants; (e) variant effect in refAllele, unknown, missense, nonframeshiftInsertion, nonframeshiftDeletion, nonframeshiftBlockSubstitution, nonsense, stoploss, frameshiftInsertion, frameshiftDeletion, frameshiftBlockSubstitution; (f) 0.0 ≤ minor allele frequency ≤ 0.5; (g) SIFT score less than 0.05 or PolyPhen-2 score greater than 0.85 were used as reference scores to understand the importance of deleterious variants. Filtered variants were manually reviewed in ExAc databases, ClinVar, and gnomAD tool as well as in the scientific literature to exclude polymorphisms or nonpathogenic variants. To detect 1p/19q LOH, we firstly applied a quality criterion based on the SNP coverage. The test was considered optimal, suboptimal, or noninformative according to the number of SNPs that were covered by fewer than 250 reads. Secondly, the allelic frequencies (AF) for each SNP (with more than 250×) were annotated. SNPs were defined homozygous when the AF was approximately 100% and led the same nucleotide as that of the reference genome or having an AF of approximately 0% with a nucleotide that differed from the reference genome. SNPs with AF of approximately 50% were considered heterozygous [30]. However, taking into consideration the semi-quantitative measurement provided by NGS, we established the following confidence intervals: 0-10% or 90-100% for homozygous markers, and 40-60% for heterozygous markers. Imbalances of 1p and 19q markers due to LOH were scored when their AFs were outside the established ranges for homozygosity or heterozygosity (i.e., 10-40% or 60-90%) [30]. To identify the 1p/19q codeletion, we used the criteria described elsewhere [35,36], which established LOH as the imbalance of all informative SNPs on both chromosomal arms. More specifically, variants with heterozygous frequencies in these chromosomal arms were considered not aberrant, and a 1p/19q codeletion was defined only when no heterozygous markers were present in each arm [30]. No codeletion was scored if at least one heterozygous marker was present in 1p or 19q. A whole chromosome arm LOH was observed when no heterozygous markers were present in either arm.

Results
In this work, we developed and validated a targeted NGS approach to analyze variants or chromosomal aberrations correlated with tumorigenesis and response to treatment in gliomas.

Study Population
The clinical and molecular parameters of patients enrolled in the study are collected in Table 2. The median age of patients was 64 years (IQR: 54-72), with the same percentage of females and males. The median MGMT promoter methylation was 23% (IQR: 7-42), six patients carried out IDH1 mutations, and seven patients showed a LOH. Notably, a value of MGMT promoter methylation > 9% is considered a favorable prognostic predictor, associated with a better response to therapy and therefore clinical outcome [37].

Target Gene Coverage
Targeted NGS analysis covered 13 genes (ACVR1, ATRX, BRAF, CDKN2A, EGFR, H3F3A, HIST1H3B, HIST1H3C, IDH1, IDH2, P53, PDGFRA, PTEN), a 125 bp region of the TERT promoter, and 54 SNPs located along chromosomes 1 and 19. Sequencing minimum depth was set to 100X, and, in this configuration, 97% of targeted regions achieved the minimum required cut-off. Amplicons with at least 100 reads were 98.7% and with 500 reads were 93.9%; total average coverage was 98.71% with 130 bp as mean read length and 99% of aligned reads to the hg19 reference genome. Adequate amplification efficiency was observed for 435/442 amplicons (98.4%) (mean assigned reads per amplicon Log10 ranging from three to five), while only seven amplicons (one for HIST1H3B, one for EGFR, one for BRAF, one for PTEN, one for IDH2, two for CDKN2A) had a low coverage value ( Figure 2). Targeted NGS analysis covered 13 genes (ACVR1, ATRX, BRAF, CDKN2A, EGFR,  H3F3A, HIST1H3B, HIST1H3C, IDH1, IDH2, P53, PDGFRA, PTEN), a 125 bp region of the TERT promoter, and 54 SNPs located along chromosomes 1 and 19. Sequencing minimum depth was set to 100X, and, in this configuration, 97% of targeted regions achieved the minimum required cut-off. Amplicons with at least 100 reads were 98.7% and with 500 reads were 93.9%; total average coverage was 98.71% with 130 bp as mean read length and 99% of aligned reads to the hg19 reference genome. Adequate amplification efficiency was observed for 435/442 amplicons (98.4%) (mean assigned reads per amplicon Log10 ranging from three to five), while only seven amplicons (one for HIST1H3B, one for EGFR, one for BRAF, one for PTEN, one for IDH2, two for CDKN2A) had a low coverage value ( Figure 2).

Figure 2.
Mean coverage of amplicons of target genes. The graph shows the mean coverage of each amplicon (indicated by "×") of 13 targeted genes in all analyzed samples. The # shows amplicons that do not satisfy a 100X depth cut-off (log10 = 2): AMPL7170424214 for HIST1H3B; AMPL7170319287 for EGFR, AMPL7170361379 for BRAF; AMPL7170348270 for PTEN; AMPL7170423161 for IDH2; AMPL7170325832 and AMPL7170330366 for CDKN2A (please refers to Supplementary Material S1, Table S2).

Variant Interpretation
The variant caller listed 2644 variants, including synonymous, unknown, missense, and indels. All variants were filtered to remove intronic variants and non-coding variants. Somatic synonymous mutations and those in the 5′ and 3′ untranslated regions (UTRs) were considered as functionally associated with cancer risk [38,39]. Using these filters, we identified 420 exonic, 1 splicesite, and 34 UTR variants, for a total of 455 alterations consisting of missense, synonymous, non-sense, nonframeshiftDeletion, and frameshiftInsertion. Based on variant allele frequency, disease prevalence, and in silico prediction score, variants were further filtered for "Pathogenic" and "Likely benign/Conflicting interpretations". These clinically relevant variants (48) are shown in the oncoplot (Figure 3) and in Supplementary Material S1, Table S3. The frequencies of genetic variants, together with those in the TCGA database, are shown in Supplementary Material S1, Table S4.
The presence of 1p/19q codeletion was assessed by a comprehensive analysis of 54 SNPs distributed across these chromosomes. As shown in Supplementary Material S2 Figure 2. Mean coverage of amplicons of target genes. The graph shows the mean coverage of each amplicon (indicated by "×") of 13 targeted genes in all analyzed samples. The # shows amplicons that do not satisfy a 100X depth cut-off (log10 = 2): AMPL7170424214 for HIST1H3B; AMPL7170319287 for EGFR, AMPL7170361379 for BRAF; AMPL7170348270 for PTEN; AMPL7170423161 for IDH2; AMPL7170325832 and AMPL7170330366 for CDKN2A (please refers to Supplementary Material S1, Table S2).

Variant Interpretation
The variant caller listed 2644 variants, including synonymous, unknown, missense, and indels. All variants were filtered to remove intronic variants and non-coding variants. Somatic synonymous mutations and those in the 5 and 3 untranslated regions (UTRs) were considered as functionally associated with cancer risk [38,39]. Using these filters, we identified 420 exonic, 1 splicesite, and 34 UTR variants, for a total of 455 alterations consisting of missense, synonymous, non-sense, nonframeshiftDeletion, and frameshiftInsertion. Based on variant allele frequency, disease prevalence, and in silico prediction score, variants were further filtered for "Pathogenic" and "Likely benign/Conflicting interpretations". These clinically relevant variants (48) are shown in the oncoplot ( Figure 3) and in Supplementary Material S1, Table S3. The frequencies of genetic variants, together with those in the TCGA database, are shown in Supplementary Material S1, Table S4.
The presence of 1p/19q codeletion was assessed by a comprehensive analysis of 54 SNPs distributed across these chromosomes. As shown in Supplementary Material S2 ( Figures S1 and S2), the amplicons fully covered the target regions of interest, with a good coverage value and a sufficient number of correctly mapped reads. Codeletion testing was based on the allelic frequency of each SNP, and the presence of multiple markers in 1p and 19q in the heterozygous form ruled out the presence of codeletion in these chromosomal arms. In total, we detected six samples with 1p/19q codeletion. Representative 1p/19q results are shown in Figure 4. Raw data of the tNGS are available at NCBI's Sequence Read Archive (SRA) with the accession number SUB11581843.    Table S3.

Assessment of Sensitivity and Specificity
In order to check the performance of both the panel and the analytic bioinformatic pipeline, we measured the diagnostic accuracy, sensitivity, and specificity of the method starting from the previously known genetic information obtained by MassARRAY ( Table 2). Sequencing of IDH1 revealed 6 correctly called TP samples, 28 TN samples, 0 FN, and 0 FP calls by comparing our results with expected data; thus, IDH1 sequencing resulted in an accuracy, sensitivity, and specificity of 100%. With regard to 1p/19q co-deletion, there were 6 correctly called TP samples, 27 TN samples, 1 FN (ID233), and 0 FP calls by comparing our results with expected data from MassArray, thus showing an accuracy of 97%, a sensitivity of 85%, and a specificity of 100% (please note that samples carrying only 1p or 19q deletion were considered as TN). With regard to variants located on the TERT promoter, both the two samples previously assessed by MassARRAY as carriers of TERT mutations (−146 C > T in the ID169 and the −124 C > T in the ID143, respectively) were confirmed by our workflow analysis. Moreover, the NGS panel detected 11 additional positive TERT promoter mutations in samples not previously analyzed with other technologies. For all called variants, the high proportion of concordance highlighted the high specificity of our method.
To evaluate the sensitivity of sequencing depth and derivate the LOD, DNA tumor samples with different genotype patterns were diluted with wild-type samples to create allelic frequencies between 1 and 10%. The sensitivity of variant detection was evaluated as the lowest frequency calling the somatic variant. As shown in Supplementary Material S2, for IDH1 (p.Arg132His variant), we observed a linear decrease of the allelic frequency with a lowest VAF of 2.99%. In similar experiments, we observed a VAF of 3.32% for TERT promoter (c.−146 C > T), 3.25% for PTEN (p.Arg335Ter variant), and 3.8% for TP53 (pArg282Trp variant). The variants were not called at a VAF < 2%. With regard to CNA detection limit, our method was able to identify the variant up to a 25% dilution (an example of Chr 19 deletion analysis is shown in Supplementary Material S2). Notably, in all conditions tested, the sequencing depth did not show significant variations (mean depth 1980X). Altogether, our results support the ability of our method to analyze samples with low tumor purity.

Discussion
The main objective of this study was to develop and validate a custom targeted NGS approach to identify genetic variations in key genes related to glioma tumorigenesis and treatment response. The preliminary results in terms of sequencing performances (Supplementary Material S1, Table S5), time, and costs provide encouraging evidence for the use of this approach in the clinical pipeline. This small panel size allowed for multiplexing of a large number of samples with adequate read depth. DNA input requirements were low (less than 1 ng of low-quality DNA), and the run time was considerably short (5-7 h for preparation of a pool of eight libraries with an automatic workflow and 8 h for data processing). We obtained high sequencing coverage, measured LOD, and accuracy of the method, and we detected a substantial number of clinically relevant variants in each analyzed sample. The amplicon-based method provided a detailed genetic characterization of genes currently required for molecular testing in gliomas. The method used in this study confirmed and extended previous molecular profiling performed by the MassARRAY Analyzer 4 system (Agena Bioscience, San Diego, CA, USA) ( Table 2 and Section 2).
Below, we introduce each of the genes included in the custom panel for their relevance in the diagnosis and management of patients with gliomas and briefly describe the identified genetic variants that are associated with adverse outcomes, disease progression, or drug resistance according to ClinVar (Supplementary Material S1, Table S3).
Isocitrate dehydrogenase genes (IDH1 and IDH2) are frequently associated with anaplastic astrocytic, oligodendroglial, and lower grade tumors [40]. IDH enzymes are involved in several cellular processes, including mitochondrial oxidative phosphorylation, glutamate metabolism, lipogenesis, regulation of glucose, DNA repair, and cellular energy balance [41]. The main heterozygous somatic point mutations are localized in the active sites of the both IDH enzymes. Using our panel, we confirmed the presence of pathogenic missense variant p.R132H in IDH1 in six patients. This point mutation in exon 4 involves arginine in position 132 and its homolog at codon 172 of IDH2 [42,43]. Deregulated IDH activity results in the inhibition of key enzymes involved in demethylation of histones and DNA, with hypermethylation, altered gene expression, and proliferation of mutant cells [44,45]. The prognostic significance of the IDH1 R132H mutation is inconsistent and was associated with tumor WHO grade. Several studies demonstrated an association between an IDH1 mutation and improved prognosis OS, whereas others showed no significant difference [46]. Several pieces of evidence have shown that IDH1/IDH2 mutational status is an important diagnostic and prognostic factor to predict the response to anti-IDH treatment [47]. In clinical trials, glioblastoma patients with IDH1 R132H treated with vandetanib or ivosidenib in combination with temozolomide showed a significant increase in overall survival (OS) compared to glioblastoma patients without IDH1 R132H [48][49][50].
EGFR belongs to the ErbB receptor family and is one of the major factors in glioma tumorigenesis [51]. EGFR plays a role in cell growth and proliferation, autophagy, and cell metabolism. Located on chromosomal band 7p12, known EGFR genetic alterations include point mutations, rearrangement, and copy number variations [52]. In gliomas, EGFR amplification promotes invasion, proliferation, and failure of radiotherapy and chemotherapy [53]. In our cohort, we identified the missense variant p.Gly719Ala in exon 18 with pathogenic/likely pathogenic clinical significance according to ClinVar (Supplementary Material S1, Table S3). In My Cancer Genomics, the presence of an EGFR-G719A mutation is cited as an inclusion criterion for glioblastoma clinical trials. In particular, a phase II trial involving patients with EGFR-activated glioblastoma is evaluating the effect of osimertinib on tumor cell growth by blocking enzymes required for cell growth [54,55].
Genetic alterations in the TERT promoter are the most frequent non-coding mutations and the earliest genetic events that occur in GBM development [56]. Germline and somatic alterations in the promoter region of this gene are point mutations, genomic rearrangements, and DNA amplifications that affect telomerase activity [57]. TERTp mutations, although not reported in ClinVar, are associated with a better prognosis in IDH-mutated gliomas, and a worse prognosis in patients with primary GBM [58]. In our cohort, two hotspot mutations, c.−124C > T and c.−146C > T, were identified (Supplementary Material S1, Table S3). These somatic mutations are known to increase TERT expression, telomerase enzyme activity, telomere length, and tumor formation. Although TERT mutation prevalence is over 80% in gliomas, their role as a prognostic/predictive biomarker is still largely controversial [59]. Numerous studies indicated the TERT promoter mutation as an independent factor for poor prognosis, and various approaches to target TERT activity, such as inhibitors, immunotherapy, and vaccines, are currently being explored [60]. Several clinical trials are evaluating the efficacy of surgery and temozolomide chemotherapy in combination with TERT inhibitors such as eribulin, imetelstat, and BIBR1532 [56,60].
TP53 is one of the most deregulated genes in neoplasms and is involved in cell cycle and invasion; response to DNA damage and immune response; cell differentiation; and proliferation, apoptosis, and genomic stability [61]. The high expression of gain of function (GOF) oncogenic variants of the p53 protein is associated with survival outcomes and with a more invasive and aggressive phenotype [61]. In our cohort, we identified several missense variants in TP53 with pathogenic/likely pathogenic clinical significance (Supplementary Material S1, Table S3) and some previously associated with gliomas [62][63][64].
ATRX encodes a chromatin-remodeling protein whose main function correlates with gene expression regulation [65]. ATRX mutations are strongly linked to DNA damage and apoptosis. Its alterations includes point mutations or deletions and are correlated with the Alternative Lengthening of Telomeres (ALT) phenotype, TP53 mutations, MGMT methylation status, and 1p/19q co-deletion [66]. In particular, loss of ATRX leads to increased genome instability and micronuclei formation, which can influence the response to therapy and the prognosis [67,68].
PTEN is a tumor suppressor gene that plays an important role in regulating cell cycle progression, cell migration, apoptosis, and DNA damage repair [69]. In gliomas, loss of PTEN expression has been associated with increased malignancy, metastasis, radiotherapy, and chemotherapy response [70]. In our cohort of patients, we identified four PTEN variants (c.697 C > T, c.1003 C > T, c.822G > A, and c.972delT), which are predicted to cause loss of normal protein function and therefore are classified as pathogenic and disease-causing (Supplementary Material S1, Table S3).
PDGFRA encodes for a cell-surface transmembrane receptor that controls many important cellular processes such as normal proliferation, growth, proliferation, and survival of glial cells [71,72]. Several genetic rearrangements in PDGFRA have been identified in gliomas, but their diagnostic or prognostic significance is unclear [73].
Activin A receptor type I (ACVR1) encodes a bone morphogenetic receptor. Its function has been linked to many biological processes and, recently, to tumorigenesis [74]. Somatic mutations in ACVR1 have been associated with high-grade pediatric brainstem gliomas, but their specific role remains to be elucidated [75].
BRAF plays key roles in proliferation, differentiation, cell motility, and apoptosis [76]. Based on these findings, the use of BRAF as a therapeutic target has been a milestone [77]. One of targetable mutations is BRAF, i.e., p.V600E, which occurs in roughly half of all epithelioid gliomas but is rare in classic GBM [78]. Detection of the BRAF V600E mutation represents a positive prognostic indicator [78].
Mutations in H3F3A and in HIST1H3B/C genes play a pivotal role in development of diffuse midline glioma, as well as a subset of adult glioblastoma multiforme [79]. The recurrent K27M mutations in both genes arising in the thalamus, brainstem, and spinal cord are associated with poor prognosis regardless of histologic grade [80]. In contrast, p.G34R or p.G34V mutations in H3F3A are found in the cerebral cortex of adolescent and young adults and are associated with a better prognosis [81].
Complete deletion of both the short arm of chromosome 1 (1p) and the long arm of chromosome 19 (19q) (1p/19q-codeletion) is a genetic signature occurring in the pathogenesis of oligodendrogliomas [82]. This 1p/19q co-deletion is a favorable predictive biomarker for response to alkylating chemotherapy, either alone or in combination with radiotherapy. Several pieces of evidence show a correlation between the LOH in 1p/19q and the mutational status of IDH1/2, TERT, and MGMT. Here we used an approach to identify the codeletion in these chromosomal arms based on the analysis of 54 SNPs, which enabled the identification of LOH in six patients.
The present study, using a targeted sequencing method, enabled a comprehensive genetic profile of each sample tested. The diagnostic and therapeutic benefits for the use of this approach are numerous. The possibility to consider more effective therapies associated with molecular alterations represents one of the most important innovations in oncology. Several tumors are characterized by a high degree of heterogeneity, and their primitive molecular profile may also change in response to ongoing treatments. Compared to other forms of brain tumors, gliomas are the most heterogeneous in terms of molecular profile and prognosis [83]. Therefore, the diagnostic approach should include high-throughput methods, such as NGS, to investigate the presence of multiple molecular alterations [84]. In addition to the benefits for patients, the advantages of NGS also rely on the reduction of time to obtain results and reduced costs. The economic advantage of the NGS-based approach versus traditional diagnostic approaches (Sanger sequencing or Real-Time PCR) has been already discussed [85]. Targeted sequencing, in particular, may allow a significant reduction in costs and simplification of clinical interpretation of the identified genotypes. Systematic profiling of the comprehensive set of mutations related to glioma may offer new clinical opportunities for tailored treatments, shorter turnaround times, and less labor intensive and more affordable testing.

Conclusions
The goal of the genomic approach proposed here is to better characterize the molecular portrait of gliomas, enabling faster tumor classification and better genotype-phenotype characterization. This high-throughput approach may enable personalized medicine and contribute to a better healthcare with a cost reduction.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/life12070956/s1, Figure S1: Amplicons' coverage of chr1. The graph shows the 29 amplicons used to assess copy number alterations of chr1; Figure S2: Amplicons' coverage of chr19. The graph shows the 25 amplicons used to assess copy number alterations of chr19; Table S1: On-demand NGS panel designed using the custom IonAmpliSeq Designer software; Table S2: Designed BED file; Table S3: Summary of clinically relevant variants identified in our cohort; Table S4: Frequencies of genetic variants according TCGA database; Table S5: NGS summary report.

Conflicts of Interest:
The authors declare no conflict of interest.