Case of multifocal glioblastoma with four fusion transcripts of ALK, FGFR2, NTRK2, and NTRK3 genes stresses the need for tumor tissue multisampling for transcriptomic analysis.

Glioblastoma multiforme (GBM) is the most malignant brain tumor with patient mortality rate close to 100%, 5-yr survival rate of ∼5%, and a median survival of 14 mo. GBMs have notorious histomorphologic and molecular heterogeneities thus giving hope for development of future personalized therapies. We describe here a case of a 48-yr-old male patient with three-nodular GBM. To address the question of intratumoral molecular heterogeneity, a comparative analysis of gene expression was performed by using multiple samples collected from different tumor sites with the aid of intraoperative magnetic resonance imaging (MRI). Sixteen GBM biosamples from parietal, temporal, and temporo-polar localizations were collected from primary, recurrent, and second recurrent tumors and were obtained and investigated by RNA sequencing. Our investigations revealed that biosamples derived from different tumor sites differ in their gene expression profiles with classical or mesenchymal signatures associated with clinically distinct molecular subtypes of GBM found within the same tumor. The results also showed significant differences in the expression of genes specific for targeted therapeutics. Our investigations have enabled the identification of four novel fusion transcripts-KIF5C-NTRK3, AC016907.2-ALK, CNTNAP3-NTRK2, and ZNF135-FGFR2-each present in only one sample. We found no differences between untreated and recurrent stages in the expression levels of genes involved in fusion transcripts, suggesting the lack of association between fusion transcript and treatment response. In contrast, longitudinal changes in the expression of VEGF and MGMT genes were concordant with the tumor response to bevacizumab and temozolomide. Our study underscores the importance of integrating a multisampling approach and RNA sequencing and demonstrates the predictive merit of an integrated approach for differentiating genomic aberrations associated with untreated or post-treatment recurrent GBMs.

INTRODUCTION

Glioblastoma multiforme (GBM) is the most aggressive form of brain tumors with almost 100% lethal outcome and a 5-yr survival rate of ∼5% (Delgado-López and Corrales-García 2016). The current standard of care for GBM comprises surgery followed by combined radiochemotherapy (Stupp et al. 2009). Despite aggressive treatment regimens, a median survival of 14 mo remained unimproved since the introduction of the current standard of care a decade ago (Young et al. 2015). Disease recurrence is inevitable with the lack of effective therapies for recurrent GBMs posing a major challenge to improving clinical outcomes (Weller et al. 2013). The elevated level of intrinsic and acquired resistance to the standard therapy is the major cause of short survival in patients with GBM. Advances in understanding the molecular basis of GBM call for rethinking the prospects of the “one-treatment-fits-all” approach and necessitates new concepts for diagnostics, patient stratification, and treatment of GBMs (Kalasauskas et al. 2017).

Recent investigations revealed that characteristic genomic alterations and gene expression signatures that were previously thought to be associated with different molecular subtypes of GBM can exist within the same tumor (Sottoriva et al. 2013; Patel et al. 2014; Mandel et al. 2016; Kim et al. 2020). Sampling strategies, therefore, gain a particular importance in molecular diagnostics of GBM. Previous attempts were mostly based on the assumption that molecular evaluation of one tissue specimen per tumor is sufficient for molecular stratification of GBMs and predicting treatment response (Mandel et al. 2016). However, there is a growing realization that traditional diagnostics based on analysis of a single tumor biopsy is not adequate for GBMs (Sottoriva et al. 2013; Delgado-López and Corrales-García 2016; Morrissy et al. 2017; GLASS Consortium et al. 2018). Although it is not clear how many specimens from the same tumor might be needed (Buzdin et al. 2018) the degree of intratumoral heterogeneity varying across different tumor types seems to be the determining factor. Mathematical predictions suggest that on average at least five samples are required to compensate for the impact of molecular intratumoral heterogeneity and increase the level of confidence in molecularly based stratifications for GBMs but not for medulloblastomas (Morrissy et al. 2017).

