An atlas of causal and mechanistic drivers of interpatient heterogeneity in glioma

Glioblastoma, causal and mechanistic network model, precision oncology 28 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. 2

1 Grade IV glioma, formerly known as glioblastoma multiforme (GBM) is the most 2 aggressive and lethal type of brain tumor, and its treatment remains challenging in part 3 due to extensive interpatient heterogeneity in disease driving mechanisms and lack of 4 prognostic and predictive biomarkers. Using mechanistic inference of node -edge 5 relationship (MINER) , we have analyzed multiomics profiles from 516 patients and 6 constructed an atlas of causal and mechanistic drivers of interpatient heterogeneity in 7 GBM (gbmMINER) . The atlas has delineated how 30 driver mutations act in a 8 combinatorial scheme to causally influence a network of regulators ( 306 transcription 9 factors and 73 miRNAs ) of 179 transcriptional “programs”, influencing disease 10 progression in patients across 23 disease states. Through extensive testing on 11 independent patient cohorts, we share evidence that a machine learning model trained 12 on activity profiles of programs within gbmMINER significantly augments risk stratification, 13 identifying patients who are super -responders to standard of care and those that would 14 benefit from 2nd line treatments. In addition to providing mechanistic hypotheses regarding 15 disease prognosis, the activity of programs containing targets of 2 nd line treatments 16 accurately predicted efficacy of 28 drugs in killing glioma stem-like cells from 43 patients. 17 Our findings demonstrate that interpatient heterogeneity manifests from differential 18 activities of transcriptional programs, providing actionable strategies for mechanistically 19 characterizing GBM from a systems perspective and developing better prognostic and 20 predictive biomarkers for personalized medicine. 21

22 Glioblastoma multiforme (GBM) is the most commonly diagnosed malignant brain tumor, with an 23 incidence of 5 per 100,000 people in North America1. GBM is among the most difficult tumor types 24 to treat2. Complete tumor resection is only rarely possible due to the highly invasive nature of 25 GBM cells, which spread to the surrounding brain tissue. The current standard of care for GBM 26 includes surgical resection, radiotherapy, and concomitant and maintenance temozolomide 27 (TMZ) chemotherapy3,4. The prognosis of GBM patients who receive standard of care remains 28 dismal, with a median survival time of approximately 15 months, a five -year survival rate of less 29 than 10%, and a recurrence rate of ~90%5. In a randomized Phase III clinical trial, GBM patients 30 treated with radiation, TMZ, and tumor treating fields (TTF) had a median overall survival (mOS) 31 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 3 of 20.5 months, compared to 15.6 months for patients receiving standard treatment without TTF6. 1 Despite the marginal improvements, GBM tumors remain difficult to treat with a poor prognosis2. 2 3 One major challenge for the lack of clinical improvement of GBM tumors is interpatient 4 heterogeneity2. Interpatient heterogeneity can be addressed by subtyping brain tumors to stratify 5 patients and identify the best -matched drugs and therapies for a particular patient or cohort of 6 patients. Molecular profiling analysis of a growing number of patients has revealed the full extent 7 of GBM tumor heterogeneity. The genetic characteristics of GBM , such as the co -deletion of 8 chromosomes 1p and 19q 7–9, mutations in the genes coding tumor protein P53 ( TP53)10–12, 9 isocitrate dehydrogenase 1/2 ( IDH1/2)8,13, epidermal growth factor receptor ( EGFR)14,15, and 10 deletion of the phosphatase and tensin homolog ( PTEN)14,16 have refined the classification of 11 GBM subtypes. Based on newer molecular and genetic data , the 2021 WHO classification of 12 adult gliomas now defines GBM as any IDH wil d-type diffuse glioma with presence of 13 microvascular proliferation or necrosis and presence of TERT promoter mutations, or, gain of 14 chromosome 7 and loss of chromosome 10 (i.e., 7+/10-), or, EGFR amplification. For sake of 15 consistency with prior work especially the TCGA nomenclature , we have decided to use “GBM” 16 with the caveat that some of the tumors in this cohort may not meet the criteria of how GBM is 17 currently defined. To help the reader relate the findings we have provided a mapping of the new 18 classification across all relevant results throughout the paper by using previous annotations17 19 (Supplementary Table 9). Analysis of gene expression and copy number alterations of GBM 20 cataloged in The Cancer Genome Atlas (TCGA) has revealed several molecular subtypes: 21 mesenchymal, proneural, and classical. Distinct clinical outcomes associated with these subtypes 22 offer a stratification scheme to select appropriate therapies, albeit not yet at the resolution of N -23 of-1 personalized medicine. For instance, the proneural subtype is associated with a more 24 favorable outcome, partly due to the high frequency of mutations in the IDH1/2 genes. In contrast, 25 the mesenchymal subtype is associated with worse survival than other subtypes 18–21. Although 26 IDH mutations have been linked to metabolic reprogramming, epigenetic reprogramming , and 27 redox imbalance, the mechanisms by which IDH mutations predict favorable disease outcome s 28 are not clearly understood 22. Similarly, there is no integrated view of how diverse prognostic 29 markers act in combination to influence clinical outcomes, including response to therapies, in an 30 individual GBM patient. 31 32 Dysfunction across multiple pathways significantly shapes etiology and progression of GBM . 33 Consequences of these mutations manifest through a network of interactions marked by non -34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 4 linear dynamics, rendering the prediction of disease progression and therapy responsiveness 1 challenging. To address the complexity of GBM and better characterize interpatient heterogeneity 2 in GBM, our focus was on generating systems-level understanding of how dysfunction in diverse 3 genes acts causally and mechanistically through a network of interactions to ultimately influence 4 oncogenic processes and clinical outcomes. Previously, we developed SYstems Genetic Network 5 Analysis (SYGNAL) to mine multi-omics and clinical outcome datasets for 422 GBM patients in 6 TCGA to construct a causal and mechanistic tra nscriptional regulatory network (CM-TRN) for 7 GBM (gbmSYGNAL)23. gbmSYGNAL delineated putative mechanisms by which 102 somatically 8 mutated genes and pathways causally perturbed regulation by 74 transcription factors (TFs) and 9 39 miRNA regulators of approximately 5,000 genes across 500 disease-associated co-regulated 10 gene modules (regulons). Through analysis of gbmSYGNAL, we uncovered novel insights, such 11 as mechanism by which mutations in NF1 and PIK3CA increase tumor lymphocyte infiltration by 12 causally inducing the expression of IRF1, which upregulates a regulon of 27 genes associated 13 with antigen processing and presentation. Furthermore, the network model generated actionable 14 insights into the formulation of synergistic anti -proliferative and pro-apoptotic interventions with 15 combinations of targeted inhibitors, anti-miRs, miRNA overexpression, and siRNAs against TFs23. 16 Here, we report the construction of a revised TRN, called gbmMINER, which was generated using 17 an updated version of SYGNAL that leverages the Mechanistic Inference of Node-Edge 18 Relationships (MINER) algorithm24 to analyze multi-omics (whole exome sequencing, microarray, 19 and RNA-seq) and clinical (cancer subtype and overall survival) data from a cohort of 516 20 patients. By outperforming cMonkey 225 in accurately grouping 9,728 genes into 3,797 regulons, 21 MINER significantly improved upon gbmSYGNAL by expanding the recall rate of TFs (from 6 to 22 26%) and miRNAs (from 15 to 25%) based on experimental support. Importantly, using a network 23 quantizing technique to cluster regulons based on similar activity profiles, we identified 179 24 transcriptional programs that stratify GBM into 23 disease states. Using the quantized network, 25 we developed prognostic models that were significantly more accurate than a previously best 26 performing model based on gene expression correlates at predicting risk of disease progression. 27 Notably, regulons and programs that significantly contribute d to risk of disease progression 28 provided hypotheses for mechanism s by which known prognostic markers influence clinical 29 outcomes. Finally, through high-throughput ( HTP) drug screening , we demonstrate proof -of-30 concept for prioritizing FDA-approved drugs based on activity profiles of regulons and programs 31 containing relevant drug targets that accurately predicted, without any prior training, the drug 32 sensitivity of 43 patient-derived glioma stem-like cells ( PD-GSCs) to 28 anticancer drugs. 33 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 5 Together the extensive validations demonstrate that gbmMINER represents an atlas of causal 1 and mechanistic drivers of interpatient heterogeneity in GBM. 2 3

