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.

1. 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 of 60 years) [4,5]. Based on gene expression, the Cancer Genome Atlas (TCGA) classifies glioblastomas into proneural, neural, mesenchymal, and classical subtypes [6,7]. Secondary glioblastomas are essentially proneural, whereas primary glioblastomas belong to any of the indicated subtypes [8]. More specifically, in the oncogenesis of glioblastomas, biologically relevant alterations in three essential signaling pathways have been reported: RTK/RAS/PI3K (88%), p53 (87%), and retinoblastoma protein (78%) [9]. Other pathways include cell cycle, DNA damage responses, cell death and differentiation, neovascularization, and additional key players, such as Kras, G2M checkpoint, and Wnt [10,11].

The fifth edition of the WHO Classification of Tumors of the Central Nervous System (WHO CNS5), published in 2021, renamed glioblastoma as “glioblastoma, IDH-wild-type” (from now on referred to as GBM) and suggested supplementary guidelines for the histological classification of gliomas, indicating a spectrum of diagnostic, prognostic and/or predictive molecular biomarkers (Figure 1) [12]. With regard to genetic alterations, GBM is typically isocitrate dehydrogenase (IDH)-wild-type and has epidermal growth factor receptor (EGFR) amplification or mutations, phosphatase and tensin homolog (PTEN) deletion or mutations, and loss of heterozygosity on 10q [2]. Additional glioma subtypes may have recurrent point mutations in IDH1 and IDH2, platelet-derived growth factor receptor (PDGFR), cyclin-dependent kinase inhibitor 2A (CDKN2A) deletion, and tumor protein P53 (TP53) [13]. According to the WHO CNS5 classification, complemented by the Consortium to Inform Molecular and Practical Approaches to CNS Tumor Taxonomy (cIMPACT-NOW), IDH-mutated gliomas can be further subdivided into astrocytomas (grade 2, 3, or 4) and oligodendrogliomas with 1p/19q codeletion status (grade 2 or 3) [12,14,15] (Figure 1). Additional genetic markers are considered potentially relevant for tumor classification, such as the concurrent gain of whole chromosome 7 and the loss of whole chromosome 10, and mutations in the telomerase reverse transcriptase (TERT) promoter that are strongly correlated with glioblastoma IDH-wild-type, while mutations of B-Raf proto-oncogene (BRAF) and H3.3 histone A (H3F3A) associated with alterations of X linked intellectual disability (ATRX) and tumor protein P53 (TP53) are common in pediatric GBM [16,17,18]. Driver mutations (K27M and G34R/V) in H3 clustered histone 2 (HIST1H3B) are linked to pediatric and young high grade GBM with poor outcome [19].

Figure 1

Integrated histomolecular classification of diffuse gliomas according to the 2021 WHO Classification of Tumors of the Central Nervous System complemented by the Consortium to Inform Molecular and Practical Approaches to CNS Tumor Taxonomy (cIMPACT-NOW).

The combinatory use of both histopathological and molecular features by the WHO CNS5 is nowadays mandatory for the diagnostic classification of gliomas [12,20]. Molecular testing includes the evaluation of markers with predictive and/or prognostic value, such as the O6-methylguanine (O6-MeG)-DNA methyltransferase (MGMT) promoter methylation, the mutational status of IDH1/2, and the presence of chromosomes 1p/19q codeletion [21]. Typically, Karnofsky performance status (KPS) > 70, age at diagnosis < 50 years, MGMT methylation (>9%), gross/subtotal resection (>90%), and location of tumor in a non-eloquent area of the brain represent the most favorable prognostic indicators [22]. Loss of heterozygosity (LOH) of 1p/19q is another prognostic indicator of better prognosis, while MGMT promoter methylation status is the strongest prognostic factor for GBM outcome and a strong predictor of response to alkylating chemotherapy such as temozolomide (TMZ) [23]. Therefore, treatment strategies in patients with gliomas are based on tissue diagnosis and evaluation of molecular markers. The current standard protocol for treatment of all newly diagnosed gliomas is gross total resection, when feasible, followed by radio- and chemotherapy with TMZ. In addition, therapies with alkylating agents of the nitrosourea class, such as lomustine, carmustine, nimustine, or fotemustine have been used but without effective benefits [24]. Preclinical studies suggest that the combination of TMZ with bevacizumab, an anti-VEGF antibody approved in the United States, Canada, Switzerland, and several other countries outside the European Union for the treatment of recurrent GBM, temporally improves clinical outcomes and quality of life but does not impact overall survival (OS) [24]. However, despite current therapeutic treatments, GBM often inexorably recurs in the same location or in another area of the brain about 14–15 months after neurosurgery.

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 high-throughput 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].

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).

2. Materials and Methods

2.1. 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, age, MGMT methylation, IDH1/2 mutation status, 1p/19q LOH, and Ki-67 (MIB-1) labeling index. Methylation of MGMT promoter (GRCh37/hg19 chr10: 131.265.471-131.265.581) was assessed by pyrosequencing with the “Pyromark Q96 ID System”. The assessment of 1p/19q LOH, hotspot mutations in IDH1 (c.394C > T p.R132C, c.394C > G p.R132G, c.394C > A p.R132S, c.395G > A p.R132H, c.395G > T p.R132L, c.395G > C p.R132P), IDH2 (c.514A > C p.R172R, c.514A > G p.R172G, c.514A > T p.R172W, c.515G > A p.R172K, c.515G > C p.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].