To date, several molecular markers have been implemented in the routine GBM diagnostics. The World Health Organization (WHO) classification of central nervous system (CNS) tumors categorizes clinically distinct subgroups of GBM based on mutational status of the IDH1 and IDH2 genes coding for isoforms of NADP+-dependent isocitrate dehydrogenase; patients with mutated tumors have a better prognosis (Yan et al. 2009; Louis et al. 2016). We showed previously that elevated levels of FREM2 mRNA and protein are associated with favorable prognosis of IDH-WT GBMs (Jovčevska et al. 2019). A higher methylation status of MGMT gene promoter is linked with better prognosis of tumor response on chemotherapeutic temozolomide (Hegi et al. 2005; Wick et al. 2014; Dovek et al. 2019). Finally, gene expression signatures were identified that can classify GBMs into three clinically distinct subtypes: mesenchymal, proneural, and classical (Phillips et al. 2006; Wang et al. 2017). Recently it was shown that most mesenchymal GBMs arise as, and evolve from, a proneural-like precursor (Ozawa et al. 2014). Previous studies utilizing single-cell RNA sequencing pointed at significant interpatient heterogeneity of infiltrating (Darmanis et al. 2017), metabolic (Saurty-Seerunghen et al. 2019), and drug response (Reinartz et al. 2017) capacities between GBM cells, which was also associated with tumor prognosis. Furthermore, several clinically actionable fusion transcripts of ALK, ROS1, NTRK1-3, FGFR3, and BCR genes (Jones et al. 2019; Gilani et al. 2020; Nørøxe et al. 2020; Torre et al. 2020) were found in GBMs, and using the corresponding targeted drugs as the experimental therapy is currently debated (Jones et al. 2019). However, little is known about penetrance of these and other molecular markers in GBM considering tumor multisampling for individual patient.

We describe here a case of 48-yr-old male patient with multifocal GBM for whom neuronavigation and intraoperative magnetic resonance imaging (MRI)-guided multisampling of tumor tissues was performed to investigate patterns of gene expression in different tumor sites. The results showed markedly different localization-specific gene expression profiles, also significantly differing in expressions of genes specific for targeted therapeutics.

In the different tumor sites, we identified four new fusion transcripts of KIF5C-NTRK3, AC016907.2-ALK, CNTNAP3-NTRK2, and ZNF135-FGFR2 genes that potentially could be clinically actionable. However, we detected no increased expressions of the respective oncogenes, which suggests that these rearrangements could be just “passenger” results of the previous genotoxic therapies and heterogeneity of GBM cells. On the other hand, expression profiles of VEGF and MGMT genes were congruent with the response records of tumors to bevacizumab and temozolomide, respectively. Our results strongly suggest that multisampling should be applied, when possible, to molecularly characterize GBM and that RNA sequencing analysis can help differentiate clinically actionable from passenger mutations.

RESULTS

A 48-yr-old male with a history of thyroid carcinoma 12 yr earlier presented with generalized seizures. Thyroid stimulating hormone (TSH) level was 0.40 mIU/L (normal range, 0.38–5.33), preoperative positron emission tomography (PET)-computed tomography (CT) test is shown in Supplemental Figures S1and S2. The MRI in January 2018 revealed a multifocal contrast-enhancing right temporal lesion with dural and extradural extension including bone destruction that was visible on CT scan. Furthermore, there was a large right postcentral non contrast enhancing lesion hyperintense on T2 and fluid attenuated inversion recovery (FLAIR) sequences with a perifocal edema (Fig. 1; Supplemental Fig. S3). Because a thyroid carcinoma was operated on in 2006 including a neck dissection, a thorough staging was performed but revealed no recurrence of the thyroid tumor or systemic metastasis.

Figure 1.

Magnetic resonance imaging (MRI) follow-up of tumor treatment. Biopsy sites are pointed with red arrows. Green arrows indicate absence of tumor at the resection sites. (A) Temporo-polar localization. (B) Temporal localization. (C) Parietal localization.