4 gbmMINER delineates how well known and novel disease driving mutations causally and 5 mechanistically stratify GBM patients into 23 distinct states 6 We applied a new version of systems genetic network analysis (SYGNAL) that leverage d 7 Mechanistic Inference of Node Edge Relationships (MINER) to generate a TRN by integrating 8 multi-omics data and clinical outcomes from a cohort of 516 GBM patients in the TCGA26. In 9 constructing gbmMINER, a total of 9,728 genes across the cohort were clustered by MINER into 10 3,797 modules of genes (henceforth “regulons”) with significant co-expression across subsets of 11 patients and co -regulated by a common transcription factor or a miRNA (see Methods). The 12 authenticity of regulons was ascertained by reproducing significant co-expression of member 13 genes across at least one of three independent cohorts of 285 27, 252 28 and 80 29 patients 14 (Supplementary Table 1 ). Each regulon was quantized into three network activity states --15 overactive, neutral, or underactive-- based on statistical assessment of whether the expression 16 levels of regulon member genes in a given patient w ere in the upper, middle, or lower thirds of 17 the distribution of expression levels of each respective gene across all patients in the TCGA 18 cohort. A notable advancement over the previous gbmSYGNAL model 23, network quantization 19 (see Methods) in gbmMINER maintained large-scale patterns of gene expression across the 20 cohort, while significantly improving signal to noise24. In contrast to gbmSYGNAL, which identified 21 500 disease -relevant biclusters (avg. regulon size: 36 genes) , gbmMINER identified 1,083 22 disease-relevant regulons (avg. regulon size = 11 genes) based on co-regulation of member 23 genes in at least one independent dataset, and their association with patient survival (Cox HR p-24 value ≤ 0.05) or enrichment for a hallmark of cancer (enrichment p-value ≤ 0.05). (Table 1). To 25 reduce the dimensionality of the TRN, the 3,797 regulons were clustered based on similar activity 26 profiles across all patients into 179 distinct transcriptional programs, of which 58 programs 27 contained the 1,083 disease -relevant regulons ( Figure 1A and Supplementary Table 1 ). 28 Similarly, based on correlated activity profiles of the 179 programs, the 516 patients were stratified 29 into sub-populations of 23 distinct transcriptional states , potentially reflecting different subtypes 30 of GBM disease-relevant expression signatures. 31 32 Using the transcription factor binding site database (TFBS_db)23 and the Framework for Inference 33 of Regulation by miRNAs (FIRM)30, 306 TFs and 83 miRNAs were implicated as likely regulators 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 6 of the 1,083 disease-relevant regulons (Table 1, Methods). Further, somatic mutations within 30 1 genes and 49 pathways were causally associated with both the altered expression of regulators 2 and concordant downstream consequences on mechanistic regulation of 2,049 genes within 3 disease-relevant regulons associated with patient survival. Specifically, the causal influences of 4 the mutations on disease -relevant regulons were confirmed to act through their mechanistic 5 regulators by assessing differential expression of regulons across wild type and mutated patient 6 samples, with and without perturbed expression of the regulators ( Figure 1B , Table 1 and 7 Supplementary Table 2 ). The causal influences of 21 somatic gene mutations in gbmMINER 8 were also previously modeled in gbmSYGNAL (overlap p -value:4.4e-32). Notably, gbmMINER 9 modeled the causal influences of nine well-known driver mutations in GBM , including EGFR, 10 IDH1, NF1, PIK3CA, PIK3R1, PTEN, RB1, TP53, and ATRX (Figure 1C and Supplementary 11 Table 5) as well as 15 mutations with some known association with GBM. gbmMINER implicated 12 causal and mechanistic diseases associations for 6 mutated genes ( DNAH2, APOB, TCHH, 13 HSD17B7P2, KEL, KRTAP4-11) that were not previously associated with GBM. Whereas most 14 prognostic mutations for GBM were previously identified using GWAS, which does not provide 15 mechanistic insights into their role in disease etiology or progression, gbmMINER delineated 16 causal and mechanistic pathways through which mutations in these genes perturb the regulation 17 of TFs and miRNAs, ultimately impacting the expression of downstream genes associated with 18 the disease (Figure 1C and gbmMINER Portal). 19 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 7 1 2 Figure 1. gbmMINER: a systems scale genotype to phenotype map of causal and mechanistic drivers of clinical outcomes in GBM. A). Heatmap of regulon activity across patients reveals distinct transcriptional states wherein patients have very similar regulon activity across the entire network. Disease relevant programs, which are group of regulons with similar activity across patients, are indicated on the left while transcriptional states are shown on the top. Overactive, neutral and underactive regulons are colored in red, white and blue, respectively. B). The inferred gbmMINER TRN is a predictive map that implicates specific somatic mutations in causally modulating the expression of a TF(s) or miRNA(s) that in turn regulates genes within disease-relevant regulons. A summary of the counts for each feature in the gbmMINER TRN is shown. C) Causal and mechanistic relationships for known and novel GBM prognostic markers are well-represented in the model. D) gbmMINER transcriptional programs and E) states stratify risk of disease progression in TCGA cohort patients. Kaplan-Meier plot of survival probabilities for all programs and all states together with comparison of high-risk programs/states (red) versus low-risk programs/states (blue) are shown. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 8 Table 1. Comparison of gbmSYGNAL and gbmMINER Model Features gbmSYGNAL gbmMINER Total regulons 1,830 3,797 Average number of genes per regulon 36 11 Number of disease relevant regulons 500 1,083 Number of disease relevant programs 58 NA Number of somatic gene/pathway mutations 33/69 30/49 Number of disease-relevant TF regulators 74 306 Number of disease-relevant genes 5,193 2,049 Number of miRNA regulators 39 73 TFs validated with CRISPR-Cas9 screening of patient- derived glioma stem-like cells31 29 (p-value=0.24) 137 (p-value=3.1e-5) TFs validated with CRISPR-Cas9 screening of patient- derived glioma stem-like cells32 26 (p-value=0.01) 91 (p-value=0.002) Overlap with GBM-relevant TFs in DisGeNET33 16 (p-value=5.2e-4) 114 (p-value=5.4e-21) miRNAs with experimental -support for association with GBM (HMDD34) 13 (p-value=2.6e-9) 22 (p-value=6.2e-14) Overlap with miRNAs dysregulated in GBM (miR2Disease35) 7 (p-value=4.5e-7) 7 (p-value=4.2e-5) 1 gbmMINER implicates novel transcription factors and miRNAs in GBM. 2 Altogether, 306 TFs were implicated by gbmMINER in the regulation of disease-relevant regulons, 3 which was a significant improvement over gbmSYGNAL, which had previously identified 74 TFs, 4 of which 38 TFs were in both models (p-value = 7.4e-11). Phenotype data for 1,543 TF knockouts 5 (94% of known TFs) from a genome -wide CRISPR -Cas9 screen 31 identified 568 TFs had 6 significantly altered proliferation in at least one of 10 patient-derived glioma stem-like cell lines 7 (PD-GSCs). Of these 568 anti -proliferative TFs, only 29 TFs were in gbmSYGNAL (p = 0.24) , 8 whereas gbmMINER identified a total of 137 TFs (p = 3.1e-5), including 48 TFs that also altered 9 the proliferation of human fetal neural stem cells (Table 1 and Supplementary Table 3). Another 10 CRISPR-Cas9 screen study on 2 PD-GSCs32 also validated 26 and 91 TFs identified in 11 gbmSYGNAL (p=0.01) and gbmMINER (p=0.002) , respectively . In addition, according to the 12 DisGeNET database 33 of disease -to-gene associations, 114 of the 306 TFs identified by 13 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 9 gbmMINER (p = 5.4e-21), as compared to 16 of the 74 TFs identified in gbmSYGNAL (p = 5.2e-1 4), have important functions in GBM (See Supplementary Table 3 for validation results) . In 2 summary, the gbmMINER network implicated 306 TFs in the regulation of 2,049 GBM-relevant 3 genes, recapitulating 114 TFs that had been previously implicated in GBM , and 192 TFs with 4 potentially novel and yet to be characterized roles in modulating disease outcomes in GBM. 5 6 gbmMINER also implicated 73 miRNAs as likely disease drivers based on their enriched binding 7 sites in the 3' UTRs of genes within disease -relevant regulons and their negatively correlated 8 expression levels (R ≤ -0.2, p-value ≤0.05)36, which represents a significant improvement over the 9 39 miRNAs identified by gbmSYGNAL. Two lines of evidence supported the biological and 10 disease relevance of miRNAs in the gbmMINER network. First, 22 of the 73 miRNAs (p-11 value=6.2e-14) were implicated in GBM in HMDD 34, a curated database of evidence -based 12 disease associations of human miRNAs (Table 1). Second, 7 miRNAs were implicated in GBM in 13 the miR2Disease database 35, which documents evidence for miRNAs that are dysregulated in 14 human diseases (p-value=4.2e-5). In total, 25 of the 73 miRNAs had been previously implicated 15 as dysregulated or causally associated with GBM, indicating the potential for discovering novel 16 biology associated with the additional 48 miRNAs identified by gbmMINER. 17 18 A network of master transcriptional regulators governs GBM transcriptional states, whose 19 association with higher disease risk is correlated with increased immune evasion. 20 Based on s imilarity of transcriptional program activities in gbmMINER, the 516 patients in the 21 TCGA cohort were grouped into 23 transcriptional states. Notably, Kaplan-Meier survival analysis 22 demonstrated that both programs and transcriptional states were associated with distinct overall 23 survival outcomes (Figure 1E & Figure 2B). For transcriptional states we also mapped subtype 24 information that was available for a significant number of tumor samples from 343 patients 25 revealing that most of the astrocytomas were included in low to moderate risk states (Figure 2A). 26 Using the ESTIMATE algorithm37, we calculated the ImmuneScore for each transcriptional state 27 and investigated the relationship between immune cell infiltration and transcriptional states. The 28 percentage of tumor-infiltrating immune cells correlated directly with the increased risk association 29 across states (Figure 2C). Additionally, we tested each state for distinct mechanisms of tumor 30 immune evasion by using the TIDE algorithm38. In general, higher risk states correlated with 31 higher immune dysfunction in contrast to immune exclusion indicating that even though there was 32 increased immune cell infiltration in higher risk states, those cells were in a dysfunctional state 33 (Supplementary Figure 9). Based on estimation of the immune cell fractions from the RNASeq 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 10 data39, there was no evidence of overrepresentation of any particular immune cells across states. 1 No correlation was observed also between known prognostic markers, subtypes, and states, 2 indicating that disease-associated pathway and gene mutations alone were not sufficient to 3 determine the transcriptional state of a patient. For example, well known GBM prognostic markers 4 such as PTEN, TP53, and EGFR were enriched across most states, while other markers were 5 mainly enriched in states with a moderate risk association (data not shown) . Th ese findings 6 suggested that the transcriptional states likely manifest from combinations of mutations acting 7 through a complex network of regulatory interactions with system wide consequences on the 8 activity levels of multiple disease-associated transcriptional programs. 9 10 To better understand the transcriptional drivers of the 23 distinct states, we built a TF–TF network 11 and shortlisted master TFs based on the ratio of outdegree/indegree edges (See Methods). We 12 identified a network of 67 master regulators that were implicated in driving or suppressing each 13 state (See Methods; Figure 2D). Most master regulators were specific to a given state, and only 14 a few master regulators , including MAFB, FOX3, RUNX1, STAT4, SOX2, SOX9, E2F1, were 15 shared across two or more states. Notably, knockdowns in 38 of 67 master regulators altered the 16 proliferation of at least one PD-GSC (p= 1.4e-4) in a genome-wide CRISPR-Cas9 screen on 10 17 PD-GSCs31. In addition, 32 of the 67 master regulators (p = 5.4e-9) have important functions in 18 GBM, such as in GSC self-renewal and maintenance (e.g., SOX2 and SOX940,41 and MYC42, and 19 repression of GBM cell differentiation (e.g., E2F143, based on disease-to-gene associations in 20 DisGeNET33. Despite extensive combinatorial control , including the influence of TFs that were 21 uniquely associated as master regulators of some states, there was some concordance between 22 type of influence (activator or repressor) of some master regulators on states with respect to their 23 disease risk association. 24 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 11 1 2 Figure 2. Association of disease risk of each transcriptional state with immune cell filtration and an elaborate network of master transcriptional regulators. A) Distribution of patients with different disease subtypes according to 2021 WHO classification 17 across each state. B) Boxplots of GuanRank risk scores for patients within transcriptional states, rank ordered from low to high median risk (from L to R, respectively). C) Boxplots of ImmuneScores indicating relative immune cell infiltration in tumor s of patients within each state . D) Sixty-seven master TFs act in a combinatorial scheme to drive distinct transcriptional programs that characterize 20 of the 23 transcriptional states of GBM. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 12 Program activities delineate known and novel biological processes underlying disease 1 prognosis 2 Of the 179 programs in gbmMINER, Cox HR analysis implicated 58 as disease -associated, of 3 which 7 programs were associated with low -risk (i.e., programs with negative hazard ratios that 4 predicted good prognosis) and 51 programs were associated with high-risk (i.e., a positive hazard 5 ratio) ( Figure 3A , top panel). While the low -risk programs were enriched for oxidative 6 phosphorylation (Pr -118), G2 -M checkpoint (Pr -61), and Myc targets (Pr -144), the high -risk 7 programs were enriched for hypoxia (15 programs, including Pr -172, Pr-6, Pr-43), epithelial-to-8 mesenchymal transition (20 programs including Pr -172, Pr -32, Pr -6), TNF -α signaling (28 9 programs including Pr-172, Pr-32, Pr-6), and other immune-related processes (Supplementary 10 Figures 1 and 3). These results are consistent with significantly better survival of a newly 11 discovered mitochondrial subtype of GBM that is characterized by oxidative phosphorylation44. In 12 contrast, the mesenchymal subtype was associated with worse survival than other subtypes and 13 was characterized by hypoxia and epithelial-mesenchymal transition18–21. The low-risk programs 14 were significantly associated with GBM markers (IDH1, ATRX, TP53 and PDGFRA) implicated in 15 good prognosis, in contrast to high-risk programs that were enriched for mutations in NF1, a well-16 known negative prognostic marker for GBM ( Figure 3A “Mutations” panel). Notably, relative to 17 low-risk programs a s ignificant number of high -risk programs also had higher ImmuneScores , 18 which was consistent with the established association of increased immune cell infiltration and 19 bad prognosis in GBM45 (Figure 3A, “ImmuneScore” panel). Notably, patients with overactive 20 high-risk programs with high immune score also had a higher immune dysfunction score indicating 21 that their bad prognosis was more likely associated with an impaired immune response . In 22 contrast, we observed that patients with overactive low risk programs had lower immune 23 dysfunction score and they had slightly increased CD4+ (non -regulatory) and CD8+ T-cells with 24 decreased M1 and M2 macrophage levels (Supplementary Figure 10). 25 Through delineation of causal -mechanistic information flow ( CM-flow) maps, gbmMINER 26 uncovered insights into how specific mutations influence clinical outcomes by acting causally 27 through TFs and miRNAs that mechanistically regulate activities of distinct biological processes. 28 For example, 9 regulons within the high-risk program 21 (Pr-21; enriched for genes associated 29 with TNF-a signaling, IL6-JAK-STAT3 signaling and inflammatory response), including regulons 30 R-383, R-1984 and R-3754, were overactive across patients in most high-risk states. 31 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 13 Figure 3. Casual and mechanistic underpinnings of high and low risk association of programs. A) Enrichment of disease hallmarks (upper heatmap) and mutations (lower heatmap) across disease associated programs. For each program, risk association and size (numbers of genes) are shown as barplots (top panels). Additionally, the level of immune cell infiltration for each program is indicated as a boxplot at the bottom (red: median ImmuneScore > 0, gray: median ImmuneScore <= 0). Detailed view of a high -risk programs 21 (B) and 118 (C). The heatmap of the network activity of regulons within each program, with red for overactivity, white for neutral and blue for under activity (top panel). Summary of causal mechanistic flows associated with each program is shown in lower left panel. A dot plot of disease hallmark enrichments is shown in the lower right panel. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 14 gbmMINER predicted that mutations in NF1 causally activated TFs SATB1, ELF1, FOSL1, 1 RUNX3 and ZNF232, which in turn upregulated 35 genes within regulons R-1684, R-332, R-1984, 2 R-196 and R -1761 that were associated with TNF-a signaling, IL6-JAK-STAT3 signaling and 3 inflammatory response , providing a causal and mechanistic hypothesis for the elevated 4 ImmuneScore in these patient samples (Figure 3B). Conversely, the positive prognosis of IDH1 5 mutations could be explained by their predicted causal inhibition of ETV7, which was implicated 6 as an activator of Pr-118 regulons with genes of oxidative phosphorylation (Figure 3C and Figure 7 5A). 8 9 Transcriptional program activities predict survival risk 10 Molecular markers such as IDH1/2 mutations and MGMT promoter methylation sub -stratify 11 patients with significantly better clinical outcomes and response to TMZ ( IDH mutant vs WT: 31 12 months vs 15 months 22, MGMT promoter methylation vs WT: 21.7 months vs 12.7 months 46). 13 However, 8% of IDH1/2 WT patients in the TCGA cohort have had longer survival times (mOS: 14 44.6 months, range: 31 -89 months) than patients with IDH1/2 mutations and 23% of MGMT 15 unmethylated patients survived longer (mOS: 42.7 months, range: 22 -117 months) than MGMT 16 methylated patients. Based on the observation that expression of individual regulons, 17 transcriptional programs, and states were significantly associated with distinct survival outcomes, 18 we explored if these features could serve as better prognostic markers of GBM (Methods). We 19 used ridge regression to predict the risk of median overall survival with program activity (See 20 Methods). First, we transformed patient survival data into GuanRank scores between 0 and 147, 21 with values > 0.5 indicating high-risk and values < 0.5 indicating low-risk. Using the concordance 22 index (C-index) metric48,48,49 we evaluated the p erformance of each model on independent test 23 datasets of 113 patients in TCGA dataset (20% of data not used to train the model) and 150 24 patients from an independent study (“Gravendeel dataset” 28) and 80 patients from a nationwide 25 observational study (“XCELSIOR study”, see Methods). We also compared the gbmMINER risk 26 prediction models to the performance of a 4-gene-panel model that stratifies patients based on a 27 risk score calculated using the sum of weighted expression levels of four autophagy genes : 28 DIRAS3, LGALS8, MAPK8, and STAM49. The performance of the gene panel in risk stratifying 29 patients was significantly lower with C-index values of 0.58, 0.57 and 0.56 for the TCGA, 30 Gravendeel and XCELSIOR datasets, respectively (Figure 4A-B and Table 2, Supplementary 31 Table 8). Notably, the KM plot for TCGA and XCELSIOR patients stratified by the gene panel was 32 not statistically significant. By contrast, performance of the program model was consistently better 33 across all datasets. T he predicted low - and high -risk classes by the program model also 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 15 effectively stratified patients based on their survival outcome s, as shown by significant 1 separations in the survival curves of these individual risk groups in the KM plot ( Figure 4A-B). 2 Furthermore, we assessed performance of the program risk model relative to risk stratification by 3 IDH mutations and MGMT promoter methylation status within the XCELSIOR cohort. As 4 expected, MGMT promoter methylation status stratified IDH Wild type patients into low and high-5 risk groups (Figure 4C) . Remarkably, the gbmMINER program risk model also effectively 6 stratified patients into low and high -risk groups, regardless of their IDH mutation and MGMT 7 promoter methylation status (Figure 4D). Importantly, the program risk model also sub-stratified 8 MGMT methylated patients, identifying low risk patients who were likely super -responders to 9 standard of care, and high risk patients who may benefit from combining standard of care with 2nd 10 line treatments (Figure 4E). 11 12 Better prognosis of IDH1, ATRX and TP53 mutations likely manifests from modulation of 13 ferroptosis through their causal influences on ETV7 and CTCF regulons 14 gbmMINER uncovered many of the well-known positive and negative prognostic markers for GBM 15 that act through a combinatorial scheme to modulate risk-associated programs. For instance, the 16 network model revealed a complex combinatorial scheme in which 15 mutations causally perturb 17 the expression of 12 TFs that mechanistically co-regulate 34 genes (Figure 5). Notably, this CM-18 flow map captured how well-known mutations, including IDH1, TP53, and ATRX, might causally 19 influence clinical outcomes through seven TFs (ETV7, POU3F2, POU3F3, CTCF, HOXA3, 20 MEF2A, and NR3C1) that combinatorially regulate eight genes. Of the seven regulators, POU3F2, 21 a neurodevelopmental TF 41,50, and HOXA3, an activator of aerobic glycolysis, have been 22 implicated in promoting GBM propagation51. The expression of these seven TFs was predicted to 23 be causally influenced by nine somatically mutated genes, of which mutations in IDH1 (p-24 value:2e-10), ATRX (p-value:5e-6), and TP53 (p-value:0.006) were associated with Pr -118 25 overexpression (Figure 5A). 26 27 From the CM flows, we extracted the novel insight that IDH1 and ATRX mutations repress ETV7, 28 thereby activating a super -regulon of eight genes within Pr -118 (Figure 5A -C). Note that by 29 design MINER reconstructs multiple regulatory influences on the same set of co-regulated genes 30 by discovering redundant instances of the same regulon , each associated with a unique 31 regulatory influence. Combining these redundant regulons gives a super -regulon, which is so 32 called because it is implicated to be under the control of a large number of regulators -in this 33 example, 7 TFs were assigned as regulators of the super-regulon of 8 genes. 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 16 Figure 4. gbmMINER-based risk prediction models outperform previously best performing gene panel and performs as good as IDH WT and MGMT Promoter Methylation based stratification. The stratification of predicted low-risk (blue) and high-risk (red) groups by Kaplan-Meier curves for TCGA, Gravendeel, and XCELSIOR datasets across A) A gene panel-based risk prediction model 49 (yellow panel header background and B) The gbmMINER based program risk prediction model (blue panel header background). For each plot C-index, number of patients, mOS and log-rank test p-values for the significance of the survival probabilities between the high and risk groups are shown . (See Supplementary Table 8 for risk model outputs and survival analysis results ). Stratification of IDH Wild Type XCELSIOR cohort patients by C) MGMT Promoter Methylation status and D) gbmMINER program risk model. E) Sub-stratification of IDH Wild type and MGMT Promoter Hypermethylated patients by program risk model. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 17 Across significant numbers of patients in the TCGA dataset, the expression of these 7 TFs was 1 inferred to be causally perturbed by at least 10 somatically mutated genes. While not previously 2 implicated in GBM, ETV7 has been shown to regulate breast cancer stem -like cell features by 3 repressing IFN-response genes52 and activating mTORC3 assembly to promote tumor growth in 4 a mouse model 53. In addition, the model predicted that IDH1 and TP53 mutations act through 5 CTCF to activate the same super-regulon (Figure 5A, C and D and Supplementary Figure 4). 6 This finding explains why loss of CTCF binding sites, which modulates communication between 7 enhancers and promoters, is associated with IDH and TP53 mutated GBM tumors54. Within this 8 super-regulon, two of the eight genes, NCOA4 and CISD1, are related to ferroptosis (enrichment 9 p-value = 0.002)55–58. Specifically, whereas NCOA4 increases free iron levels in cells and 10 promotes ferroptosis, CISD1 reduces mitochondrial iron uptake 57. A third super-regulon gene, 11 MSRB2, encodes methionine-R-sulfoxide reductase a multifunctional enzyme that both 12 scavenges reactive oxygen species and is involved in the biosynthesis of cysteine , a key 13 component of GSH, which regulates ferroptosis57. All genes in the ETV7 regulon were significantly 14 upregulated in IDH1 mutants ( Figure 5E ), a key finding that was also reproduced in an 15 independent cohort59 (Supplementary Figure 5). The mechanistic hypotheses were validated by 16 the CRISPR-Cas9 screen31, which showed that ETV7 significantly altered proliferation in one 17 glioblastoma stem cell, POU3F2 altered proliferation in three glioblastoma stem cells, and CTCF 18 altered proliferation in seven glioblastoma stem cell lines. A number of mechanisms have been 19 proposed for the improved prognosis of patients with IDH1/2 mutated GBM22. Our findings are 20 consistent with these prior explanations and extend our mechanistic understanding of how IDH1/2 21 mutations may causally inhibit ETV7 to activate ferroptosis and improve prognosis (i.e., median 22 survival of ~31 months as compared to ~15 months for IDH wildtype GBM22). 23 24 Further, t he combinatorial influences of IDH1, TP53 and ATRX mutations on regulation of 25 ferroptosis offered a plausible explanation for why IDH1 mutant g liomas are genetically 26 associated with TP53 mutations or ATRX mutations60–62. To investigate these influences on a 27 systems level, we mined gbmMINER for the entire causal-mechanistic network linking the three 28 mutations to ferroptosis -related TFs that were implicated in regulating at least two genes 29 implicated as drivers or suppressors of ferroptosis as per FerrDb database 63. We also required 30 ferroptosis-related TFs to have significantly altered proliferation of at least one PD-GSC in the 31 CRISPR-Cas9 screen 31. Based on these criteria, IDH1, TP53 and ATRX mutations were 32 implicated in causally influencing the expression of 29 ferroptosis -related TFs of which 14 were 33 previously implicated in regulating ferroptosis (Supplementary Table 4). 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 18 Figure 5. IDH1, ATRX and TP53 mutations causally impact ETV7 and CTCF to modulate ferroptosis leading to better prognosis. A) Network diagram of causal and mechanistic influences of mutations and TFs on Program 118. Regulons comprising program 118 are represented as blue boxes with their associated regulators (red triangles) and putative causal genetic abnormalities (green chevrons). Red edges denote activation and blue edges denote inhibition. Regulons with significant overlap of member genes are grouped into super-regulons (yellow, purple and green rectangles). B) ETV7 mediates the causal influence of somatically mutated IDH1 on downstream genes in regulon R-843. ETV7 expression (left boxplot) and expression of genes in R-843 (middle boxplot) increases when IDH1 is mutated. Conditioning R-843 expression on ETV7 abolishes the increase in expression indicating that the causal influence of IDH1 is mediated by ETV7. C) ETV7 mediates the causal influence of somatically mutated ATRX on downstream genes in regulon R -843. D) CTCF mediates the causal influence of s omatically mutated IDH1 on downstream genes in regulon R-3529. E) ETV7 regulon genes are upregulated in IDH1 mutants. F) Combinatorial scheme through which IDH1, ATRX, and TP53 mutations modulate 29 TFs implicated in mechanistically regulating 53 ferroptosis genes. The 23 ferroptosis driver genes are grouped within the yellow shaded box and 30 ferroptosis suppressor genes in the green shaded box. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 19 The 29 TFs were implicated in mechanistically regulating 53 ferroptosis-related genes across 18 1 super-regulons. 23 of the 5 3 genes were implicated as suppressors of ferroptosis and the 2 remainder were implicated as ferroptosis drivers (Figure 5F). For example, a known ferroptosis 3 pathway, HSF1-HSPB1, increases the resistance of cancer cells to ferroptosis through the 4 inhibition of iron uptake 58. In addition to the 14 known regulators of ferroptosis , gbmMINER 5 predicted that this process was also putatively regulated by an additional 15 TFs, including ETV7, 6 HOXA11, LEF1, LHX1, LHX3, NFATC4, BPTF, POU3F3, ELF3, SOX11, CTCF, SOX3, SOX5, 7 SRF, and ELF5 (Supplementary Table 4). Consistent with their gbmMINER predicted roles in 8 activating ferroptosis-driver genes, ELF5, SOX11 and HOXA11 were previously demonstrated to 9 act as tumor suppressor s in many cancers64–68. Conversely, consistent with its gbmMINER -10 predicted suppressor role, NFATC4 is known to promote oncogenesis69. Finally, the integrated 11 analysis demonstrated that mutations in IDH1 (11 ferroptosis-activating and 5 ferroptosis-12 suppressing CM-flows), ATRX (8 ferroptosis-activating and 7 ferroptosis-repressing CM-flows) 13 and TP53 (8 ferroptosis-activating and 5 ferroptosis-suppressing CM-flows) were the major 14 modulators of ferroptosis in GBM. Because IDH mutations can fundamentally change the biology 15 of the disease, and given the low numbers of IDH mutant glioma patients in our cohort, we sought 16 to ascertain generalizability of our findings to a larger cohort of 514 low grade glioma (LGG) 17 patients by analyzing multiomics data from the TCGA-LGG cohort and constructing a lggMINER 18 model. Analysis of the lggMINER model independently confirmed that IDH1 mutations do indeed 19 causally influence regulation of ferroptosis genes in LGG patient tumors by modulating a similar 20 set of TFs (Supplementary Figure 12). In summary, gbmMINER delineated how mutations in 21 IDH1, ATRX and TP53 causally and mechanistically modulate specific cancer related processes, 22 including ferroptosis, to have prognostic implications on clinical outcomes. 23 24 Network quantization and drug -constrained regulon activity uncovers patient -specific 25 disease-network map to enable therapy prioritization. A major hurdle in uncovering predictive 26 biomarkers of drug response for personalized medicine applications is the insufficiency of HTP 27 screening datasets for most drugs across a diversity of patient s, patient-derived xenografts or 28 tumor cells, with paired pre - and post-treatment genomic/transcriptomic profiles. Using patient 29 derived glioma stem-like cells (PD-GSCs), we investigated if activity of transcriptional programs 30 (or regulons) containing drug target(s) would accurately predict drug sensitivity of a patient’s 31 tumor cells. In brief, RNA-Seq profiles of 43 PD-GSCs were analyzed through network 32 quantization with gbmMINER to generate PD-GSC-specific disease network maps. In brief, the 33 entire disease-perturbed network within a PD-GSC was uncovered by treating each regulon as a 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 20 discrete unit and evaluating whether it was overactive, underactive, or neutral (using a p -value 1 cutoff of 0.05) based on the distribution of z -scored expression levels of member genes. Drug 2 constrained regulon activity (DCRA) for a given drug was then calculated as the mean activity of 3 all regulons containing or regulated by its target(s), and compared to the background distribution 4 of DCRA values for that drug across all patient samples in a cohort (e.g., TCGA cohort). In doing 5 so, we used DCRA to estimate the overall status of the disease-relevant networks targeted by a 6 given drug and thereby predict whether a given PD-GSC was likely to be a responder (low IC50) 7 or non-responder (high IC50) to treatment with that drug . Specifically, if DCRA was greater than 8 zero for an antagonist drug (less than or equal to zero for an agonist drug), then the PD-GSC was 9 predicted to be sensitive to that drug. By contrast, if the reverse was true then the PD-GSC was 10 predicted to be a non-responder. Using this approach, we predicted sensitivity of 43 PD-GSCs to 11 62 drugs based on patterns of over - and under -activity of regulons containing cognate drug 12 targets (Methods). Notably, this analysis revealed that the distributions of DCRAs for most drugs 13 across the PD-GSCs was significantly skewed relative to corresponding distributions in the TCGA 14 cohort (Supplementary Fig 6). This skewed distribution could be attributable to inherent 15 increased or decreased drug susceptibility of PD -GSCs or the absence of tumor 16 microenvironment influences. Nonetheless, t he predicted susceptibilities of 43 PD-GSCs to 62 17 drugs was then compared to responder/non-responder classifications based on IC50 values that 18 were experimentally determined in a HTP drug screen (Methods and Figure 6A). 19 20 While a significant number of drugs had variable efficacy across PD -GSCs, many drugs had 21 consistently high (i.e., low IC50) or low (i.e., high IC 50) efficacy across most of the 43 PD-GSCs, 22 which was consistent with their skewed distributions of DCRA values relative to the TCGA cohort. 23 Importantly, the predicted sensitivity of the PD -GSCs were significantly concordant with the 24 experimentally determined IC 50 values. Based on permutation test s, it was determined that at 25 least 30 concordances between predicted and experimentally determined response to each drug 26 would be necessary to be considered statistically significant alignment at an FDR < 0.05. For a 27 dataset of 4 3 PD-GSCs against 6 2 drugs, while we would expect by chance that just 7 drugs 28 would show ≥30 concordances with experimental results, we observed efficacy predictions for 28 29 drugs agreed with experimentally determined sensitivity of up to 43 PD-GSCs (p-value=0, Figure 30 6B, Supplementary Figure 8 ). The distinct DCRA-predicted sensitivity profiles across the se 31 drugs for each PD -GSC was consistent with the HTP screen , which demonstrated significant 32 difference in IC 50 values between drugs that were predicted to be effective or not -effective for 33 each PD-GSC (Figure 6 C). 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 21 1 Given that the MTT assay estimates cell viability, we investigated whether the concordance of 2 PD-GSCs vis-à-vis predicted and measured drug sensitivity was explained by the activity status 3 of drug-target containing regulons enriched for genes associated with the viability-relevant cancer 4 hallmark “evading apoptosis”. Indeed, responder PD-GSCs that were concordant (i.e., with high 5 DCRA and low IC50 values) had a high proportion of drug target containing regulons (59%) that 6 were overactive and enriched for genes associated with “evading apoptosis”. The converse was 7 also true, that is non -responder PD-GSCs that were concordant (i.e., low DCRA and high IC 50 8 values) were associated with a greater proportion (61%) of drug-target containing regulons that 9 were underactive and enriched for “eva ding apoptosis” genes (Supplementary Figure 11). 10 These results suggest that drugs which induce cytotoxic effects have higher concordance 11 between DCRA and MTT viability -assay-based assessment of sensitivity , and that discordant 12 drugs may target other functions that do not necessarily lead to cell death. If so, then the accuracy 13 of drug sensitivity predictions by gbmMINER could turn out to be even higher if IC50 of discordant 14 drugs are estimated using biochemical assays that probe relevant biological functions enriched 15 in regulons containing their respective target(s). 16 17