2.2. 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).

2.3. 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.

Table 1

Target regions of the NGS panel.

Target	Chromosome	Number of Amplicons	Covered Bases	% Overall Coverage	Number of Exons
ACVR1	chr2	23	1710	100	9
ATRX	chrX	106	8161	99	35
BRAF	chr7	42	2655	99	18
CDKN2A	chr9	11	1008	99	5
EGFR	chr7	63	4489	100	30
H3F3A	chr1	5	273	58	3
HIST1H3B	chr6	5	431	100	1
HIST1H3C	chr6	5	431	100	1
IDH1	chr2	20	1405	100	8
IDH2	chr15	20	1572	99	11
TP53	chr17	13	1503	100	12
PDGFRA	chr4	50	3710	100	22
PTEN	chr10	22	1901	98	10
TERT-Promoter	chr5	1	124	99	n/a
29 SNPs	chr1	30	30	100	n/a
25 SNPs	chr19	24	24	100	n/a

2.4. 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.

2.5. 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 (somatic_low_stringency_540).

2.6. 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.

2.7. 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 non-informative 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.

3. 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.

3.1. 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].

Table 2

Clinical and molecular characteristics of patients and related tumor WHO grade and location.

Sample ID	WHO Grade	Tumor Location	Sex	Age	MGMT	Idh1/2 1	1p/19q LOH 1	TERT 1	Ki-67 (MIB-1)
ID 19	IV	Right TI	M	67	4%	wt	wt	-	25%
ID 48	IV	Right F	M	72	8%	wt	wt	-	20%
ID 78	IV	Right T	F	67	22%	wt	wt	-	17%
ID 93	IV	Right F	M	62	44%	wt	wt	-	35%
ID 41	IV	Right TI	F	45	23%	wt	wt	-	37%
ID 66	IV	Right T	M	47	63%	wt	wt	-	30%
ID 117	IV	Right FTP	F	56	67%	wt	wt	-	30%
ID 169	IV	Right F	M	52	13%	wt	wt	–146 C > T	67%
ID 121	IV	Right P	M	50	26%	wt	wt	-	60%
ID 143	IV	Right F	M	61	23%	wt	wt	–124 C > T	30%
ID 31	IV	Left FP	M	75	4%	wt	wt	-	25%
ID 30	IV	Left PO	F	76	41%	wt	wt	-	-
ID 21	IV	Right P	F	54	49%	wt	wt	-	80%
ID 164	IV	Right F	F	63	2%	wt	wt	-	40%
ID 40	IV	Right T	F	76	27%	wt	wt	-	35%
ID 166	IV	Left T	F	73	4%	wt	wt	-	-
ID 153	III	Right F	F	67	72%	p.R132H	LOH	-	4%
ID 209	IV	Left P	M	74	56%	p.R132H	LOH	-	70%
ID 232	II	Left F	F	48	53%	p.R132H	LOH	-	5%
ID 233	IV	Right T	M	42	60%	p.R132H	LOH	-	20%
ID 258	II	Left T	M	47	33%	p.R132H	LOH	-	4%
ID 260	IV	Left CC	M	41	5%	wt	LOH	-	-
ID 278	II	Left F	F	52	35%	p.R132H	LOH	-	1%
ID 49	IV	Left F	M	51	32%	wt	wt	-	50%
ID 79	IV	Right P	F	71	28%	wt	wt	-	20%
ID 38	IV	Right PO	M	54	5%	wt	wt	-	20%
ID 192	IV	Left T	F	91	-	wt	wt	-	-
ID 199	IV	Right F	M	53	-	wt	wt	-	-
ID 176	IV	Right F	M	62	3%	wt	wt	-	-
ID 104	IV	Right F	F	70	3%	wt	wt	-	85%
ID 09	IV	Left FP	F	81	4%	wt	1p	-	15%
ID 10	IV	Right FI	F	63	9%	wt	19q	-	35%
ID 85	IV	Left P	M	59	4%	wt	19q	-	30%
ID sc21	IV	Right P	M	46	41%	wt	1p	-	60%

T = temporal lobe; I = insular lobe; F = frontal lobe; P = parietal lobe; O = occipital lobe; - = test not performed; M = Male; F = Female. 1 The molecular profile was performed by the MassARRAY Analyzer 4 system (Agena Bioscience, CA, USA).

3.2. 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).

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).

3.3. 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.

Figure 3

Oncoplot showing the distribution of pathogenic variants detected by tNGS in tumor samples. Each column represents a sample and each row the genes. The upper barplot shows the number of mutations in each sample, while the right barplot shows the frequency and type of mutations in each gene. The central plot shows the gene mutations in each sample. Mutation types are indicated in the legend below. The detailed list of variants is reported in Supplementary Material S1, Table S3.

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.

Figure 4

Analysis of the 1p/19q codeletion detected by the assessment of SNP allelic frequencies. The presence of all markers in 1p and 19q in homozygous form and the presence of at least one marker in heterozygous form in the opposite arm indicate codeletion of these chromosomal regions. The blue line indicates the homozygous frequency; orange lines show the established ranges of heterozygosity.

3.4. 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.

4. 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.

5. 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.