In January 2018, another conventional MRI was performed that now showed two new small ring-enhancing space occupying lesions in the temporal lobe posterior to the initial lesions (Fig. 1). According to the recommendation of the interdisciplinary institutional tumor board the temporal lesions were operated using 5-aminolevulinic acid (5 ALA) and intraoperative MRI guidance. The temporal lesions could be removed completely, and although the temporal polar lesion did not show any 5 ALA uptake, the new ring-enhancing lesions including a perilesional rim clearly showed a 5 ALA–positive signal. Because the post central lesion showed heterogeneous MR characteristics, biopsies documented by neuronavigation and intraoperative MRI were taken accordingly from different areas of the tumor (Fig. 1). Among them, biopsies of the parietal lesion were obtained from various regions of the tumor that did not show any 5 ALA–positive signal. The final intraoperative MRI confirmed a total resection of the temporal lesions and documented the localization of the biopsy's sites. Postoperatively the patient had an unremarkable course with no new neurological deficit.

In temporal lesion, histology showed unusual dural and bone invasion. Neuropathological diagnostic has classified the tumor as IDH1-wild type GBM WHO grade IV without overexpression of p53 and no loss of nuclear ATRX level, with no MGMT promoter methylation, and with C250T TERT promoter mutation (Table 1). The temporal tumor had glomerular vessels with multilayered endothelial cells and round nuclei (Fig. 2A, middle). Staining of the temporo-polar lesion evidenced a malignant astrocytic tumor forming band-like necrosis with palisading tumor cells (Fig. 2A, top) and thrombotic vessels (Fig. 2A, top).

Table 1.

Immunohistochemical and molecular pathological parameters of different lesions

	Temporal	Temporo-polar	Parietal
p53 IHC	Weak+	Weak+	Weak+
IDH1/2	WT	n.d.	WT
1p/19q	n.d.	n.d.	n.d.
ATRX	Nuclear ATRX expression	n.d.	n.d.
TERT	C250T	n.d.	C250T
CDKN2A/B	n.d.	n.d.	n.d.
MGMT	Unmethylated	n.d.	n.d.
PDGFRa	WT	n.d.	n.d.

(IHC) Immunohistochemistry, (WT) wild type, (n.d.) not determined.

Figure 2.

(A) Hematoxylin and eosin staining of (top) temporo-polar, (middle) temporal, and (bottom) parietal glioblastoma multiforme (GBM); example of atypically configured mitosis is highlighted with a circle. (B) Distribution of the tumor in the cranium, including annotation of biopsy sites. (C) Sample collection dates and volume kinetics of tumor regions (fold change for sum of dimensions for target lesions) for different tumor localizations.

In the parietal lesion no solid tumor tissue was found, but there was a diffuse infiltration zone with increased density of tumor cells and residual elements of normal brain. Many of the tumor cells had pleomorphic nuclei and there were detectable sporadic mitoses (Fig. 2A, bottom). No IDH1 R132H mutation but the C250T TERT mutation identical to the temporal lesion were identified (Table 1).

Four weeks postoperatively, the patient started receiving 1 mo of bevacizumab treatment because of the rapid increase in cerebral blood volume (CBV) in the parietal lesion and thereafter a fractionated 60-Gy radiation therapy. A follow-up MRI in May 2018 showed some increase of the CBV for the parietal tumor with no contrast enhancement (Figs. 1 and 2C). Simultaneously, tumor progression/recurrence was observed in the temporo-polar localization (Figs. 1 and 2C). Based on the institutional tumor board recommendation gross-total resection of the tumor volume was performed in May 2018. The patient had an uneventful postoperative course.

In June 2018, the patient suffered from a cerebrospinal fluid fistula with secondary wound healing disorder and received surgical revision. Thereafter the patient received four cycles of temozolomide in a 5/23 regimen (340 mg for days 1–5, no temozolomide for days 6–28). MRI in October revealed progression of the tumor in both temporal and temporo-polar sites, with more rapid tumor increase in the temporal localization. In December 2018, the patient had reradiation (10 × 3.5 Gy) and started treatment with palbociclib, an inhibitor of cyclin-dependent kinases 4 and 6 (125 mg for 21 d every 28 d). Because of progressive hydrocephalus, the patient received the ventricular peritoneal shunt in January 2019. The patient passed away in March 2019.