18 Figure 6. Network quantization uncovers patient-specific disease networks enabling N-of-1 prediction of therapy responsiveness. A) Schematic illustration of the workflow for validating network-based therapy response predictions through concordance analysis against PD -GSC high throughput drug screening results. B) Statistical significance of concurrency was determined through a permutation test. The inset plot shows distribution of the single drug FDR adjusted p-values up to 43 concurrencies. At least 30 concurrencies are needed (red dashed lines) to achieve FDR = 30 PD-GSCs that have concurrent predicted and experimentally-measured sensitivities (green dashed lines), while 28 (red dashed lines) drugs were observed to have >=30 concurrencies in the actual comparison (p-value = 0e-1 as determined with permutation testing ). C) Comparison of the experimentally -determined IC50 values for predicted responders and nonresponders for significantly concordant drugs. p-values were determined with a Mann-Whitney U test. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 22 We have generated a model for characterizing interpatient heterogeneity in GBM. The model was 1 demonstrated to have high mechanistic accuracy based on extensive performance evaluation of 2 independent patient cohorts, disease databases, literature, CRISPR-cas9 phenotypic screens 3 from two independent studies on 12 PD-GSCs and normal neuronal cell lines, and through HTP 4 screening of 62 drugs against 43 PD -GSCs. Findings from these studies demonstrated the 5 remarkable mechanistic accuracy with which gbmMINER has mapped modular architecture of 6 gene-gene co-regulation within regulons, identifying TFs and miRNAs that differentially regulate 7 cancer processes across sub-groups of GBM patients. In so doing, gbmMINER delineated 8 networks through which well-known and novel somatic mutations in genes and pathways 9 influence clinical outcomes and response to therapy . For instance, gbmMINER delineated how 10 IDH1, ATRX and TP53 mutations act through a network of 29 regulators to modulate 5 3 11 ferroptosis genes across three super -regulons, affecting treatment response and prognosis. In 12 this manner, we can now interpret how 30 driver mutations modulate expression of 306 TF and 13 83 miRNA regulators to alter activity patterns of 179 transcriptional programs, 58 of which were 14 determined to be disease-relevant. In so doing, gbmMINER has stratified GBM patients into 23 15 states with varied disease risk, which is mechanistically governed by a distinct TRN. 16 17 A p rognostic model developed using TRN status (i.e., activity patterns of programs) was 18 significantly more accurate and robust at predicting clinical outcomes (mOS) , particularly in 19 comparison to a 4-gene panel70. The higher performance of the transcriptional program -based 20 risk model is attributable to the network quantization approach, which likely overcomes issues of 21 variability in gene expression that may manifest from multiple sources, including technical issues 22 such as sample processing and variabilities across profiling technologies, and biological issues 23 such as disease stage of sampling, and composition of cells in patient sample. Specifically, since 24 the program model averages expression values of co-regulated genes within modules, it likely 25 overcomes noise and missing values in expression patterns for individual genes. Furthermore , 26 the program model also provided insights into the biology of disease -driving mechanisms that 27 distinguish low-risk patients from high-risk patients, regardless of their IDH mutation and MGMT 28 methylation status . For instance, we discovered that IDH1, TP53 and ATRX mutations act 29 combinatorially to modulate up to 29 TFs and activate ferroptosis, which explain s the better 30 prognosis of IDH mutants (31 months of mOS compared to 15 months mOS for IDH WT)22. 31 32 Three principal mechanisms have been put forth to explain why IDH mutations predict favorable 33 disease outcomes22. First, IDH mutations drive metabolic reprogramming by increasing oxidative 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 23 metabolism and suppressing glutamine metabolism 22. Indeed, we observed genes of oxidative 1 phosphorylation were significantly enriched in program 118 modulated by IDH1 mutations (p-2 value= 8.1e-07). Second, IDH mutations epigenetically reprogram tumor cells by inhibiting DNA 3 and histone demethylation, leading to hyp ermethylation globally or at specific targets 22. 4 Consistent with this hypothesis, HM27 methylation profiling revealed that 3 ferroptosis-related 5 TFs (LHX3, ETV7, PPRX2) were significantly hypermethylated at their promoters in IDH1 mutant 6 GBM71 (Supplementary Figure 7). At least in the case of ETV7, the hypermethylation of this TF 7 was consistent with the repression of its expression levels in IDH1 mutated GBM patients. Third, 8 IDH mutations activate ferroptosis through accumulation of excessive reactive oxygen species22. 9 Our findings integrate and extend the three mechanisms into a unified view that explains how IDH 10 mutations perturb a network of TFs to modulate oxidative phosphorylation and DNA methylation, 11 resulting in activation of ferroptosis. 12 13 We inferred CM-flows in the ferroptosis network from just 14 IDH-mutated tumor samples, within 14 the TCGA-GBM cohort, that should be reclassified as low-grade gliomas. We reproduced these 15 findings across risk-stratified patient subgroups in the 514 LGG patient dataset, including the 16 causal influence of IDH1 mutation on ETV7 and its downstream regulons. In addition to 17 showcasing the sensitivity of gbmMINER in delineating disease -driving mechanisms in a small 18 subgroup of patients within a larger cohort, this finding also underscores generalizability of the 19 mechanisms by which IDH1 acts through a combinatorial network to activate ferroptosis. While 20 there is evidence in literature that both IDH1 and TP53 regulate ferroptosis , the primary 21 mechanism was unknown 72. In fact , a recent study discovered that although t he acetylation -22 defective mutant TP53 3KR lost its ability to induce cell senescence, apoptosis, and cell -cycle 23 arrest, it was still able to inhibit tumorigenesis by inducing ferroptosis 73. Activating the ferroptosis 24 pathway is being explored as a promising therapeutic strategy for many cancers, including GBM, 25 as cancer cells tend to have higher metabolic demand for iron74. The r egulatory network of 26 ferroptosis delineated in gbmMINER could facilitate such efforts. We demonstrated that quantized 27 DCRA values of regulons and programs in gbmMINER accurately predict drug sensitivity based 28 on the distribution of DCRA values across a reference cohort, without the need for machine 29 learning on a large training dataset. This approach is also generalizable towards repurposing 30 FDA-approved drugs with targets within disease-relevant regulons and programs in gbmMINER, 31 for matching treatment regimen to a patient based on their DCRA profiles, and for identifying and 32 prioritizing new drug targets. 33 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 24

