Publisher Correction: Genomic analyses of early responses to radiation in glioblastoma reveal new alterations at transcription, splicing, and translation levels

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

Glioblastoma is the most common intracranial malignant brain tumor with an aggressive clinical course. Standard of care entails maximally safe resection followed by radiotherapy with concomitant and adjuvant temozolomide as described in the landmark European Organization for Research and Treatment of Cancer (EORTC) Brain Tumor and Radiotherapy Groups and the National Cancer Institute of Canada study 1 . Adjuvant radiation therapy is commonly delivered 4-6 weeks after maximally safe resection.. Nonetheless, the median overall survival remains approximately 16 months 1,2 , and the recent addition of tumor-treating fields to the standard of care has only increased median overall survival to 20.5 months 1 . Recurrence occurs in part because glioblastoma uses sophisticated cellular mechanisms to repair DNA damage from double-stranded breaks caused by ionizing radiation, specifically homologous recombination and non-homologous end-joining. Thus, the repair machinery confers a mechanism for resistance to radiation therapy. Ionizing radiation can also cause base damage and single-strand breaks, which are repaired by base excision and single-strand break repair mechanisms, respectively 3 . A comprehensive analysis of molecular mechanisms driving resistance to chemotherapy and radiation is required to surpass major barriers and advance treatments for glioblastoma.
Changes in global transcriptome profile in response to radiation. We first conducted an integrated analysis to evaluate the early impact of radiation on the expression profile of U251 and U343 GBM lines. Using this approach, we focus on genes undergoing up regulation or downregulation post irraditaion after adjusting for cell line differences (See Methods). A relatively small number of genes displayed altered expression at T1 (Supplementary Fig. S4B,C). Downregulated genes are mainly involved in transcription regulation and include 18 zinc finger transcription factors displaying high expression correlation in glioblastoma samples from TCGA (Supplementary Table S1). Upregulated sets contain genes implicated in cell cycle arrest, apoptosis, and stress such as ZFP36, FBXW7, SMAD7, BTG2, and PLK3 (Supplementary Table S1).
Since many alterations were observed when comparing the T1 vs. T24 time points (Supplementary Table S1 and Supplementary Fig. S4C), we opted to focus on genes showing the most marked changes (log 2 fold-change >1.0 or <−1.0 and adjusted p-value < 0.05) to identify biological processes and pathways most affected at the T24 time point. Top enriched GO terms and pathways among downregulated genes include chromatin remodeling, cell cycle, DNA replication, and repair (Fig. 1A). Additionally, we identified several GO terms associated with mRNA metabolism, decay, translation, and ncRNA processing, suggesting active participation of RNA-mediated processes in radio-response (Fig. 1B). Network analysis indicated the set of genes in these categories is highly interconnected (Fig. 1C and Supplementary Table S2).
To expand the expression analysis, we employed Weighted Gene Co-expression Network Analysis (WGCNA) 22 to identify gene modules with significant co-expression variation in response to radiation. All identified modules, along with the complete list of genes in each module, are shown in Supplementary Fig. S5A Supplementary Fig. S5B). Among modules with the highest significant correlation (0.8, p-value < 1e−7), module 2 contains genes downregulated in T24, with many involved in cell cycle, metabolism mRNA metabolism, processing, splicing, and transport (Supplementary Table S3), corroborating results described above.
Next, we investigated downregulated genes with the gene set enrichment analysis (GSEA) tool Enrichr 23 and conducted expression correlation analysis with Gliovis 24 . Based on their genomic binding profiles and effect of gene expression, FOXM1 and the E2F family of transcription factors emerged as potential regulators of a large group of cell cycle/DNA replication-related genes in the affected set ( Fig. 2A, Supplementary Table S4). In agreement, E2F1, E2F2, E2F8, and FOXM1 displayed a significant decrease upon radiation. FOXM1 and E2F factors have been previously implicated in chromatin remodeling, cell cycle regulation, DNA repair, and radio-resistance 25,26 . All four factors are highly expressed in glioblastoma with respect to low-grade glioma.
www.nature.com/scientificreports www.nature.com/scientificreports/ Importantly, they display high expression correlation with a large set of downregulated genes implicated in cell cycle and DNA replication and among themselves in glioblastoma samples in TCGA (Fig. 2B,C). We corroborated the changes in expression of these transcription factors and some of their potential targets by qRT-PCR in U343 cells and also observed similar changes in the glioblastoma line A172 and the glioma stem cell line 3565 ( Supplementary Fig. S6).
Upregulated genes at T24 are preferentially associated with the extracellular matrix receptor interaction pathway, extracellular matrix organization, axonogenesis, and response to type I interferon (Fig. 3A, and  Supplementary Table S2). With respect to the extracellular matrix, we observed changes in the expression levels of several collagens (types II, IV, V, and XI), glycoproteins of the laminin family (subunits α β γ , , and ), and also integrins (subunits α, and β) (Fig. 3B, and Supplementary Table S1). Collagen type IV is highly expressed in glioblastoma and implicated in tumor progression 27 . In addition, it has been observed that the activation of two integrins, ITGB3 and ITGB5, contributes to radio-resistance 28 .
Radiation treatment also induced the expression of genes involved in neuronal differentiation and axonogenesis. Some key genes in these categories include SRC, VEGFA, EPHA4, DLG4, MAPK3, BMP4, and several semaphorins. These genes can have very different effects on glioblastoma development, with some factors activating oncogenic programs and others behaving as tumor suppressors. Similarly, type I interferon's effects on treatment are varied. For instance, interferon inhibited proliferation of glioma stem cells and their sphere-forming capacity and induced STAT3 activation 29 . On the other hand, chronic activation of type I IFN signaling has been linked to adaptive resistance to therapy in many tumor types 30 .
Activation of oncogenic signals post-radiation could counteract treatment effects and later contribute to relapse. We searched the set of highly up-regulated genes post-radiation for previously identified radio-resistance genes in glioblastoma, oncogenic factors and genes whose high expression is associated with poor prognosis (Supplementary Table S5). In Table 1, we list these genes according to their molecular function. Since several of these genes have never been characterized in the context of glioblastoma, our results open new opportunities to prevent radio-resistance and increase treatment efficiency. Importantly, there are inhibitors available against several of these proteins (Supplementary Table S4).
Changes in lncRNA profile in response to radiation. lncRNAs have been implicated in the progression of glioblastoma 31 , but their role in response to ionizing radiation is still poorly understood. We identified 161 lncRNAs with expression alterations in T1 vs. T24 comparisons. Analysis of this set with LnC2Cancer 32 identified several lncRNAs aberrantly expressed in cancer and with relevance to prognosis (Supplementary Table S1). We also detected significant downregulation of MIR155HG, whose high expression is associated with glioma www.nature.com/scientificreports www.nature.com/scientificreports/ progression and poor survival 33 . Another downregulated lncRNA with relevance to prognosis is linc000152, whose increased expression has been observed in multiple tumor types 34 . On the other hand, we observed a significant upregulation of two "oncogenic" lncRNAs, NEAT1 and FTX. NEAT1 is associated with tumor growth, grade, and recurrence rate in gliomas 35 , while FTX promotes cell proliferation and invasion through negatively regulating miR-342-3p 36 . Thus, if further studies corroborate NEAT1 and FTX as players in radio-resistance, targeting these lncRNAs should be considered to improve treatment response.  www.nature.com/scientificreports www.nature.com/scientificreports/ Effect of radiation on splicing. Alternative splicing impacts genes implicated in all hallmarks of cancer 37 and is an important component of changes in expression triggered by ionizing radiation 38 . All types of splicing events (exon skipping, alternative donor, and acceptor splice sites, multiple exclusive exons, and intron retention) were affected similarly upon exposure to radiation (Supplementary Table S6). At T24, we observed that transcripts associated with RNA-related functions (especially translation), showed the most splicing alterations. Affected transcripts encode ribosomal proteins, translation initiation factors, regulators of translation, and genes involved in tRNA processing and endoplasmic reticulum. Other enriched GO terms include mRNA and ncRNA processing, mRNA degradation, and modification. Catabolism is another process associated with several enriched terms, suggesting that splicing alterations in genes involved in catabolic routes could ultimately contribute to apoptosis (Fig. 4A,B, and Supplementary Table S7). Changes in the splicing profile are likely driven by an alteration in the expression of splicing regulators. In Table 2, we show a list of splicing factors displaying strong expression alterations. Among those previously connected to glioblastoma development, LGALS3 is the most extensively characterized.
LGALS3 is a galactosidase-binding lectin and non-classic RNA binding protein implicated in pre-mRNA splicing and regulation of proliferation, adhesion, and apoptosis; LGALS3 also is a marker of the early stage of glioma 39 . Differential translational efficiency. We used Ribo-seq 40 to identify changes in translation efficiency triggered by radiation (Supplementary Table S8). Translation, protein localization, and metabolism appear as top enriched terms among downregulated genes in T1 vs. T24 comparisons (Supplementary Table S9). In particular, several ribosomal proteins, along with translation initiation factors and mTOR, showed a significant decrease in translation efficiency (Fig. 5A,B). Overall, these results indicate repression of the translation machinery post-radiation exposure and its strong auto-regulation. Since changes in components of the translation machinery are occurring at all levels (transcription, splicing, and translation) at T24, we expect that major translational alterations take place in later stages of post-radiation.
In the upregulated set, we highlight three genes FTH1, APIP, and LRIG2 that could potentially counteract the impact of radiation (Supplementary Table S10). FTH1 encodes the heavy subunit of ferritin, an essential component of iron homeostasis 41 . Pang et al. 42 , reported that H-ferritin plays an important role in radio-resistance in glioblastoma by reducing oxidative stress and activating DNA repair mechanisms. The depletion of ferritin causes down-regulation of ATM, leading to increased DNA sensitivity towards radiation. APIP is involved in the methionine salvage pathway and has a key role in various cell death processes. It can inhibit mitochondria-mediated apoptosis by directly binding to APAF-1 43 . LRIG2 is a member of the leucine-rich and immunoglobulin-like domain family 44 , and its expression levels are positively correlated with the glioma grade and poor survival. LRIG2 promotes proliferation and inhibits apoptosis of glioblastoma cells through activation of EGFR and PI3K/ Akt pathway 45 . crosstalk between regulatory processes. Parallel analyses of transcription, splicing, and translation alterations in the early response to radiation provided an opportunity to identify crosstalk between different regulatory processes. The datasets showed little overlap, with just a few genes showing alterations in two different regulatory processes. However, we identified several shared GO terms when comparing the results of alternative splicing, mRNA levels (transcription), and translation efficiency (Supplementary Table S10). These terms show two main groups of biological processes. The first group indicates that the expression of genes involved in DNA and RNA synthesis and metabolism is particularly compromised. The second group is related to translation initiation. Ribosomal proteins were particularly affected (Figs. 4 and 5). There is growing support for the concept of specialized ribosomes. According to this model, variations in the composition of the ribosome due to the presence or absence of certain ribosomal proteins or alternative isoforms could ultimately dictate which mRNAs get preferentially translated 46 . Therefore, these alterations could later lead to translation changes of a specific set of genes.