The RNA sequencing analyses were performed for 10 biopsy samples obtained during the first surgery in January 2018 from temporal (n = 4), temporo-polar (n = 3), and parietal (n = 3) tumor sites, for three samples of recurrent and three samples of secondary recurrent parietal tumor. Recurrent tumor samples from parietal site were obtained during surgical interventions in May 2018 (R1) and January 2019 (R2). The molecular profiles obtained were compared with each other and with the reference human brain samples from other healthy individuals killed in road accidents obtained using the same RNA sequencing protocol (Suntsova et al. 2019).

Intratumoral heterogeneity of gene expression patterns was evident within the set of 10 primary biopsy samples. Different GBM samples fell into two major clusters, of which one corresponded mostly to temporal and temporo-polar, and another mostly to parietal localizations (except secondary recurrent tumors; Fig. 3A). Comparing the expression profiles of selected 162 cancer drug target genes (Law et al. 2014) resulted in a similar clustering pattern (Fig. 3B). We then looked at the GBM diagnostic 150-gene signature specific for classical, proneural, and mesenchymal subtypes according to Wang et al. (2017). A fraction of temporal (T1 and T2) and temporo-polar (TP2) tumor samples could be attributed to the mesenchymal GBM subtype associated with a more malignant phenotype, whereas all parietal samples and rest of temporal and temporo-polar samples belonged to the classical subtype (Fig. 3C) having a better prognosis (Wang et al. 2017). However, not all genes were informative for predicting mesenchymal and classical subtype according to unsupervised clustering of the 150 genes (Supplemental Fig. S4). To assess factors other than tumor heterogeneity that could explain such clustering pattern, we calculated tumor purity using the ESTIMATE bioinformatical tool (Yoshihara et al. 2013). Tumor purity estimates were in the range 0.6–0.94 across the studied biosamples (Supplemental Table S1). However, the clustering pattern could not be explained by the tumor purity: among the top five samples with the highest purity (>80%) (P3, R2–P2, TP3, P2, and P1), only P2 and P3 fall within one cluster, whereas other samples did not (Fig. 3A). Thus, tumor purity unlikely explains the observed clustering pattern. We further analyzed RNA-seq profiles using CIBERSORT tool to assess immune infiltration (Newman et al. 2015). No infiltration was detected in all samples except T1 and T2, where CIBERSORT indicated infiltration with macrophages (Supplemental Table S1). Interestingly, T1 and T2 samples were clustering together in all dendrograms on Figure 3. Thus, macrophage infiltration detected by CIBERSORT may at least partly explain T1 and T2 coclustering.

Figure 3.

Characteristic of temporal, temporo-polar, parietal, recurrent, and secondary recurrent parietal experimental glioblastoma multiforme (GBM) biosamples. Color code indicates sample type (parietal, recurrent, and secondary recurrent), parietal, temporal, temporo-polar samples of our patient, and healthy brain tissue controls. (A) Hierarchical clustering dendrogram built using expression levels of 36,596 gene transcripts. (B) Dendrogram built using expression of cancer drug target genes (Law et al. 2014) for the same biosamples. (C) Dendrogram built using diagnostic gene signature for proneural, mesenchymal, or classical GBM subtypes according to Wang et al. (2017) for the same experimental biosamples. The vertical color code indicates GBM subtype.

Oncobox tests (Poddubskaya et al. 2018, 2019; Buzdin et al. 2020) have been performed to simulate cancer drug efficiencies using molecular data from all 10 first-surgery specimens. The results obtained for bevacizumab and temozolomide were in line with the documented clinical records for drug responses of temporal/temporo-polar and parietal tumors (Table 2). In particular, temporo-polar (n = 3) samples had lower balanced efficiency scores (BESs) for bevacizumab than temporal (n = 4) samples (t-test P-value = 0.033) (Supplemental Table S2, sheet “BES”). Higher BES values indicate increased predicted drug efficiency. Concordant with predictions of treatment response based on BES values, the temporal tumor did not recur after bevacizumab treatment, whereas temporo-polar recurred (Fig. 2B). At the same time, parietal tumors that had BES values significantly lower than temporal (t-test P-value = 0.023; Supplemental Table S2, sheet “BES”) did not change in size during bevacizumab treatment (Fig. 2B). In turn, temporo-polar samples had higher BES for temozolomide than temporal samples (P = 0.0034), which was congruent with the onset of rapid regrowth of the temporal tumor upon temozolomide treatment. Temozolomide BES values were not significantly different in parietal versus temporal and parietal versus temporo-polar comparisons.