1 S.T. is supported by R01-CA259469, Washington Research Foundation Grant No: GA-00117 and 2 Andy Hill Cancer Research Endowment (CARE) Fund . W.W., K.K. and C.M. are supported by 3 R01-CA259469. M.M. is supported by R01 -CA259469 and Washington Research Foundation 4 Grant No: GA -00117. Y.H. is supported by R01 -CA259469 and Andy Hill Cancer Research 5 Endowment (CARE) Fund . A.N. is supported by NSF -2150265. A.P.P. is supported by R01 -6 NS119650, the Burroughs Wellcome Career Award for Medical Scientists, and Discovery Grant 7 from the Kuni Foundation. P.H. and H.L. were funded by the Ben and Catherine Ivy Foundation 8 and are currently supported by R01 -CA259469 and philanthropic funding from the Swedish 9 Medical Center Foundation. C.C. is supported by R01-CA259469. J.P. was funded by a fellowship 10 from the NIH (F32 -CA247445) and currently supported by NIH grant R01 -CA259469. N.S.B. is 11 supported by R01-CA259469, Washington Research Foundation Grant No: GA-00117 and Andy 12 Hill Cancer Research Endowment (CARE) Fund. 13 We thank Timothy J. Martins and the UW Quellos High-throughput Screening Core for advice and 14 services in running the HTP screens. 15 16 DECLARATION OF INTERESTS 17 NB is a founder and owns equity (founder shares) and ST owns equity in Sygnomics, which is a 18 privately owned biotechnology company focusing on commercializing the Systems Genetic 19 Network Analysis (SYGNAL) technology for personalized medicine applications. AP is a 20 consultant for and has equity interest in Sygnomics. 21 22 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 25 1 SUPPLEMENTARY INFORMATION 2 Supplementary Tables 3 Supplementary Table 1: Summary of disease relevant regulons and validations across 4 independent datasets 5 Supplementary Table 2: List and details of disease relevant causal mechanistic flows 6 Supplementary Table 3: Validation for the GBM association of 306 TFs. 7 Supplementary Table 4: TFs were revealed to be associated with ferroptosis in published 8 studies 9 Supplementary Table 5: Known and new prognostic markers identified in gbmMINER 10 Supplementary Table 6: FDR values for a drug showing n (count) number of concurrent 11 instances. 12 Supplementary Table 7: p-values for observing n (count) number of drugs showing at least 30 13 concurrent instances. 14 Supplementary Table 8: Risk prediction model and Survival analysis output 15 Supplementary Table 9: Clinical information about patient outcomes and tumor subtype s for 16 TCGA cohort patients 17 18 19 Supplementary Figures 20 Supplementary Figure 1: Enriched molecular hallmarks for significant risk-associated programs. 21 Supplementary Figure 2: Enrichment of genetic programs for Hallmarks of Cancer. 22 Supplementary Figure 3: KM stratification curves for Program 118. 23 Supplementary Figure 4: An insight that TP53 mutations activate CTCF and activate ferroptosis 24 is extracted from program 118. 25 Supplementary Figure 5: Data from Trautwein e t al., supports that ETV7 regulon genes are 26 upregulated in IDH1 mutants. 27 Supplementary Figure 6: Comparison of PDGSC and TCGA training cohort drug constrained 28 network activity (DCRA). 29 Supplementary Figure 7. LHX3, ETV7 and PPRX2 have significantly higher methylation at their 30 promoters in IDH1 mutants. 31 Supplementary Figure 8. Distribution of the concurrencies for all the drugs tested in PD-GSC 32 high throughput screen. 33 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 26 Supplementary Figure 9. Distribution of Immune Dysfunction scores and immune cells across 1 states. 2 Supplementary Figure 10. Distribution of Immune Dysfunction scores and immune cells across 3 programs 4 Supplementary Figure 11. Activity status of “evading apoptosis” -enriched regulons containing 5 targets of drugs that were concordant between DCRA-predicted and MTT-assay determined drug 6 sensitivity across relevant responder and non-responder PD-GSCs. “Concurrency percent” refers 7 to the proportion of concurrency events between predicted and observed drug sensitivity. 8 Supplementary Figure 12. Overlap of gbmMINER and lggMINER network for ferroptosis related 9 IDH mutation causal mechanistic flows. 10 11 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 27 Supplementary Figures 1 2 3 Supplementary Figure 1. Enriched molecular hallmarks for significant risk-associated programs. Low-risk programs 4 are mainly enriched for Oxidative Phosphorylation and Myc Targets. High -risk programs are mainly enriched for 5 hypoxia, EMT, TNF-alpha Signaling, and immune responses. 6 7 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 28 1 Supplementary Figure 2. Enrichment of genetic programs for Hallmarks of Cancer. Enrichment for each program 2 was calculated by using semantic similarity of the GO terms associated with genes in the programs. 3 4 High Risk Low Risk Hallmarks of Cancer Limitless_Replicative_Potential Tissue_Invasion_And_Metastasis Sustained_Angiogenesis Reprogramming_Energy_Metabolism Tumor_Promoting_Inflammation Self_Sufficiency_In_Growth_Signals Insensitivity_To_Antigrowth_Signals Genome_Instability_And_Mutation Evading_ImmuneDetection Evading_Apoptosis 172 32 16 35 6 11 43 3 13 36 20 22 18 17 12 29 30 28 0 76 85 50 10 21 15 38 177 34 168 155 39 165 1 131 137 170 174 25 27 142 104 5 24 9 140 37 124 4 7 169 55 58 61 118 128 144 161 68 Enrichment 0.8 0.85 0.9 0.95 1 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 29 1 Supplementary Figure 3. KM stratification curves for Program 118. Program 118 is a top program that stratifies GBM 2 patients in low-risk programs. 3 4 5 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 30 1 2 Supplementary Figure 4. An insight that TP53 mutations activate CTCF and activate ferroptosis is extracted from 3 program 118. A) Extracted insight for TP53 mutations from program 118. B) CTCF expression is significantly 4 upregulated when TP53 is somatically mutated. C) The expression of CTCF regulon is significantly upregulated when 5 TP53 is somatically mutated. D) Conditioning CTCF regulon expression on the expression of CTCF destroys the link 6 between CTCF regulon and mutations in TP53. 7 8 9 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 31 1 Supplementary Figure 5. Data from Trautwein et al., supports that ETV7 regulon genes are upregulated in IDH1 2 mutants. 3 4 5 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 32 Supplementary Figure 6: Comparison of PDGSC and TCGA training cohort drug constrained network activity (DC RA). A) Heatmaps of DCRA for both the TCGA and PDGSC samples for the drugs contained in the GBM drug screen. The drugs (rows) are organized for each drug based on the highest number of responder concurrencies from the top and number of highest nonresponder concurrencies from the bottom. The column between the two heatmaps identifies those drugs where gbmMINER was able to make statistically significant drug sensitivity predictions (>= 30 concurrencies based on an FDR of 0.05). The drugs were binned into 3 gr oups based on those drugs where gbmMINER correctly predicted at least either 30 responder classifications (top 6 rows), at least 30 nonresponder classifications (bottom 19 rows), and those in neither of the first two categories (middle 37 drugs). B) Density plots of the TCGA and PDGSC DC RA for the three binned groups of drugs (responders, nonresponders, and neutral). The PDGSC density for the responder bin is clearly skewed upwards in comparison to the TCGA density which explains why these PDGSCs were sens itive to these particular drugs. Similarly, the PDGSC density for the nonresponder samples was skewed lower than the TCGA density and thus shows why there was such poor sensitivity to those particular drugs. Alternatively, the neutral bin did not show much difference between the PDGSC and TCGA densities. Consequently, we did not observe a preponderance of concurrency skewed towards either responders or nonresponders. 1 TCGA PDGSC DAUNORUBICIN IDARUBICIN EPIRUBICIN HYDROCHLORIDE MITOXANTRONE HYDROCHLORIDE BORTEZOMIB METHOTREXATE DOXORUBICIN CLOFARABINE CLADRIBINE DACTINOMYCIN TOPOTECAN HYDROCHLORIDE DIGOXIN TENIPOSIDE GEMCITABINE AURANOFIN PINDOLOL PITAVASTATIN VINBLASTINE SULFATE VINCRISTINE SULFATE VINORELBINE PACLITAXEL COLCHICINE ETOPOSIDE BLEOMYCIN FLUOROURACIL DOCETAXEL MYCOPHENOLATE MOFETIL METHYLENE BLUE FLUDARABINE PHOSPHATE CYTARABINE ATORVASTATIN CALCIUM FLUVASTATIN LOVASTATIN ROSUVASTATIN CALCIUM SIMVASTATIN SORAFENIB BROMPHENIRAMINE MALEATE IRINOTECAN ETHACRYNIC ACID VORINOSTAT ERLOTINIB METFORMIN CLONIDINE HYDROCHLORIDE CLOMIPRAMINE LOMUSTINE ITRACONAZOLE CELECOXIB SODIUM TETRADECYL SULFATE SERTRALINE HALOPERIDOL TRIFLUOPERAZINE CAPTOPRIL DISULFIRAM HYDROXYCHLOROQUINE MEFLOQUINE MEMANTINE NELFINAVIR SIROLIMUS SUNITINIB FLUPHENAZINE HYDROCHLORIDE NISOLDIPINE VALPROIC ACID DRCA −1 −0.5 0 0.5 1 Significant yes no samples nonresponder (19 drugs) neutral (37 drugs) responder (6 drugs) −1.0 −0.5 0.0 0.5 1.0 0.0 0.5 1.0 0 1 2 3 4 0 1 2 3 4 DCRA density data pdgsc tcga A B All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 33 1 Supplementary Figure 7. LHX3, ETV7 and PPRX2 have significantly higher methylation at their promoters in IDH1 2 mutants. 3 4 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 34 Supplementary Figure 8. Distribution of the concurrencies for all the drugs tested . Green line marks number of concurrencies for the FDR < 0.05 and green dots indicate drugs with significant concurrencies. 1 2 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 35 Supplementary Figure 9. Distribution of Immune Dysfunction scores and immune cells across states. A) Distribution of the patient subtypes across each state identified by using 2021 WHO annotations provided by (Zakharova 2022). B) Boxplot of the transcriptional states rank ordered low to high risk from left to right, respectively) based on the distribution of the median Guan risk scores of patients included in the state. C) High risk states were associated with higher ImmuneScores indicating increased percentages of immune cell infiltration in the tumor. D) Higher risk states were with higher Immunescore are also associated with higher Immune Dysfunction score as determined by TIDE algorithm. E) Distribution of various immune cell types across states as determined by quantiSeq algorithm. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 36 1 Supplementary Figure 10. Distribution of Immune Dysfunction scores and immune cells across programs. A) Program risk scores, B) The level of immune cell infiltration for each program is indicated as a boxplot (red: median ImmuneScore > 0, gray: median ImmuneScore 0, gray: median Immune Dysfunction <= 0). D) Distribution of main immune cell types across patients with overactive program activity. 2 3 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 37 Supplementary Figure 11. Activity status of “evading apoptosis”-enriched regulons containing targets of drugs that were concordant between DCRA-predicted and MTT-assay determined drug sensitivity across relevant responder and non-responder PD-GSCs. “Concurrency percent” refers to the proportion of concurrency events between predicted and observed drug sensitivity. 1 2 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 38 Supplementary Figure 12. Overlap of gbmMINER and lggMINER network for ferroptosis related IDH mutation causal mechanistic flows. 1 2 3 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 39