Discussion
We performed the first integrated analysis to define global changes associated with the early response to radiation in glioblastoma. Our approach allowed the identification of "conserved" alterations at the transcription, splicing and translation levels and defined possible crosstalk between different regulatory processes. Alterations at the level of transcription were dominant, but changes affecting genes implicated in RNA mediated regulation were ubiquitous; they indicate that these processes are important components in radio-response and suggest that more robust changes in splicing and translation might take place later. E2F1, E2F2, E2F8, and FOXM1 as major drivers of transcriptional responses upon radiation. We observed marked changes in the mRNA levels of genes implicated in cell cycle, DNA replication, and repair 24 hours (T24) after radiation. Downregulation of several transcription factors, most of them members of the zinc finger family, was observed at one hour post-radiation. This group displays high correlations in expression within glioblastoma samples from TCGA, suggesting that they might work together to regulate gene expression. Unfortunately, most are poorly characterized, and the lack of information has prevented establishing further connections to changes in the cell cycle and DNA replication that we observed at T24.
GSEA and expression correlation analysis suggested that the downregulation of members of the E2F family is likely responsible for several of the expression changes we observed at T24. E2Fs have been defined as major   www.nature.com/scientificreports www.nature.com/scientificreports/ transcriptional regulators of the cell cycle. The family has eight members that could act as activators or repressors depending on the context, and are known to regulate one another. They are upregulated in many tumors due to overexpression of cyclin-dependent kinases (CDKs), inactivation of CDK inhibitors, or RB Transcriptional Corepressor 1 (RB1) and are linked to poor prognosis. Alterations in E2F genes can induce cancer in mice 47,48 . Specifically, we found that three E2F members showed decreased expression upon radiation: E2F1, E2F2, and E2F8, all of which have been previously implicated in glioblastoma development.
E2F1 is probably the best-characterized member of the E2F family. Besides its known effect on cell cycle regulation and DNA replication, it is also a positive regulator of telomerase activity, binding the TERT promoter 49 . Recent studies show that lncRNAs and miRNAs function in an antagonistic fashion to regulate E2F1 expression, ultimately affecting cell proliferation, glioblastoma growth, and response to therapy [50][51][52] . E2F2 has been linked to the maintenance of glioma stem cell phenotypes and cell transformation 53,54 . Several tumor suppressor miRNAs (let7b, miR-125b, miR-218, and miR-138) decrease the proliferation and growth of glioblastoma cells by targeting E2F2 53,55-57 . Although still poorly characterized in the context of glioblastoma, E2F8 drives an oncogenic phenotype in glioblastoma. Its expression is modulated by HOXD-AS1, which serves as a sponge and prevents the binding of miR-130a to E2F8 transcripts 58 . FOXM1 is another potential regulator of the group of cell cycle and DNA replication genes affected by radiation. FOXM1 is established as an important player in chemo-and radio-resistance and a contributor to glioma stem cell phenotypes 5,6,59-64 . FOXM1 and E2F protein have a close relationship and share target genes 65 . Additionally, FOXM1-and E2F2-mediated cell cycle transitions are implicated in the malignant progression of IDH1 mutant glioma 66 .
E2F and FOXM1 targeting could be considered as an option to increase radio-sensitivity. Since the development of transcription factor inhibitors is very challenging, an alternative to be considered is the use of BET (bromodomain and external) inhibitors. BET is a family of proteins that function as readers for histone acetylation and modulates the transcription of oncogenic programs 67 . Recent studies in glioblastoma with a new BET inhibitor, dBET6, showed promising results and established that its effect on cancer phenotypes comes via disruption of the transcriptional program regulated by E2F1 68 .
RnA processing and regulation as novel categories in radio-response. Besides the expected changes in expression of cell cycle, DNA replication and repair genes, radiation affected preferentially the expression of genes implicated in RNA processing and regulation. Additionally, we identified a co-expression module containing multiple genes associated with translation initiation, rRNA and snoRNA processing, RNA localization, and ribonucleoprotein complex biogenesis.
Many regulators of RNA processing are implicated in glioblastoma development, and splicing alterations affect all hallmarks of cancer 69 . Radiation-induced changes in the splicing patterns of oncogenic factors and tumor suppressors such as CDH11, CHN1, CIC, EIF4A2, FGFR1, HNRNPA2B1, MDM2, NCOA1, NUMA1, RPL22, SRSF3, TPM3, APC, CBLB, FAS, PTCH1, and SETD2. We also observed changes in expression of four RNA processing regulators previously identified in genomic/functional screening for RNA binding proteins contributing to glioblastoma phenotypes: MAGOH, PPIH, ALYREF, and SNRPE 70 . potential new targets to increase radio-sensitivity and prevent relapse. Activation of oncogenic signals is an undesirable effect of radiation that could influence treatment response and contribute to relapse. We observed increased expression or translation and splicing alterations of a number of pro-oncogenic factors, genes whose high expression is associated with poor survival and genes previously implicated in radio-resistance. www.nature.com/scientificreports www.nature.com/scientificreports/ Among genes with the most marked increase in expression upon radiation, we identified members of the Notch pathway (HES2, NOTCH3, MFNG, and JAG2). Notch activation has been linked to radio-resistance in glioblastoma, and Notch targeting improves the results of radiation treatment 71,72 . We also identified several genes associated with the PI3K-Akt, Ras, and Rap1 signaling pathways that increased expression levels upon radiation exposure. Targeting these pathways has been explored as a therapeutic option in glioblastoma 71,73 . Other oncogenic factors relevant to glioblastoma that had increased expression after radiation exposure include SRC, MUC1, LMO2, PML, PDGFR β, BCL3, and BCL6.
Anti-apoptotic genes (BCL6, RRM2B, and IDO1) also showed increased expression upon radiation. BCL6 is a member of the ZBTB family of transcription factors, which functions as a p53 pathway repressor. The blockage of the interaction between BCL6 and its cofactors has been established as a novel therapeutic route to treat glioblastoma 74 . RRM2B is an enzyme essential for DNA synthesis and participates in DNA repair, cell cycle arrest, and mitochondrial homeostasis. The depletion of RRM2B resulted in ADR-induced apoptosis, growth inhibition, and enhanced sensitivity to chemo-and radiotherapy 75 . IDO1 is a rate-limiting metabolic enzyme involved in tryptophan metabolism that is highly expressed in numerous tumor types 76 . The combination of radiation therapy and IDO1 inhibition enhanced therapeutic response 77 .
Among genes whose high expression correlates with decreased survival in glioblastoma, we identified several components of the "matrisome" and associated factors (FAM20C, SEMA3F, ADAMTSL4, ADAMTS14, SERPINA5, and CRELD1). The core of the "matrisome" contains ECM proteins, while associated proteins include ECM-modifying enzymes and ECM-binding growth factors. This complex of proteins assembles and modifies extracellular matrices, contributing to cell survival, proliferation, differentiation, morphology, and migration 78 . In addition, several genes of the proteinase inhibitor SERPIN family (SERPINA3, SERPINA12, SERPINA5, and SERPINI1) implicated in ECM regulation 79 were among those with high levels of expression upon radiation.
In conclusion, our results generated a list of candidates for combination therapy. Contracting the effect of oncogenic factors and genes linked to poor survival could increase radio-sensitivity and treatment efficiency. Importantly, there are known inhibitors against several of these proteins (Supplementary Table S5). Moreover, RNA processing and translation were determined to be important components of radio-response. These additional vulnerable points could be explored in therapy, as many inhibitors against components of the RNA processing and translation machinery have been identified 80,81 . Methods cell culture and radiation treatment. A172 was obtained from ATCC. U251 and U343 cells were obtained from the University of Uppsala (Sweden) and maintained in Dulbecco's Modified Eagle Medium (DMEM, Hyclone) supplemented with 10% fetal bovine serum, 1% Penicillin/Streptomycin at 37 °C in 5% CO 2humidified incubators and were sub-cultured twice a week. Cells were plated after appropriate dilution, and ionizing radiation treatment was performed on the next day at a dose of 5 Gray (Gy). 5 Gy was chosen since in current hypofractionated regimen, radiation therapy is delivered over the course of 5 fractions at 5 Gy per fraction. A cabinet X-ray system (CP-160 Cabinet X-Radiator; Faxitron X-Ray Corp., Tucson, AZ) was used. Glioma stem line 3565 was established in Dr. Jeremy Rich;s lab; cells were cultured in Neurosphere Media consisting of Neurobasal-A supplemented with B-27, 50 ng/ml EGF, and 50 ng/ml hFGF. qRt-pcR analysis. After exposure to ionizing radiation, cells were cultured for 14, 24, and 48 hours (hrs).
Total RNA was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. Reverse transcription was performed using High-Capacity cDNA reverse transcription kit (Applied Biosystems). Quantitative PCR was performed using TaqMan Universal PCR Master Mix (Applied Biosystems) or PowerUp SYBR Green Master Mix (Applied Biosystems) and reactions were performed on ViiA TM 7 Real-Time PCR System (Applied Biosystems). Data were acquired using ViiA 7 RUO software (Applied Biosystems) and analyzed using the 2-ΔΔ CT method with GAPDH as an endogenous control. Probes and primers used in qRT-PCR are listed in Supplementary Table S11. RNA preparation, RNA-seq and ribosome profiling (Ribo-seq). RNA was purified using a GeneJet RNA kit from Thermo Scientific. The TruSeq Ribo Profile (Mammalian) kit from Illumina was used to prepare material for ribosome profiling (Ribo-seq). RNA-seq and Ribo-seq samples were prepared according to Illumina protocols and sequenced at UTHSCSA Genome Sequencing Facility. overall strategy to identify gene expression alterations upon radiation. To identify the most relevant expression alterations in the early response to radiation, we analyzed samples from U251 and U343 cells collected at 0 (T0), 1 (T1), and 24 (T24) hours post-radiation. The time of one hour and 24 hours were chosen as this reflects the immediacy of the DNA damage response and resolution of double stranded DNA breaks. The production of double stranded DNA breaks occurs nearly immediately with the subsequent resolution of the breaks nearly immediately (within 15-60 minutes after exposure to ionizing radiation) with near resolution of all breaks by 24 hours. To capture the progressive dynamics of expression alterations, we compared T0 to T1 samples and T1 to T24 samples. Our strategy to identify the most relevant alterations in expression with maximal statistical power was to combine all samples and use a design matrix with cell type defined as a covariate with time points (Supplementary Fig. S1).
Sequence data pre-processing and mapping. The quality of raw sequences reads from RNA-Seq and Ribo-Seq datasets were assessed using FastQC 82 . Adaptor sequences and low-quality score (phred quality score <5) bases were trimmed from RNA-Seq and Ribo-Seq datasets with TrimGalore (v0.4.3) 83 . The trimmed reads were then aligned to the human reference genome sequence (Ensembl GRCh38.p7) using STAR aligner www.nature.com/scientificreports www.nature.com/scientificreports/ (v.2.5.2b) 84 with GENCODE 85 v25 as a guided reference annotation, allowing a mismatch of at most two positions. All the reads mapping to rRNA and tRNA sequences were filtered out before downstream analysis. Most reads in the Ribo-seq samples mapped to the CDS. The periodicity analysis was performed using ribotricer 86 . The number of reads assigned to annotated genes included in the reference genome was obtained by htseq-count 87 . Differential gene expression analysis. For differential expression analysis, we performed counting over exons for the RNA-seq samples. For translational efficiency analyses, counting was restricted to the CDS. Differential gene expression analysis was performed by employing the DESeq2 package 88 , with read counts from both U251 and U343 cell samples as inputs. Both the cell line and time were treated as covariates along with their interaction term. To find differentially expressed genes that show changes between two time points, we used time specific contrasts. The p-values were adjusted for controlling for the false discovery rate (adjusted p-value) using the Benjamini and Hochberg (BH) procedure 89 . Differentially expressed genes were defined with an adjusted p-value < 0.05.
Weighted gene co-expression network analysis. WGCNA 22 uses pairwise correlations on expression values to identify genes significantly co-expressed across samples. We used this approach to identify gene modules with significant co-expression variations as an effect of radiation. The entire set of expressed genes, defined here as those with one or higher transcripts per million higher (TPM), followed by variance stabilization) from U251 and U343 samples were clustered separately using the signed network strategy. We used the Z summary 90 statistic as a measure of calculating the degree of module preservation between U251 and U343 cells. Z summary is a composite statistic defined as the average of the density and connectivity based statistic. Thus, both density and connectivity are considered for defining the preservation of a module. Modules with > Z 5 summary were considered as significantly preserved. The expression profile of all genes in each co-expression module can be summarized as one "eigengene". We used the eigengene-based connectivity (kME) defined as the correlation of a gene with the corresponding module eigengene to assess the connectivity of genes in a module. The intramodular hub genes were then defined as genes with the highest module membership values (kME > 0.9). All analysis was performed using the R package WGCNA. The protein-coding hub genes were then selected for gene ontology enrichment analysis.
Translational efficiency analysis. We used Riborex 91 to perform differential translational efficiency analysis. The underlying engine selected was DESeq2 88 . DESeq2 estimates a single dispersion parameter per gene. However, RNA-Seq and Ribo-Seq libraries can have different dispersion parameters owing to different protocols. We estimated the dispersion parameters for RNA-Seq and Ribo-Seq samples separately and found them to be significantly different (mean difference = 0.04, p-value < . − e 2 2 16). This leads to a skew in translational-efficiency p-value distribution since the estimated null model variance for the Wald test is underestimated. To address this issue, we performed a p-value correction using fdrtool 92 that re-estimates the variance using an empirical bayes approach. Normalized read counts for both Ribo-seq and RNA-seq samples are provided in Supplementary Table S12.
Alternative splicing analysis. Alternative splicing analysis was performed using rMATS 93 . All reads were trimmed using cutadapt 94 with parameters (-u -13 -U -13) to ensure trimmed reads had equal lengths (138 bp). rMATs was run with default parameters in paired end mode (-t paired) and read length set to 138 bp (-len 138) using GENCODE GTF (v25) and STAR index for GRCh38.
Gene ontology (GO) and pathway enrichment analysis. To classify the functions of differentially enriched genes, we performed GO enrichment, and the Reactome pathway 95 analysis using Panther 96 . For both analyses, we considered terms to be significant if BH adjusted p-values weree <0.05, and fold enrichment is >2.0. Further, we used REVIGO 97 to reduce redundancy of the enriched GO terms and visualize the semantic clustering of the identified top-scoring terms. We used STRING database (v10) 98 to construct protein-protein interaction networks and determine associations among genes in a given dataset. The interactions are based on experimental evidence procured from high-throughput experiments text-mining, and co-occurrence. Only high-confidence (0.70) nodes were retained. expression correlation analysis. Gene expression correlation analysis was done using Gliovis 24 using glioblastoma samples (RNAseq) from the TCGA. To select correlated genes, we used Pearson correlation, > .
R 0 3, and p-value < 0.05. A list of genes affecting survival in glioblastoma was downloaded from GEPIA 99 . A list of long non-coding RNAs (lncRNAs) implicated in glioma development was obtained from Lnc2Cancer 32 . Drug-gene interactions were identified using the Drug-Gene Interaction Database 100 .