Table 2.

Characteristic of samples from different localizations

Sample code	Bevacizumab BES	Temozolomide BES	Tumor site	Fusion transcript
P1	6.27	−0.20	Parietal	-
P2	9.69	−0.08	Parietal	-
P3	9.42	−0.24	Parietal	ZNF135-FGFR2
T1	13.53	−0.32	Temporal	-
T2	13.48	−0.33	Temporal	-
T3	15.86	−0.27	Temporal	-
T4	10.96	−0.21	Temporal	AC016907.2-ALK
TP1	6.71	−0.09	Temporo-polar	-
TP2	10.64	−0.11	Temporo-polar	-
TP3	9.14	−0.13	Temporo-polar	-
R1–P1	n.d.	n.d.	Parietal, recurrent	KIF5C-NTRK3
R1–P2	n.d.	n.d.	Parietal, recurrent	-
R1–P3	n.d.	n.d.	Parietal, recurrent	-
R2–P1	n.d.	n.d.	Parietal, second recurrent	CNTNAP3-NTRK2
R2–P2	n.d.	n.d.	Parietal, second recurrent	-
R2–P3	n.d.	n.d.	Parietal, second recurrent	-

(BES) Balanced efficiency score, (n.d.) not determined.

In addition, gene expression analysis identified 695 differentially expressed genes (DEGs) between the temporal and temporo-polar localizations with a false discovery rate (FDR)-adjusted P-value of <0.05 (Supplemental Table S2). We then calculated DEGs for each tumor site versus normal samples. A Venn diagram illustrating the DEGs observed between the temporal, temporo-polar, and parietal sites is shown in Supplemental Figure S5: 426/235, 947/1062, and 620/987 up/down regulated genes were specific for initial parietal, temporal, and temporo-polar samples, respectively. The DEGs observed between the samples taken from the primary, first-recurrent, and second-recurrent parietal tumors are shown in Supplemental Figure S6: 290/322, 650/995, and 926/667 up/down regulated genes were specific for initial parietal, first-recurrence parietal, and second-recurrence parietal samples, respectively. A principal component analysis (PCA) plot for all studied biosamples is shown in Supplemental Figure S7.

Finally, in the different tumor sites we identified four new fusion transcripts of the KIF5C-NTRK3, AC016907.2-ALK, CNTNAP3-NTRK2, and ZNF135-FGFR2 genes (Table 3; Fig. 4). The presence of novel fusions was confirmed by reverse transcription polymerase chain reaction (RT-PCR) and Sanger sequencing (data not shown). Each of these fusion transcripts was identified in only one sample from either newly diagnosed or recurrent tumors (Table 2). Notably, none of the fusion transcripts identified in the untreated tumor could be found in recurrent samples and vice versa. Functional cancer fusion transcripts are known to enhance expressions of oncoproteins, where the upstream fusion partner serves as the driver (Gao et al. 2018). It was concluded previously that fusions drive the development of ∼17% of all cancer cases and function as the sole driver in ∼1% of them. Furthermore, ∼6% of tumors have druggable fusions (Gao et al. 2018).

Table 3.

Details for the four fusion transcripts found in the patient

Gene	Position 1	Position 2	Variant type	Predicted effect	ClinVar ID
KIF5C-NTRK3	Chr 2:148994538	Chr 15:87940753	Translocation	Fusion transcript	SCV001745862
AC016907.2-ALK	Chr 2:30119772	Chr 2:29717697	Translocation	Fusion transcript	SCV001745863
CNTNAP3-NTRK2	Chr 9:67201219	Chr 9:84751986	Translocation	Fusion transcript	SCV001745864
ZNF135-FGFR2	Chr 19:58061706	Chr 10:121551459	Translocation	Fusion transcript	SCV001745865