1 Glioblastoma somatic mutations, miRNA and mRNA expression data: Training data for 2 GBM-MINER network construction was acquired from The Cancer Genome Atlas (TCGA). Gene 3 expression microarray data were normalized and Z-scored. Gene expression RNA-seq data were 4 TMM-normalized and Z -scored. A total of 536 patient tumors and four post -mortem controls 5 containing gene expression data for 10,263 genes combined with both RNASeq and microarray 6 Z-scored data were used to train the MINER model ( Supplementary Table 9). The miRNA 7 expression data for 536 miRNA samples were quantified by microarray and used for causal and 8 mechanistic inference of miRNA regulation of co-expressed gene clusters identified by MINER as 9 described before 23. A gene was considered to be mutated in a patient’s tumor if a somatic 10 mutation was observed that modified the gene coding sequence (missense, nonsense, frameshift, 11 in-frame insertion or deletion, splice site, modified translation start site, introduced new start site, 12 or removed stop codon). An NCI -nature pathway was considered to be mutated in a patient’s 13 tumor if at least one gene member of the pathway had a nonsynonymous somatic mutation (all 14 TCGA mutations were not annotated as silent). Mutated genes and pathways used in causal 15 network inference were required to be mutated in at least 5% of the patient tumors. 16 17 Glioblastoma clinical data for survival analyses: Glioblastoma clinical phenotypes for survival 18 analyses were acquired from The Cancer Genome Atlas (TCGA) using Broad Firehose. Clinical 19 information about patient outcomes and tumor subtype was used in Regulon post -processing 20 (Supplementary Table 9). 21 22 Independent glioblastoma validation gene expression cohorts: Raw Affymetrix HGU133 23 Plus 2.0 CEL files were either downloaded from the NCBI Gene Expression Omnibus (GEO; 24 GSE7696 and GSE16011)28,29 or EMBL-EBI ArrayExpress (E-MTAB-3073)27. The CEL files for 25 each of the three studies used in our analysis were background subtracted and quantile 26 normalized using the justRMA function from the affy package of Bioconductor 27 (www.bioconductor.org) using an alternative CDF file75. The 18,023 non redundant probes were 28 selected by collapsing highly correlated probes (correlation coefficient ≥ 0.8) that mapped to the 29 same Entrez ID as they are likely measuring the same transcript isoform. Only grade II, III and IV 30 gliomas were included in the analyses. In addition to these three microarray datasets, RNA -seq 31 data generated as part of the Ivy Glioblastoma Atlas Project ( 32 https://glioblastoma.alleninstitute.org/static/home) was TMM -normalized and used as an RNA -33 seq expression-based validation cohort. All validation data were Z-scored after normalization. 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 40 1 MINER disease relevant regulon membership prediction validation: Regulons were 2 considered significantly coexpressed if the variance explained by first principal component within 3 coherent members was greater than or equal to 0.3 and was significantly larger than random 4 samples (empirical p-value ≤ 0.05 and replication in at least one independent dataset empirical 5 p-value ≤ 0.05; Each of the 3,764 regulons were post -processed to discover: (i) coexpression 6 quality via variance explained by first principal component within coherent members (empirical p-7 value < 0.05 and variance explained ≥ 0.3), (ii) validation of co-expression quality in at least one 8 independent GBM cohort (empirical p -value ≤ 0.05); (iii) putative TF or miRNA (via the FIRM 9 pipeline) regulators associated with regulon with a minimum spearman’s correlation of 0.3 10 between regulators and eigen gene within coherent members (only negative correlation between 11 miRNA and regulon) in the training dataset; (iv) survival analysis with regulon eigen genes (cox 12 hazards ratio (HR) p-value ≤ 0.05); (v) independent validation of survival association in one of the 13 three microarray datasets (cox HR p -value ≤ 0.05); (vi) association with hallmarks of cancer 14 (Jiang-Conrath Semantic Similarity Score ≥ 0.8) 30,76; and (vii) putative somatic gene mutations 15 and pathways with T-test p-value <= 0.05 for regulators and regulon eigen gene between mutated 16 and wild -type samples. The regulons were filtered using the above -mentioned criteria for 17 validation of co -expression, mechanistic inference, causal inference and disease relevance 18 (Supplementary Table 1). Recall rate of TFs between gbmSYGNAL and gbmMINER was 19 calculated as the number of TFs supported by CRISPR and DisGeNet divided by total number of 20 TFs. The recall rate for miRNA was calculated as the number of miRNAs supported by HMDD 21 divided by the total number of miRNAs. 22 23 Network quantization: Calculation of discrete regulon activity. Calculation of discrete regulon 24 activity was performed as described before24. Briefly, for each standardized patient sample, gene 25 expressions were organized from the least to the most expressed and divided into three equal 26 segments: lower, middle, and upper thirds. Assuming effective normalization of gene expression, 27 we can propose a null hypothesis suggesting that a random assortment of genes lacking any co-28 expression relationship will distribute evenly across these thirds. To assess this, we employ a 29 binomial distribution with a probability parameter (p) of 1/3 to model the likelihood of 'k' genes 30 falling within the same third, given a selection of 'N' genes, where 'N' is greater than or equal to 31 'k'. A standard p-value of 0.05 is applied as a threshold for rejecting the null hypothesis, indicating 32 that the selected gene set lacks co-expression. Genes surpassing this threshold in the lower third 33 are designated as "under -active," while those exceeding it in the upper third are labeled "over -34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 41 active." Instances falling outside these categories are denoted as "neutral." Consequently, we 1 generate a matrix assigning values of { -1, 0, 1} to represent the discrete activity of all regulons 2 across all samples. This methodology for calculating discrete regulon activity has previously been 3 employed in the MINER analysis for patients diagnosed with multiple myeloma 24, with 4 adaptations made accordingly. 5 6 Regulon clustering to identify super-regulons. For the 12 regulons in Program 118 and the 55 7 ferroptosis-related regulons, genes in each regulon were first converted into 1 (appearing) or 0 8 (non-appearing) using bag of words encoding method 77. Then k -means clustering was used to 9 cluster the encoded regulons. The optimal number of clusters (super-regulons) is determined by 10 the combination of Elbow method and Silhouette method78. 11 12 Risk prediction. We used ridge regression models to predict risk scores. TCGA dataset was 13 randomly split into a training set (70%), validation set (10%), and test set (20%). The models were 14 learned from the training set and evaluated on the test set and two independent datasets. The 15 validation set was used to select the optimal hyperparameter in the ridge regression models. 16 Based on three types of input data, gene expression (TPM), regulon activity (-1,0,1), and program 17 activity (-1,0,1), we learned three different models and referenced them as gene, regulon, and 18 program models, respectively. An extra parameter related to the number of features 19 (genes/regulons/programs) was learned from the validation set, similar to the hyperparameter in 20 ridge regression. The models output GuanScore, where a score 0.5 is predicted to be high risk. Models were evaluated using the concordance 22 index as well as Kaplan-Meier curves to estimate patient survival for the predicted risk groups in 23 all datasets. The 4-gene-panel model used for comparison calculates the risk score using the sum 24 of weighted expression levels of four autophagy genes [(0.1052 × [DIRAS3]) + (0.2152 × 25 [LGALS8]) + (-0.3603 × [MAPK8]) + (-0.2851 × [STAM])]. 26 27 Identification of master regulators from TF-TF network 28 For each state, we first identify a list of TFs causally associated with the state and then build TF–29 TF network in three steps. First, Least Absolute Shrinkage and Selection Operator (LASSO) 30 regression models are generated for each TF in the list to predict its expression using a subset of 31 the other TFs in the list as predictors. Specifically, the TF list is a subset to include only TFs with 32 a binding site for the target TF in TFBSDB or a CHiP- seq database. The second step is to prune 33 the LASSO models to minimize the number of TF predictors necessary to maintain the same level 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 42 of predictive accuracy. Finally, the LASSO coefficients of each predictor TF for each target TF are 1 defined as weighted edges to connect the TFs in a network. After building TF-TF network, we rank 2 all TFs by the ratio of outdegree/indegree and pick the top TFs as master regulators, because 3 master regulator is defined as the gene at the top of the regulatory hierarchy, which should not 4 be affected by the regulation of any other genes. A script for the whole process is provided. 5 6 Drug mapping pipeline: OpenTargets (with 6515 drugs in the database) and the Cancer Cell 7 Line Encyclopedia (CCLE) (with 24 drugs as part of drug screening) were the two main sources 8 used to map therapies to our updated GBM network. Drugs were mapped to genes belonging to 9 all levels of causal-mechanistic flows, namely, genes containing somatic mutations, transcription 10 factors associated with regulons, and individual regulon genes, if the gene is targeted by the drug. 11 12 Patient samples, cultivation of PD-GSCs, and high throughput drug screening. WHO grade 13 IV glioblastoma tumors were obtained from surgeries performed at Swedish Medical Center 14 (Seattle, WA) according to institutional guidelines . The freshly resected tumor tissues were 15 processed to generate PD-GSCs as described79. For the 43 PD-GSCs used in this study, 39 were 16 established as neurosphere cultures and maintained in ultra -low attachment dishes (Corning). 17 The other 4 PD-GSCs grew as adherent monolayers in T75 flasks pre-treated with laminin (1:100; 18 Sigma). All 43 PD-GSCs were cultured in serum-free media consisting of Neurobasal Medium-A 19 (GibcoTM) with 2.0% (v/v) B-27 serum-free supplement minus vitamin A (GibcoTM), 20 ng/mL EGF 20 (PeproTech Inc.), 20 ng/mL FGF -2 (PeproTech Inc.), 20 ng/mL insulin (Sigma), 1 mM sodium 21 pyruvate (Corning), 2 mM L -glutamine (GibcoTM) and 1% Antibiotic -Antimycotic (GibcoTM). PD-22 GSC cultures were maintained at 37°C, 5% CO2, 1% O2, with culture pH monitored with phenol 23 red. Cultures were refed every 2-3 days. High throughput (HTP) screening assays (IC50 studies) 24 were conducted essentially as described79. In brief, PD-GSCs were added to either laminin coated 25 384-well plates (adherent monolayer cultures) or polyhema coated 384-well plates 26 (neurospheres) at a density of ~2000 cells per well using a Thermo Scientific Matrix WellMate . 27 Following incubation of PD-GSCs overnight at 37 oC, drugs were added to individual wells using 28 the CyBi-Well vario liquid handler (Analytik Jena) and plates were incubated at 37oC for 96 hours. 29 Drugs were added at concentrations ranging from 0.0005 – 10 µM to generate 10 -point dose 30 response curves for each PD-GSC culture (n=43) and each candidate drug (n=62). Following the 31 incubation period, CellTiter-Glo (Promega) was added to individual wells per the manufacturer’s 32 recommendations and luminescence was measured on a Perkin Elmer EnVision plate reader. 33 HTP measurements were collected in duplicate or triplicate (depending on PD-GSC cell counts), 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 43 and data was corrected for background luminescence. PD-GSC cultures tested were within 6 1 passages from the initial GSC enrichment from the original tumor biopsy. 2 3 Determination of IC50 values for each drug-PD-GSC combination. The IC50 values used in our 4 analysis were generated by applying the ‘drm’ (“drc” R package80) log-logistic 4 parameter model 5 (LL4) to the dose -response curve data (dose values and associated percent viability across 3 6 replicates) where all 4 parameters were estimated. We then estimated the IC 50 values manually 7 based on the LL4 model produced by the ‘drm’ function. The percent viabilities across the curve’s 8 dose range were calculated by inputting values that covered the drug screen dose range into the 9 LL4 regression model which then returned the associated percent viabilities at each point in the 10 estimated dose response curve. The IC 50 value was then identified by finding the dosage that 11 returned a percent viability closest to the 50% viability value. The 50% viability value used to 12 identify the IC50 dose value was calculated as half the mean percent viability of the 3 replicates at 13 the lowest dose range. 14 15 Calculating the concurrency between the gbmSYGNAL predicted and experimentally 16 determined drug sensitivities of PD-GSCs. We first classified each PD-GSC/drug combo as to 17 whether they responded to a drug (responder or nonresponder ) for the GBM drug screen, and 18 the gbmMINER drug sensitivity prediction. The IC 50 was used to classify the PD-GSC’s in the 19 drug screen while the drug constrained network activity was used to predict drug sensitivity 20 (DCRA) for gbmMINER. Drugs are mapped to gbmMINER by identifying all regulons that either 21 contain, or are regulated by, at least one of a drug’s targets and those regulons are then 22 considered to be mapped to that drug. The gbmMINER DCRA for each PD-GSC/drug 23 combination is calculated as the mean of all drug-mapped regulons and will range from -1 to 1. If 24 both the GBM drug screen and the gbmMINER predictions show the same response (responder 25 or nonresponder) then they are concurrent between the two measurements. 26 27 The statistical significance of the number of concurrency instances for each drug (x/43) and the 28 entire drug screen was generated with permutation testing described in the following steps. 29 30 Single drug FDR adjusted p-values for 1 to 43 number of concurrencies. 31 1. Permute the classification labels (responder, nonresponder) for both the drug screen and 32 the gbmMINER predictions. This creates a null distribution of responses based on the original 33 results. 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 44 2. Calculate the concurrencies of the permuted null distribution in the same manner as with 1 the original data for each drug/PD-GSC combination. 2 3. Repeat the permutation 10000 times. 3 4. For each drug (62 drugs) calculate the probability of seeing at least n (n=1 to 43) instances 4 of concurrency across the 43 PD-GSCs exposed to that drug. 5 a. Sum number of times each drug showed at least n concurrent instances and divide by 6 10000. 7 b. This produces a probability for the number of concurrencies ranging from 1 to 43 for each 8 drug and can be interpreted as a p-value. This created a 62 by 43 matrix where the rows are the 9 62 drugs, and the columns are the number of potential number of concurrencies out of the 43 PD-10 GSCs. Each cell has the p-value associated with the number of concurrencies of greater than or 11 equal than the nth number of concurrencies. 12 c. FDR was applied across each column (62 drugs) to account for multiple testing. 13 5. The FDRs were highly consistent across all drugs for all n number of concurrent instances 14 so the mean of the FDRs across the 62 drugs were used as the final FDR for each number of 15 concurrencies (Supplementary Table 6). 16 6. A drug showing at least 30 (from Supplementary Table 6 ) instances of concurrency 17 passed an FDR cutoff of 0.05. 18 19 Significance of entire drug screen 20 1. We found 28 drugs that showed at least 30 instances of concurrency in our original 21 analysis. 22 2. We calculated the number of times we found a drug showing at least 30 instances of 23 concurrency in each iteration of our permutation process and divided by 10,000. 24 a. This gives us a probability, interpretable as a p-value, that shows how likely we are to see 25 at least n (1 to 62) number of drugs that showed at least 30 instances of concurrency 26 (Supplementary Table 7). 27 3. A literal p-value of 0 was generated for the likelihood of seeing at least 28 drugs showing 28 30 instances of concurrency. In fact, the highest number of drugs that showed at least 30 29 instances of agreement out of 10000 iterations was 7 out of 62 drugs. 30 31 N-of-1 therapy mapping for a patient. Causal mechanistic (CM) flows associated with disease-32 relevant regulons are shortlisted from all CM flows inferred by the gbmMINER. Therapy 33 candidates are also shortlisted as those used as anti-cancer drugs in the past and/or used in GBM 34 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 45 trials in the past according to clinicaltrials.gov. Regulon activity is calculated for a particular patient 1 to show which regulons are over - or under-active among the disease -relevant regulons. Their 2 associated CM flows are also shortlisted. Short-listed drugs were mapped to patient-specific CM 3 flows. If a disease-relevant regulon is associated with a positive RMST difference (RMST higher 4 in under-active patients compared to over -active patients), and if the regulon is over -active, as 5 well as if the patient shows a mutation that upregulates regulon activity, shutting down the regulon 6 genes might inhibit the cancer -causing pathway associated with the regulon genes. Hence, if a 7 regulon gene is a target of a shortlisted drug that is an inhibitor of the gene, and if the gene is 8 highly expressed, the drug is a good therapeutic candidate for the patient and may be used for 9 further testing. Similar conclusions of potential agonists or antagonists can be drawn with different 10 combinations for RMST differences, regulon activity, and mutations affecting regulon activity. 11 Similar inferences have been made for the TF regulator target of a drug. In the case of mutations, 12 if mutations present in a patient cause gene amplification and the mutation is causally associated 13 with upregulation of a regulon that shows poor survival in upregulated patients, an inhibitor acting 14 against the mutated gene might be a promising target. 15 16 XCELSIOR Observational Study . The XCELSIOR real -world registry is a pan-cancer 17 observational study (NCT03793088) that started in 2019 . Patients may enroll via electronic 18 consent and sign a blanket HIPAA release to permit aggregation of electronic medical records 19 and structured EMR data from all sites of care. Unstructured text from clinic narratives were 20 utilized as source documents for annotation in an electronic database. Raw genomics and 21 transcriptomics data (FASTQ files) were collected from commercial labs. Dates of diagnosis were 22 abstracted from pathology reports. Patient death dates were verified for accurate overall survival 23 calculations. For patients without known death date, censoring was based on the most recent of 24 the document date or last medication record available in EMR. In the current study, clinical and 25 omics data from 80 patients were used. 26

27 1. Oronsky, B., Reid, T.R., Oronsky, A., Sandhu, N., and Knox, S.J. (2020). A Review of Newly 28 Diagnosed Glioblastoma. Frontiers Oncol 10, 574012. 10.3389/fonc.2020.574012. 29 2. Park, J.H., Lomana, A.L.G. de, Marzese, D.M., Juarez, T., Feroze, A., Hothi, P., Cobbs, C., 30 Patel, A.P., Kesari, S., Huang, S., et al. (2021). A Systems Approach to Brain Tumor Treatment. 31 Cancers 13, 3152. 10.3390/cancers13133152. 32 3. Lau, D., Magill, S.T., and Aghi, M.K. (2014). Molecularly targeted therapies for recurrent 33 glioblastoma: current and future targets. Neurosurg Focus 37, E15. 34 10.3171/2014.9.focus14519. 35 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 46 4. Thakkar, J.P., Dolecek, T.A., Horbinski, C., Ostrom, Q.T., Lightner, D.D., Barnholtz-Sloan, 1 J.S., and Villano, J.L. (2014). Epidemiologic and molecular prognostic review of glioblastoma. 2 Cancer Epidem Biomar 23, 1985–1996. 10.1158/1055-9965.epi-14-0275. 3 5. Weller, M., Cloughesy, T., Perry, J.R., and Wick, W. (2013). Standards of care for treatment 4 of recurrent glioblastoma--are we there yet? Neuro-oncology 15, 4–27. 10.1093/neuonc/nos273. 5 6. Stupp, R., Taillibert, S., Kanner, A., Read, W., Steinberg, D.M., Lhermitte, B., Toms, S., 6 Idbaih, A., Ahluwalia, M.S., Fink, K., et al. (2017). Effect of Tumor-Treating Fields Plus 7 Maintenance Temozolomide vs Maintenance Temozolomide Alone on Survival in Patients With 8 Glioblastoma: A Randomized Clinical Trial. JAMA 318, 2306–2316. 10.1001/jama.2017.18718. 9 7. Bello, M.J., Leone, P.E., Nebreda, P., Campos, J.M. de, Kusak, M.E., Vaquero, J., Sarasa, 10 J.L., Garcia-Miguel, P., Queizan, A., Hernandez-Moneo, J.L., et al. (1995). Allelic status of 11 chromosome 1 in neoplasms of the nervous system. Cancer Genet Cytogen 83, 160–164. 12 10.1016/0165-4608(95)00064-v. 13 8. Sabha, N., Knobbe, C.B., Maganti, M., Omar, S.A., Bernstein, M., Cairns, R., Cako, B., 14 Deimling, A. von, Capper, D., Mak, T.W., et al. (2014). Analysis of IDH mutation, 1p/19q 15 deletion, and PTEN loss delineates prognosis in clinical low-grade diffuse gliomas. Neuro-16 oncology 16, 914–923. 10.1093/neuonc/not299. 17 9. Smith, J.S., Alderete, B., Minn, Y., Borell, T.J., Perry, A., Mohapatra, G., Hosek, S.M., 18 Kimmel, D., O’Fallon, J., Yates, A., et al. (1999). Localization of common deletion regions on 1p 19 and 19q in human gliomas and their association with histological subtype. Oncogene 18, 4144–20 4152. 10.1038/sj.onc.1202759. 21 10. Mellai, M., Monzeglio, O., Piazzi, A., Caldera, V., Annovazzi, L., Cassoni, P., Valente, G., 22 Cordera, S., Mocellini, C., and Schiffer, D. (2012). MGMT promoter hypermethylation and its 23 associations with genetic alterations in a series of 350 brain tumors. J Neuro-oncol 107, 617–24 631. 10.1007/s11060-011-0787-y. 25 11. Wakimoto, H., Tanaka, S., Curry, W.T., Loebel, F., Zhao, D., Tateishi, K., Chen, J., Klofas, 26 L.K., Lelic, N., Kim, J.C., et al. (2014). Targetable Signaling Pathway Mutations Are Associated 27 with Malignant Phenotype in IDH-Mutant Gliomas. Clin Cancer Res 20, 2898–2909. 28 10.1158/1078-0432.ccr-13-3052. 29 12. Chang, J., Zhong, R., Tian, J., Li, J., Zhai, K., Ke, J., Lou, J., Chen, W., Zhu, B., Shen, N., et 30 al. (2018). Exome-wide analyses identify low-frequency variant in CYP26B1 and additional 31 coding variants associated with esophageal squamous cell carcinoma. Nat Genet 50, 338–343. 32 10.1038/s41588-018-0045-8. 33 13. Lewandowska, M.A., Furtak, J., Szylberg, T., Roszkowski, K., Windorbska, W., Rytlewska, 34 J., and Jóźwicki, W. (2014). An Analysis of the Prognostic Value of IDH1 (Isocitrate 35 Dehydrogenase 1) Mutation in Polish Glioma Patients. Mol Diagn Ther 18, 45–53. 36 10.1007/s40291-013-0050-7. 37 14. Smith, J.S., Tachibana, I., Passe, S.M., Huntley, B.K., Borell, T.J., Iturria, N., O’Fallon, J.R., 38 Schaefer, P.L., Scheithauer, B.W., James, C.D., et al. (2001). PTEN Mutation, EGFR 39 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 47 Amplification, and Outcome in Patients With Anaplastic Astrocytoma and Glioblastoma 1 Multiforme. Jnci J National Cancer Inst 93, 1246–1256. 10.1093/jnci/93.16.1246. 2 15. Gil-Benso, R., Lopez-Gines, C., Benito, R., López-Guerrero, J.A., Callaghan, R.C., Pellín, 3 A., Roldán, P., and Cerdá-Nicolás, M. (2007). Concurrent EGFR amplification and TP53 4 mutation in glioblastomas. Clin Neuropathol 26, 224–231. 10.5414/npp26224. 5 16. Sasaki, H., Zlatescu, M.C., Betensky, R.A., Ino, Y., Cairncross, J.G., and Louis, D.N. (2001). 6 PTEN Is a Target of Chromosome 10q Loss in Anaplastic Oligodendrogliomas and PTEN 7 Alterations Are Associated with Poor Prognosis. Am J Pathology 159, 359–367. 10.1016/s0002-8 9440(10)61702-6. 9 17. Zakharova, G., Efimov, V., Raevskiy, M., Rumiantsev, P., Gudkov, A., Belogurova-10 Ovchinnikova, O., Sorokin, M., and Buzdin, A. (2022). Reclassification of TCGA Diffuse Glioma 11 Profiles Linked to Transcriptomic, Epigenetic, Genomic and Clinical Data, According to the 2021 12 WHO CNS Tumor Classification. Int. J. Mol. Sci. 24, 157. 10.3390/ijms24010157. 13 18. Phillips, H.S., Kharbanda, S., Chen, R., Forrest, W.F., Soriano, R.H., Wu, T.D., Misra, A., 14 Nigro, J.M., Colman, H., Soroceanu, L., et al. (2006). Molecular subclasses of high-grade 15 glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in 16 neurogenesis. Cancer Cell 9, 157–173. 10.1016/j.ccr.2006.02.019. 17 19. Verhaak, R.G.W., Hoadley, K.A., Purdom, E., Wang, V., Qi, Y., Wilkerson, M.D., Miller, 18 C.R., Ding, L., Golub, T., Mesirov, J.P., et al. (2010). Integrated Genomic Analysis Identifies 19 Clinically Relevant Subtypes of Glioblastoma Characterized by Abnormalities in PDGFRA, 20 IDH1, EGFR, and NF1. Cancer Cell 17, 98–110. 10.1016/j.ccr.2009.12.020. 21 20. Huse, J.T., Phillips, H.S., and Brennan, C.W. (2011). Molecular subclassification of diffuse 22 gliomas: Seeing order in the chaos. Glia 59, 1190–1199. 10.1002/glia.21165. 23 21. Wang, Q., Hu, B., Hu, X., Kim, H., Squatrito, M., Scarpace, L., deCarvalho, A.C., Lyu, S., Li, 24 P., Li, Y., et al. (2017). Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes 25 Associates with Immunological Changes in the Microenvironment. Cancer Cell 32, 42-56.e6. 26 10.1016/j.ccell.2017.06.003. 27 22. Han, S., Liu, Y., Cai, S.J., Qian, M., Ding, J., Larion, M., Gilbert, M.R., and Yang, C. (2020). 28 IDH mutation in glioma: molecular mechanisms and potential therapeutic targets. Brit J Cancer 29 122, 1580–1589. 10.1038/s41416-020-0814-x. 30 23. Plaisier, C.L., O’Brien, S., Bernard, B., Reynolds, S., Simon, Z., Toledo, C.M., Ding, Y., 31 Reiss, D.J., Paddison, P.J., and Baliga, N.S. (2016). Causal Mechanistic Regulatory Network for 32 Glioblastoma Deciphered Using Systems Genetics Network Analysis. Cell Syst 3, 172–186. 33 10.1016/j.cels.2016.06.006. 34 24. Wall, M.A., Turkarslan, S., Wu, W.J., Danziger, S.A., Reiss, D.J., Mason, M.J., Dervan, A.P., 35 Trotter, M.W.B., Bassett, D., Hershberg, R.M., et al. (2021). Genetic program activity delineates 36 risk, relapse, and therapy responsiveness in multiple myeloma. Npj Precis Oncol 5, 60. 37 10.1038/s41698-021-00185-0. 38 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 48 25. Reiss, D.J., Plaisier, C.L., Wu, W.-J., and Baliga, N.S. (2015). cMonkey2: Automated, 1 systematic, integrated detection of co-regulated gene modules for any organism. Nucleic Acids 2 Res 43, e87–e87. 10.1093/nar/gkv300. 3 26. Brennan, C.W., Verhaak, R.G.W., McKenna, A., Campos, B., Noushmehr, H., Salama, S.R., 4 Zheng, S., Chakravarty, D., Sanborn, J.Z., Berman, S.H., et al. (2013). The Somatic Genomic 5 Landscape of Glioblastoma. Cell 155, 462–477. 10.1016/j.cell.2013.09.034. 6 27. Madhavan, S., Zenklusen, J.C., Kotliarov, Y., Sahni, H., Fine, H.A., and Buetow, K. (2009). 7 Rembrandt: helping personalized medicine become a reality through integrative translational 8 research. Mol Cancer Res 7, 157–167. 10.1158/1541-7786.mcr-08-0435. 9 28. Gravendeel, L.A., Kouwenhoven, M.C., Gevaert, O., Rooi, J.J. de, Stubbs, A.P., Duijm, J.E., 10 Daemen, A., Bleeker, F.E., Bralten, L.B., Kloosterhof, N.K., et al. (2009). Intrinsic gene 11 expression profiles of gliomas are a better predictor of survival than histology. Cancer Res 69, 12 9065–9072. 10.1158/0008-5472.can-09-2307. 13 29. Murat, A., Migliavacca, E., Gorlia, T., Lambiv, W.L., Shay, T., Hamou, M.F., Tribolet, N. de, 14 Regli, L., Wick, W., Kouwenhoven, M.C., et al. (2008). Stem cell-related “self-renewal” signature 15 and high epidermal growth factor receptor expression associated with resistance to concomitant 16 chemoradiotherapy in glioblastoma. J Clin Oncol 26, 3015–3024. 10.1200/jco.2007.15.7164. 17 30. Plaisier, C.L., Pan, M., and Baliga, N.S. (2012). A miRNA-regulatory network explains how 18 dysregulated miRNAs perturb oncogenic processes across diverse cancers. Genome Res 22, 19 2302–2314. 10.1101/gr.133991.111. 20 31. MacLeod, G., Bozek, D.A., Rajakulendran, N., Monteiro, V., Ahmadi, M., Steinhart, Z., 21 Kushida, M.M., Yu, H., Coutinho, F.J., Cavalli, F.M.G., et al. (2019). Genome-Wide CRISPR-22 Cas9 Screens Expose Genetic Vulnerabilities and Mechanisms of Temozolomide Sensitivity in 23 Glioblastoma Stem Cells. Cell Reports 27, 971-986 e9. 10.1016/j.celrep.2019.03.047. 24 32. Toledo, C.M., Ding, Y., Hoellerbauer, P., Davis, R.J., Basom, R., Girard, E.J., Lee, E., 25 Corrin, P., Hart, T., Bolouri, H., et al. (2015). Genome-wide CRISPR-Cas9 Screens Reveal Loss 26 of Redundancy between PKMYT1 and WEE1 in Glioblastoma Stem-like Cells. Cell Reports 13, 27 2425–2439. 10.1016/j.celrep.2015.11.021. 28 33. Pinero, J., Sauch, J., Sanz, F., and Furlong, L.I. (2021). The DisGeNET cytoscape app: 29 Exploring and visualizing disease genomics data. Comput Struct Biotechnology J 19, 2960–30 2967. 10.1016/j.csbj.2021.05.015. 31 34. Huang, Z., Shi, J., Gao, Y., Cui, C., Zhang, S., Li, J., Zhou, Y., and Cui, Q. (2019). HMDD 32 v3.0: a database for experimentally supported human microRNA-disease associations. Nucleic 33 Acids Res 47, D1013–D1017. 10.1093/nar/gky1010. 34 35. Jiang, Q., Wang, Y., Hao, Y., Juan, L., Teng, M., Zhang, X., Li, M., Wang, G., and Liu, Y. 35 (2009). miR2Disease: a manually curated database for microRNA deregulation in human 36 disease. Nucleic Acids Res 37, D98–D104. 10.1093/nar/gkn714. 37 36. Baek, D., Villen, J., Shin, C., Camargo, F.D., Gygi, S.P., and Bartel, D.P. (2008). The impact 38 of microRNAs on protein output. Nature 455, 64–71. 10.1038/nature07242. 39 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 49 37. Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., 1 Treviño, V., Shen, H., Laird, P.W., Levine, D.A., et al. (2013). Inferring tumour purity and stromal 2 and immune cell admixture from expression data. Nat Commun 4, 2612. 10.1038/ncomms3612. 3 38. Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., and Liu, X.S. (2020). Large-scale 4 public data reuse to model immunotherapy response and resistance. Genome Med. 12, 21. 5 10.1186/s13073-020-0721-z. 6 39. Finotello, F., Mayer, C., Plattner, C., Laschober, G., Rieder, D., Hackl, H., Krogsdam, A., 7 Loncova, Z., Posch, W., Wilflingseder, D., et al. (2019). Molecular and pharmacological 8 modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. 9 Genome Med. 11, 34. 10.1186/s13073-019-0638-6. 10 40. Wang, Z., Xu, X., Liu, N., Cheng, Y., Jin, W., Zhang, P., Wang, X., Yang, H., Liu, H., and Tu, 11 Y. (2018). SOX9-PDK1 axis is essential for glioma stem cell self-renewal and temozolomide 12 resistance. Oncotarget 9, 192–204. 10.18632/oncotarget.22773. 13 41. Suva, M.L., Rheinbay, E., Gillespie, S.M., Patel, A.P., Wakimoto, H., Rabkin, S.D., Riggi, N., 14 Chi, A.S., Cahill, D.P., Nahed, B.V., et al. (2014). Reconstructing and reprogramming the tumor-15 propagating potential of glioblastoma stem-like cells. Cell 157, 580–594. 16 10.1016/j.cell.2014.02.030. 17 42. Tateishi, K., Iafrate, A.J., Ho, Q., Curry, W.T., Batchelor, T.T., Flaherty, K.T., Onozato, M.L., 18 Lelic, N., Sundaram, S., Cahill, D.P., et al. (2016). Myc-Driven Glycolysis Is a Therapeutic 19 Target in Glioblastoma. Clin Cancer Res 22, 4452–4465. 10.1158/1078-0432.ccr-15-2274. 20 43. Godoy, P., Donaires, F.S., Montaldi, A.P.L., and Sakamoto-Hojo, E.T. (2021). Anti-21 Proliferative Effects of E2F1 Suppression in Glioblastoma Cells. Cytogenet Genome Res 161, 22 372–381. 10.1159/000516997. 23 44. Garofano, L., Migliozzi, S., Oh, Y.T., D’Angelo, F., Najac, R.D., Ko, A., Frangaj, B., Caruso, 24 F.P., Yu, K., Yuan, J., et al. (2021). Pathway-based classification of glioblastoma uncovers a 25 mitochondrial subtype with therapeutic vulnerabilities. Nat Cancer 2, 141–156. 10.1038/s43018-26 020-00159-4. 27 45. Huang, S., Song, Z., Zhang, T., He, X., Huang, K., Zhang, Q., Shen, J., and Pan, J. (2020). 28 Identification of Immune Cell Infiltration and Immune-Related Genes in the Tumor 29 Microenvironment of Glioblastomas. Front. Immunol. 11, 585034. 10.3389/fimmu.2020.585034. 30 46. Thon, N., Kreth, S., and Kreth, F.-W. (2013). Personalized treatment strategies in 31 glioblastoma: MGMT promoter methylation status. OncoTargets Ther. 6, 1363–1372. 32 10.2147/ott.s50208. 33 47. Huang, Z., Zhang, H., Boss, J., Goutman, S.A., Mukherjee, B., Dinov, I.D., Guan, Y., and 34 Consortium, P.R.O.-A.A.C.T. (2017). Complete hazard ranking to analyze right-censored data: 35 An ALS survival study. Plos Comput Biol 13, e1005887. 10.1371/journal.pcbi.1005887. 36 48. Harrell, F.E., Lee, K.L., and Mark, D.B. (1996). MULTIVARIABLE PROGNOSTIC MODELS: 37 ISSUES IN DEVELOPING MODELS, EVALUATING ASSUMPTIONS AND ADEQUACY, AND 38 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 50 MEASURING AND REDUCING ERRORS. Stat. Med. 15, 361–387. 10.1002/(sici)1097-1 0258(19960229)15:43.0.co;2-4. 2 49. Wang, Y., Zhao, W., Xiao, Z., Guan, G., Liu, X., and Zhuang, M. (2020). A risk signature 3 with four autophagy-related genes for predicting survival of glioblastoma multiforme. J Cell Mol 4 Med 24, 3807–3821. 10.1111/jcmm.14938. 5 50. Fiscon, G., Conte, F., Licursi, V., Nasi, S., and Paci, P. (2018). Computational identification 6 of specific genes for glioblastoma stem-like cells identity. Sci Rep-uk 8, 7769. 10.1038/s41598-7 018-26081-5. 8 51. Yang, R., Zhang, G., Dong, Z., Wang, S., Li, Y., Lian, F., Liu, X., Li, H., Wei, X., and Cui, H. 9 (2022). Homeobox A3 and KDM6A cooperate in transcriptional control of aerobic glycolysis and 10 glioblastoma progression. Neuro-oncology 25, 635–647. 10.1093/neuonc/noac231. 11 52. Pezzè, L., Meškytė, E.M., Forcato, M., Pontalti, S., Badowska, K.A., Rizzotto, D., 12 Skvortsova, I.-I., Bicciato, S., and Ciribilli, Y. (2021). ETV7 regulates breast cancer stem-like cell 13 features by repressing IFN-response genes. Cell Death Dis 12, 742. 10.1038/s41419-021-14 04005-y. 15 53. Harwood, F.C., Geltink, R.I.K., O’Hara, B.P., Cardone, M., Janke, L., Finkelstein, D., Entin, 16 I., Paul, L., Houghton, P.J., and Grosveld, G.C. (2018). ETV7 is an essential component of a 17 rapamycin-insensitive mTOR complex in cancer. Sci Adv 4, eaar3938. 10.1126/sciadv.aar3938. 18 54. Sesé, B., Ensenyat-Mendez, M., Iñiguez, S., Llinàs-Arias, P., and Marzese, D.M. (2021). 19 Chromatin insulation dynamics in glioblastoma: challenges and future perspectives of precision 20 oncology. Clin Epigenetics 13, 150. 10.1186/s13148-021-01139-w. 21 55. Hong, Y., Lin, M., Ou, D., Huang, Z., and Shen, P. (2021). A novel ferroptosis-related 12-22 gene signature predicts clinical prognosis and reveals immune relevancy in clear cell renal cell 23 carcinoma. Bmc Cancer 21, 831. 10.1186/s12885-021-08559-0. 24 56. Liu, H., Hu, H., Li, G., Zhang, Y., Wu, F., Liu, X., Wang, K., Zhang, C., and Jiang, T. (2020). 25 Ferroptosis-Related Gene Signature Predicts Glioma Cell Death and Glioma Patient 26 Progression. Frontiers Cell Dev Biology 8, 538. 10.3389/fcell.2020.00538. 27 57. Shi, Z., Zhang, L., Zheng, J., Sun, H., and Shao, C. (2021). Ferroptosis: Biochemistry and 28 Biology in Cancers. Frontiers Oncol 11, 579286. 10.3389/fonc.2021.579286. 29 58. Tang, D., Chen, X., Kang, R., and Kroemer, G. (2021). Ferroptosis: molecular mechanisms 30 and health implications. Cell Res 31, 107–125. 10.1038/s41422-020-00441-1. 31 59. Trautwein, C., Zizmare, L., Mäurer, I., Bender, B., Bayer, B., Ernemann, U., Tatagiba, M., 32 Grau, S.J., Pichler, B.J., Skardelly, M., et al. (2021). Tissue metabolites in diffuse glioma and 33 their modulations by IDH1 mutation, histology and treatment. Jci Insight 7, e153526. 34 10.1172/jci.insight.153526. 35 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 51 60. Sledzinska, P., Bebyn, M.G., Furtak, J., Kowalewski, J., and Lewandowska, M.A. (2021). 1 Prognostic and Predictive Biomarkers in Gliomas. Int J Mol Sci 22, 10373. 2 10.3390/ijms221910373. 3 61. Murnyak, B., and Huang, L.E. (2021). Association of TP53 Alteration with Tissue Specificity 4 and Patient Outcome of IDH1-Mutant Glioma. Cells 10, 2116. 10.3390/cells10082116. 5 62. Liu, X.Y., Gerges, N., Korshunov, A., Sabha, N., Khuong-Quang, D.A., Fontebasso, A.M., 6 Fleming, A., Hadjadj, D., Schwartzentruber, J., Majewski, J., et al. (2012). Frequent ATRX 7 mutations and loss of expression in adult diffuse astrocytic tumors carrying IDH1/IDH2 and 8 TP53 mutations. Acta Neuropathol 124, 615–625. 10.1007/s00401-012-1031-3. 9 63. Zhou, N., Yuan, X., Du, Q., Zhang, Z., Shi, X., Bao, J., Ning, Y., and Peng, L. (2023). FerrDb 10 V2: update of the manually curated database of ferroptosis regulators and ferroptosis-disease 11 associations. Nucleic Acids Res 51, D571–D582. 10.1093/nar/gkac935. 12 64. Xu, X., Chang, X., Li, Z., Wang, J., Deng, P., Zhu, X., Liu, J., Zhang, C., Chen, S., and Dai, 13 D. (2015). Aberrant SOX11 promoter methylation is associated with poor prognosis in gastric 14 cancer. Cell. Oncol. 38, 183–194. 10.1007/s13402-015-0219-7. 15 65. Qu, X., Li, Q., Tu, S., Yang, X., and Wen, W. (2021). ELF5 inhibits the proliferation and 16 invasion of breast cancer cells by regulating CD24. Mol. Biol. Rep. 48, 5023–5032. 17 10.1007/s11033-021-06495-7. 18 66. Wang, L., Cui, Y., Sheng, J., Yang, Y., Kuang, G., Fan, Y., Jin, J., and Zhang, Q. (2017). 19 Epigenetic inactivation of HOXA11, a novel functional tumor suppressor for renal cell 20 carcinoma, is associated with RCC TNM classification. Oncotarget 8, 21861–21870. 21 10.18632/oncotarget.15668. 22 67. QU, Y., ZHOU, C., ZHANG, J., CAI, Q., LI, J., DU, T., ZHU, Z., CUI, X., and LIU, B. (2014). 23 The metastasis suppressor SOX11 is an independent prognostic factor for improved survival in 24 gastric cancer. Int. J. Oncol. 44, 1512–1520. 10.3892/ijo.2014.2328. 25 68. Se, Y.-B., Kim, S.H., Kim, J.Y., Kim, J.E., Dho, Y.-S., Kim, J.W., Kim, Y.H., Woo, H.G., Kim, 26 S.-H., Kang, S.-H., et al. (2016). Underexpression of HOXA11 Is Associated with Treatment 27 Resistance and Poor Prognosis in Glioblastoma. Cancer Res. Treat. 49, 387–398. 28 10.4143/crt.2016.106. 29 69. Zhong, Q.-H., Zha, S.-W., Lau, A.T.Y., and Xu, Y.-M. (2022). Recent knowledge of NFATc4 30 in oncogenesis and cancer prognosis. Cancer Cell Int. 22, 212. 10.1186/s12935-022-02619-6. 31 70. Prasad, B., Tian, Y., and Li, X. (2020). Large-Scale Analysis Reveals Gene Signature for 32 Survival Prediction in Primary Glioblastoma. Mol Neurobiol 57, 5235–5246. 10.1007/s12035-33 020-02088-w. 34 71. Cerami, E., Gao, J., Dogrusoz, U., Gross, B.E., Sumer, S.O., Aksoy, B.A., Jacobsen, A., 35 Byrne, C.J., Heuer, M.L., Larsson, E., et al. (2012). The cBio Cancer Genomics Portal: An Open 36 Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov. 2, 401–404. 37 10.1158/2159-8290.cd-12-0095. 38 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint 52 72. Zhuo, S., He, G., Chen, T., Li, X., Liang, Y., Wu, W., Weng, L., Feng, J., Gao, Z., and Yang, 1 K. (2022). Emerging role of ferroptosis in glioblastoma: Therapeutic opportunities and 2 challenges. Front. Mol. Biosci. 9, 974156. 10.3389/fmolb.2022.974156. 3 73. Mou, Y., Wang, J., Wu, J., He, D., Zhang, C., Duan, C., and Li, B. (2019). Ferroptosis, a 4 new form of cell death: opportunities and challenges in cancer. J Hematol Oncol 12, 34. 5 10.1186/s13045-019-0720-y. 6 74. Mitre, A.-O., Florian, A.I., Buruiana, A., Boer, A., Moldovan, I., Soritau, O., Florian, S.I., and 7 Susman, S. (2022). Ferroptosis Involvement in Glioblastoma Treatment. Medicina 58, 319. 8 10.3390/medicina58020319. 9 75. Zhang, J., Finney, R.P., Clifford, R.J., Derr, L.K., and Buetow, K.H. (2005). Detecting false 10 expression signals in high-density oligonucleotide arrays by an in silico approach. Genomics 85, 11 297–308. 10.1016/j.ygeno.2004.11.004. 12 76. Hanahan, D., and Weinberg, R.A. (2011). Hallmarks of Cancer: The Next Generation. Cell 13 144, 646–674. 10.1016/j.cell.2011.02.013. 14 77. Juluru, K., Shih, H.-H., Murthy, K.N.K., and Elnajjar, P. (2021). Bag-of-Words Technique in 15 Natural Language Processing: A Primer for Radiologists. RadioGraphics 41, 210025. 16 10.1148/rg.2021210025. 17 78. Saputra, D.M., Saputra, D., and Oswari, L.D. (2020). Effect of Distance Metrics in 18 Determining K-Value in K-Means Clustering Using Elbow and Silhouette Method. Proc. Sriwij. 19 Int. Conf. Inf. Technol. Appl. (SICONIAN 2019), 341–346. 10.2991/aisr.k.200424.051. 20 79. Hothi, P., Martins, T.J., Chen, L., Deleyrolle, L., Yoon, J.-G., Reynolds, B., and Foltz, G. 21 (2012). High-Throughput Chemical Screens Identify Disulfiram as an Inhibitor of Human 22 Glioblastoma Stem Cells. Oncotarget 3, 1124–1136. 10.18632/oncotarget.707. 23 80. Ritz, C., Baty, F., Streibig, J.C., and Gerhard, D. (2015). Dose-Response Analysis Using R. 24 PLoS ONE 10, e0146021. 10.1371/journal.pone.0146021. 25 26 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprintthis version posted April 7, 2024. ; https://doi.org/10.1101/2024.04.05.24305380doi: medRxiv preprint