Figure 4.

Schematic representation of the fusion transcripts identified for the KIF5C-NTRK3 (A), AC016907.2-ALK (B), CNTNAP3-NTRK2 (C), and ZNF135-FGFR2 (D) chimeric RNAs. (Left panels) Gene structures upstream and downstream from fusion site. (Right panels) Gene expression levels detected for the fusion partners in all experimental biosamples. Color code indicates the current biosample with a fusion (red), other tumor biosamples (gray), and the control healthy brain tissue samples (green).

To estimate whether the fusion transcripts discovered could be functional, we compared expressions of the fusion partner genes among all the samples investigated (Fig. 4). Clinically relevant/druggable oncogenic gene fusions are very frequently accompanied by overexpression of the corresponding protein: Fluorescence in situ hybridization (FISH) is highly concordant with immunohistochemistry (IHC) for the ALK (Lin et al. 2019), ROS1 (Sholl et al. 2013), and NTRK (Hechtman et al. 2017) genes. We found that the presence of a fusion did not coincide with increased expression of the respective oncogenic fusion partner gene in all the cases (Fig. 4). Thus, this information was not used to set an experimental therapy for the current patient.

DISCUSSION AND CONCLUSION

We describe here a case of 48-yr-old male patient with a multifocal GBM (temporal/temporo-polar and parietal) for whom neuronavigation and intraoperative MRI-guided multisampling of tumor tissues was used to investigate patterns of gene expression in different tumor sites. The results showed markedly different localization-specific gene expression profiles, also significantly differing in expression of genes specific for targeted therapeutics (Supplemental Table 2). The biosamples from different localizations also related to two different molecular phenotypic groups (Wang et al. 2017). In addition, the temporal/temporo-polar and parietal parts of the patient's tumor had statistically significantly different simulated drug efficiency scores (BESs) for cancer therapeutics used in this patient (bevacizumab, temozolomide). These apparent phenotypic differences based on gene expression also reflected the patient's tumor response to treatment. These findings are in line with the results previously documented on the cohort of 128 primary and recurrent tumor samples, suggesting that a transcriptomics-based BES value can serve as the high-quality predictor of the temozolomide efficacy in GBM (Kim et al. 2020).

In four out of 26 tumor samples, we identified four previously unpublished fusion transcripts: KIF5C-NTRK3, CNTNAP3-NTRK2, AC016907.2-ALK, and ZNF135-FGFR2, one fusion per sample (Table 2; Fig. 4). Two fusions (AC016907.2-ALK, ZNF135-FGFR2) were found in the primary tumor biopsy specimens (one in temporal localization and one in parietal) and thus could represent the notable heterogeneity of a primary tumor. The other two fusions (KIF5C-NTRK3 and CNTNAP3-NTRK2) were found in the recurrent parietal tumors and thus could appear because of the mutagenic effect of radiation therapy and chemotherapy treatments (Bao et al. 2014).

That none of the four fusions coincided with an increased expression of the respective oncogenes (NTRK2-3, ALK, and FGFR2) in the same biosamples (Fig. 4) suggests that they are likely not the drivers of molecular evolution of the current tumor, but rather nonfunctional passenger mutations. Lack of overlap between fusion transcripts found in different tumor samples is concordant with such an explanation. This suggests that the presence of a fusion transcript with known oncogene per se should not be regarded as the sufficient condition for experimental therapy; instead, a gene expression assay should be carried out to confirm if a fusion significance is supported by an enhanced oncogene expression compared to the reference samples or not. On the other hand, our results do not either rule out the possibility that newly diagnosed and recurrent GBMs may differ in the spectrum of fusion transcripts or that some of the fusions identified in recurrent tumors could appear as a result of treatment-driven increase of genomic instability. Such an interpretation would be concordant with the previously found association between temozolomide treatment and increased rates of genomic aberrations/mutations in recurrent GBMs (TCGA study). Even though our study involved the samples from only one patient, the results warrant further exploration of the fusion-transcriptomic profiles in newly diagnosed and recurrent GBMs.

Thus, we conclude that in-depth RNA sequencing analyses of multisampled tumor specimens can be beneficial for personalizing current and emerging glioblastoma treatments.

METHODS

The patient provided informed written consent for molecular analyses of his cancer biosamples and for presentation of relevant clinical and molecular data and their interpretation in this paper. The 10 primary tumor tissue samples used for gene expression analyses were obtained in January 2018 during surgery from different sites of the tumor as documented by neuronavigation and intraoperative MRI (Supplemental Table 1). Recurrent tumor samples from parietal site were obtained during surgical interventions in May 2018 (R1) and January 2019 (R2). Volume kinetics of tumor regions was measured as fold change for sum of dimensions for target lesions. The samples were stored in RNA later in liquid nitrogen and then profiled by RNA sequencing (Illumina HiSeq 3000) at UCLA Center for Genomics and Bioinformatics according to Suntsova et al. (2019) following ribosomal RNA depletion protocol. The number of reads produced per library was 23.2–75.7 (average 40.8) million, and the number of gene-mapped reads was 11.7–25.2 (average 16.3) million. RNA sequencing FASTQ files were processed with STAR aligner (Dobin et al. 2013) in “GeneCounts” mode with the Ensembl human transcriptome annotation (Build version GRCh38 and transcript annotation GRCh38.89). Ensembl gene IDs were converted to HGNC gene symbols using Complete HGNC data set (https://www.genenames.org, database version of July 13, 2017). In total, expression levels were established for 36,596 annotated genes with corresponding HGNC identifiers.

The expression data of 16 biopsy specimens were deposited in the NCBI Sequencing Read Archive (SRA) repository with ID PRJNA590641. Differential gene expression analysis was done using DESeq2 (Love et al. 2014). The molecular pathway activation analyses and ranking of targeted therapeutics were performed using Oncobox cancer bioinformatic platform according to Poddubskaya et al. (2018, 2019), Borisov et al. (2020), Buzdin et al. (2020), and Tkachev et al. (2020).

Fusion transcripts were initially screened using STAR-Fusion software (Haas et al. 2019). Preliminary files containing fusion candidates were generated and the corresponding RNA sequencing reads were extracted. The output data were manually curated using UCSC BLAT and UCSC Browser (https://genome.ucsc.edu/) to interrogate fusion candidates according to the following criteria: (i) does the read cover exon junction of two different transcripts, (ii) if the junction point exactly matches exon termini of known genes with canonic splice sites, and (iii) if both transcripts are in the same orientation. The triple-positive candidates were identified as fusion transcripts and confirmed by RT-PCR with Sanger sequencing and investigated by RT-PCR in all tumor biosamples.

ADDITIONAL INFORMATION

Data Deposition and Access

The expression data of 16 biopsy specimens were deposited in the NCBI Sequencing Read Archive (SRA) repository with ID PRJNA590641. The variants were submitted to ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/) and can be found under accession numbers SCV001745862–SCV001745865.

Ethics Statement

The study was approved by the interdisciplinary institutional tumor board of the International Neuroscience Institute.

Acknowledgments

We thank International Neuroscience Institute–OmicsWay research initiative for clinical and technical support and OmicsWay Corp. research program for machine learning and digital oncology for software and pathway databases. Cloud-based computational facilities for this study were supported by Amazon and Microsoft Azure grants.

Author Contributions

A.S., S.K., C.H., and A.Gi. collected and interpreted patient data. A.S., C.H., and A.Gi. were involved in clinical management. C.H., S.K., and E.K. performed molecular analyses. A.S., L.M., and A.Gi. analyzed and interpreted the MRI data. M.S., A.Ga., A.M., and A.B. analyzed and interpreted gene expression data. A.S., S.K., M.S., A.M., A.B., and A.Gi. wrote the paper.

Funding

Financial support was provided by the ERA.Net RUS Plus program (GOTICS project), Russian Foundation for Basic Research grant 20-515-76007.

Competing Interest Statement

The authors have declared no competing interest.

Referees

Toni Rose Jue

Roel Verhaak