HSPA5 and FGFR1 genes in the mesenchymal subtype of glioblastoma can improve a treatment efficacy

ABSTRACT

Tyrosine kinase inhibitors (TKIs) have emerged as a potential treatment strategy for glioblastoma multiforme (GBM). However, their efficacy is limited by various drug resistance mechanisms. To devise more effective treatments for GBM, genetic characteristics must be considered in addition to pre-existing treatments. We performed an integrative analysis with heterogeneous GBM datasets of genomic, transcriptomic, and proteomic data from DepMap, TCGA and CPTAC. We found that poor prognosis was induced by co-upregulation of heat shock protein family A member 5 (HSPA5) and fibroblast growth factor receptor 1 (FGFR1). Co-up regulation of these two genes could regulate the PI3K/AKT pathway. GBM cell lines with co-upregulation of these two genes showed higher drug sensitivity to PI3K inhibitors. In the mesenchymal subtype, the co-upregulation of FGFR1 and HSPA5 resulted in the most malignant subtype of GBM. Furthermore, we found this newly discovered subtype was correlated with homologous recombination deficiency (HRD) In conclusion, we discovered novel druggable candidates within the group exhibiting co-upregulation of these two genes in GBM, suggest potential strategies for combination therapy.

KEYWORDS:

• Glioblastoma multiforme
• heat shock protein family A member 5
• fibroblast growth factor receptor 1
• mesenchymal subtype
• combination treatment

Introduction

Glioblastoma multiforme (GBM) is a highly malignant brain tumor with a poor prognosis, characterized by rapid growth, invasiveness, and resistance to treatment. Recent advances in high-throughput sequencing technology have significantly progressed the understanding of cancer genomes, leading to a radical advancement in the molecular classification of cancer. Prominent mutations, such as MGMT promoter methylation, p53, RB1, and IDH1 have been recognized in studies of GBM (Wakimoto et al. 2014; Chen et al. 2017; Zhang et al. 2019). Furthermore, GBM has been intricately classified based on transcriptome profiling according to cell type, enabling personalized treatment strategies for patients based on subtypes.

Based on transcriptional profiling characteristics, GBM has been classified into subtypes: classical, neural, proneural, and mesenchymal (Sidaway 2017). The classical subtype is characterized by the expression signature genes EGFR, AKT2, SMO, GAS1, GLI2, NOTCH3, JAG1, and LFNG (Verhaak et al. 2010). EGFR amplification is associated with the proliferation and survival of tumor cells through the PI3K/AKT/mTOR and RAS/RAF/MEK/ERK pathways, with mutated known genes including PTEN, CHKN2, and PDGFRA (Verhaak et al. 2010; Park et al. 2019; Saito et al. 2019). For the proneural subtype, the expression signature includes PDGFRA, OLIG2, DDL3, SOX2, and NKX2-2, reflecting characteristic expression features. Mutation information indicates alterations in key neurodevelopmental pathways due to mutations in TP53, PI3K, IDH1, and PDGFRA (Verhaak et al. 2010; Alentorn et al. 2012; Steponaitis and Tamasauskas 2021). For the neural subtype, the expression signature includes MBP/MAL, NEFL, SLC12A5, SYT1, GABRA1, and NOTCH1, among others, which are known as important signaling factors in neuron development and differentiation. Lastly, for the mesenchymal subtype, the expression signature includes YKL40, MET, CD44, MERTYK, TRADD, RELB, and TNFRSF1A (Verhaak et al. 2010). TNF-α, NF-κβ, STAT3, and CEBPB have been reported as related signaling pathways (Olmez et al. 2018; Yamini 2018; Tan et al. 2019). These signaling pathways are known to induce inflammatory responses. Mutation information has been reported for NF-κβ and NF1 (Verhaak et al. 2010).

Among these GBM subtypes, the mesenchymal subtype, characterized by traits associated with inflammatory responses within the tumor and increased vascularization, is associated with a higher correlation with tumor invasion and metastasis (Das and Marsden 2013; Zanotto-Filho et al. 2017). Treatment of the mesenchymal subtype is challenging due to its potential for malignancy (Behnan et al. 2019). The mesenchymal subtype of GBM is associated with various signaling pathways. Receptor tyrosine kinases (RTKs), which are part of these signaling pathways, are known to functionally regulate tumor cell growth, cell cycle progression, and other processes (Pearson and Regad 2017). RTK inhibitors are widely used in many cancer types. Notable uses include imatinib (Gleevec) in gastrointestinal stromal tumors, erlotinib in non-small cell lung cancer, lapatinib (Tykerb) in HER2-positive breast cancer, sorafenib (Nexavar) in hepatocellular and renal cell carcinoma, and sunitinib (Sutent) in renal cell carcinoma, gastric cancer, and neuroendocrine tumors.

According to recent reports, a clinical trial was attempted with the FDA-approved drug bevacizumab targeting VEGFR2 in patients with recurrent GBM. However, clinical trials have not yet demonstrated significant superiority of RTK-targeted therapy using RTK inhibitors alone or in combination (Qin et al. 2021). Therefore, continuous research is needed to explore drug combination strategies based on molecular and cellular biology mechanisms, considering the genetic characteristics within the tumor. At the in-silico level, integrating heterogeneous data, including transcriptomics, genomics, and proteomics, from currently available databases, such as The Cancer Genome Atlas (TCGA), the Cancer Dependency Map (DepMap), and Clinical Proteomic Tumor Analysis Consortium (CPTAC), can facilitate the proposal of novel therapeutic strategies (Liu et al. 2020; Li et al. 2022).

Here, we proposed and identified combinatorial druggable target genes based on the existing transcription-based ones, as well as their interaction with transcriptomic expressions. We collected publicly available data, including cell lines, transcriptome, genome, and drug sensitivity data at the cell line level. In particular, we analyzed the expression patterns of RTK-related genes closely associated with GBM and identified the characteristics of cellular signaling pathways interacting with these genes. Furthermore, we found out that the candidate genes were associated with Homologous Recombination Deficiency (HRD) signature using mutational signature analysis (Lee et al. 2018). Thus, we proposed more refined subtypes through integrated analysis of transcriptomics, genomics, proteomics, and other data at the transcription-based subtype level. This study provided new candidates in specific subtype of GBM and offered personalized treatment strategies for patients with GBM based on this evidence.

Materials and methods

Data collection

The following data was collected from a public database. (Figure 1). Dependency Map (DepMap, https://depmap.org/portal/): {RNA expression (n = 19,221), Mutation (n = 18,784), CNV(n = 25,368), CERES (n = 17,386), PRISM, and, Metadata}, Genotype-Tissue Expression (GTEx, https://www.gtexportal.org/home/): {RNA expression (n = 19,221)}, UCSC Xena (TCGA dataset, https://xenabrowser.net): {RNA expression (n = 20,531), Mutation (n = 40,544), CNV(n = 24,777), and, Metadata}, and Clinical Proteomic Tumor Analysis Consortium (CPTAC, https://proteomics.cancer.gov/data-portal): {protein expression (n = 11,141) and phospho-protein (n = 101,266)}.

Figure 1. Research scheme. Genomes, transcriptomes, and proteomes were downloaded from DepMap, GTEx, TCGA, and CPTAC for utilization as follows: (1) Data preprocessing, (2) Discovery of genes significantly interacting with RTK and related genes, (3) Classification of signal transduction system and genome, and (4) Drug screening.

Cell line level analysis

A gene effect dataset involving 20,252 genes across 24 glioblastoma cell lines was utilized. CERES, a measure of gene dependency, was employed through CRISPR-Cas9 screens provided by DepMap to assess their necessity for cell survival (https://github.com/cancerdatasci/ceres). Specifically, a score lower than −0.5 represents depletion in most cell lines, while a score lower than −1 represents strong killing. To narrow down candidate genes, we filtered genes focusing on those with CERES scores below −1 across all glioblastoma cell lines.

Patient level analysis

RNA sequencing data was obtained from TCGA-GBM, comprising 11 normal samples, 577 primary tumor samples, and 13 recurrent tumor samples. Integration was done with 40 randomly selected GTEx (dbGaP Accession phs000424.v8.p2) brain datasets to account for differences in sample numbers between tumor and normal samples for DEG analysis. The TCGA-GBM data was represented as in log2(x + 1) transformed RSEM normalized count. and the GTEx brain data as expected count data. Batch effects were removed by using ComBat-seq (Zhang et al. 2020) To investigate the potential relationship between 30 Cancer Gene Signatures (CGSs) and 63 well-known Receptor Tyrosine Kinases (RTKs) (Robinson et al. 2000).

The association between multi-omics data within the cell signaling pathway

The interaction between HSPA5 and FGFR1 was analyzed to investigate RNA-protein interactions through multi-omics analysis. Based on the expression levels of HSPA5 and FGFR1, four distinct groups were classified: co_up (both HSPA5 and FGFR1 up-regulation), co_down (both HSPA5 and FGFR1 down-regulation), HSPA5_up (only HSPA5 up-regulation), and FGFR1_up (only FGFR1 up-regulation). Gene expression (DEG) and protein expression (DEP) analyses were conducted to analyze the TCGA and CPTAC datasets. Genes that showed no expression across all samples were excluded. For the DEG analysis, the ‘limma’ package was utilized for the analysis (Ritchie et al. 2015), and for the CPTAC dataset, t-tests were employed to investigate activated pathways within each group (|FC| ≥ 1, p-value ≤ 0.05). Additionally, the ‘PhosR’ package (Kim et al. 2021), a part of R's Bioconductor suite, was used to analyze protein phosphorylation. This analysis was executed only when more than half of the replicates in at least one condition possessed the corresponding phosphosites. In instances of missing values for phosphosites, imputation was performed by sampling from the empirical normal distribution, constructed using the quantification values of phosphosites from the same condition.

Drug screening

The 32 glioblastoma cancer cell lines were categorized into four groups based on the expression levels of FGFR1 and HSPA5, using the expression data from CCLE (22Q2) for each cell line. We downloaded and utilized the ‘primary-screen-replicate-collapsed-logfold-change.csv’ file from the DepMap portal. We examined the expression patterns of 32 GBM cell lines based on the expression levels of HSPA5 and FGFR1, using expression data from CCLE (22Q2) for each cell line. Additionally, expression patterns were investigated based on mesenchymal and cell morphology. Drugs exhibiting statistically significant differences were also identified.

Statistical analysis

A proportional hazard model was used to identify genes that could potentially interact with RTK genes in GBM. Additionally, Pearson correlation analysis was conducted to examine the correlation between Clinically Significant Gene Sets (CSGs) and RTK genes, requiring a Pearson’s correlation coefficient of at least 0.4. Using the chi-square test, we conducted analysis on the association between the expression of HSPA5 and FGFR1 and the transcription-based subtypes (classical, neural, proneural, and mesenchymal). Additionally, Fisher's exact test was conducted when the frequency was less than 5. A significance level of p ≤ 0.05 was applied to all statistics, and R package version 4.1.2 was used.

The analysis pipeline and datasets are provided through https://github.com/Honglab-Research/SPde.

Results

FGFR1 and HSPA5 co-upregulation indicates poor prognosis

For this study, we collected data from publicly available databases, including DepMap, Genotype-Tissue Expression (GTEx), TCGA, and CPTAC (Material and Methods). Using multi-omics data, we conducted subtyping based on gene expression and performed integrated analysis of genomic, transcriptomic, and proteomic data to extract characteristics at the patient level. Finally, we analyzed drug responses according to subtype at the cell line level using PRISM database. (Figure 1).

We discovered 477 genes associated with survival across 24 GBM cell lines based on CERES dependency scores from the DepMap database (Supplementary Figure 1). Next, we analyzed differentially expressed genes (DEGs) between tumor (TGGA-GBM) and normal samples (Integrated data; TCGA-GBM normal and GTEx normal brain data, Material and Methods). We then conducted Cox analysis to confirm whether changes in the expression of these genes affected the survival of patients with GBM. As a result, we successfully identified 30 genes that are clinically significant in GBM, and we termed them Clinically Significant Gene Sets (CSGs). (Figure 2(A)).

Figure 2. Co-regulation of HSPA5 and FGFR1 and survival characteristics across GBM subtypes. (A) Forest plot with HRs. Only selected 30 Clinically Significant Gene Sets(CSGs) from Cox analysis are shown. The HRs are presented as the centers of the error bars. The error bars range from the lower to the upper 95% confidence limit. A positive association between gene expression and mortality rate is displayed in a pink color. A negative association is displayed in a blue color. (B) Correlation analysis of RTK and CSG. Color gradient indicates correlation coefficient(r). We restricted the color gradient range to – 0.4–0.4 for display.(C) Survival analysis of the entire data set compared to co_up in TCGA data. (D) Survival analysis of the entire data set compared to co_down in TCGA data. (E) Integrated analysis of the expression pattern between the transcription-based subtype and HSPA5 and FGFR1. Count indicates The number of patients corresponding to each subtype within each group. (F) Differences in survival analysis according to the HSPA5 and FGFR1 expression patterns in the mesenchymal and proneural subtypes.

We investigated the correlation between the CSG and 61 well-known RTK genes. Only HSPA5 and BUD31 were associated with prominent RTK genes in GBM. (Figure 2(B)). In particular, HSPA5 showed positive correlations with two GBM RTK genes, FGFR1 and VEGFA (R = 0.50 and 0.52, respectively). We categorized 590 TCGA-GBM patients into four groups: co_up (HSPA5 and FGFR1 upregulated), co_down (HSPA5 and FGFR1 downregulated), h_up (only HSPA5 upregulated), and f_up (only FGFR1 upregulated). Survival analysis was performed based on the expression patterns. The prognosis was worse for the co_up group, compared to other groups in TCGA GBM data, whereas the prognosis for the co_down group showed no difference (Figure 2(C and D), p = 0.0016 and p = 0.14, respectively). We validated our data using the Chinese Glioma Genome Atlas dataset. similarly, the results showed that the prognosis was worse when comparing the co_up group with other groups (Supplementary Figure 2, p < 0.0001). In addition, 49% (123/249) of cases classified as WHO Grade IV, the highest grade, were found in the co_up group (chi-square test, p-value = 1.634e-07, Supplementary Figure 3B).

We investigated the association between the classification based on expression (proneural, neural, classical, mesenchymal) in GBM and the expression patterns identified in this study. The co_up patient group exhibited the highest prevalence of the mesenchymal subtype, followed by the classical subtype. In survival analysis, the prognosis was worse in the co_up group (M_coup) compared to the other groups in the mesenchymal subtype; in the proneural subtype, the prognosis was better in the co_down group (M_codn) (Figure 2(E and F), p = 0.027 and p = 0.0014, respectively).

Cell signaling pathways resulting from the interaction between HSPA5 and FGFR1 at the transcriptome and proteome levels

Based on the four groups identified in our study, we investigated cell signaling pathways. In the co_up group, we observed the activation of the PI3 K/AKT pathway. (Supplementary Figure 3A and Figure 3). Furthermore, it consistently observed at the protein level with a high correlation. (Figure 3, r = 0.79)In addition, we identified consistent expression patterns between RNA and protein expression in genes related to cell adhesion, RTK, and other relevant genes (co_up vs co_down, p ≤ 0.05, Figure 3). Since cell adhesion is relevant in the mesenchymal subtype, we inferred that the co_up group identified in this study was affected downstream of the PI3K-AKT pathway. Additionally, we examined the activation of intracellular signaling pathways by investigating phosphoprotein activity. The co_up group contained a total of 18 regulatory protein modules regulated by four kinases. (Supplementary Figure 4A). These modules represent clusters where similar dynamic phosphorylation profiles and kinase regulations are shared among phosphorylation sites, delineating specific groups within a protein signaling network. (Ref). Notably, only two modules displayed significant differences between the co_up and co_down groups. Module 7, predominantly regulated by MAPK, exhibited enhanced activity in the co_down group, while Module 8, regulated by AKT, displayed increased activity in the co_up group. (Supplementary Figure 4B). Consequently, we discovered that the co_up group activates the AKT pathway and deactivates the MAPK pathway. Furthermore, utilizing transcriptome data from 693 CGGA, we confirmed that the upregulation of FGFR1 and H9,SPA5 leads to the activation of the PI3 K/AKT pathway. We also identified the downregulation of the MAPK pathway in the co_up group. (Supplementary Figure 3B).

Figure 3. Co-upregulation of HSPA5 and FGFR1 activates the PI3 K/AKT pathway. RNA and protein expression was schematized according to the PI3 K/AKT pathway. The expression patterns were classified as co_up, co_down, FGFR1_up, and HSPA5_up based on the expression patterns of FGFR1 and HSPA5, whose trends are displayed.

The effect of PI3 K inhibitor on subtypes separated by FGFR1 and HSPA5 gene expression patterns. We sought to validate our findings using the DepMap PRISM and CCLE datasets. Initially, we categorized 32 GBM cell lines into four groups based on FGFR1 and HSPA5 expression.(Supplementary Table 1)

Next, we compared the effect of ECM-related adhesion inhibitors between co_up and co_down group, targeting molecules upstream of the PI3 K/AKT pathway. However, PF-573228 (ECM-related adhesion inhibitor) showed higher drug sensitivity in the co_down group (Supplementary Figure 5A and Supplementary Table 2). Instead, we noticed that all PI3 K inhibitors exhibited higher drug sensitivity in the co_up group compared to the co_down group. Particularly, the co_up group exhibited high drug sensitivity when treated with PIK-93, in contrast to the co_down group (Figure 4(A), p = 0.027). Although not all MAPK inhibitors yielded similar results, SB-203580 exhibited significantly high drug sensitivity in the co_down group (Figure 4(A) and p = 0.009).

Figure 4. Verification of the response to PI3 K and MAPK inhibitors through the HSPA5 and FGFR1 interaction. (A) Box plots illustrate the extent of changes in cell viability when applying a compound within each group. A lower cell viability indicates higher drug sensitivity. The lower y-value indicates better drug sensitivity. The drug response differences to PIK-93 and SB-203580 were demonstrated according to the HSPA5 and FGFR1 expression patterns (p ≤ 0.05). (B) Drug response differences to PIK-93, CUDC-907, and XL-147 are shown according to the mesenchymal and proneural subtypes (p ≤ 0.05).

Compared to cell lines with the mesenchymal subtype, cell lines simultaneously categorized as co_up (M_coup) displayed increased sensitivity to PI3 K inhibitors (Figure 4(B) and XL-147: p = 0.006). Moreover, we confirmed that the M_coup group, in comparison to the mesenchymal subtype, exhibits higher sensitivity to PI3 K and HDAC inhibitors (Figure 4(B) and CUDC-907: p = 0.037). Based on the association between HSPA5 and FGFR1 interaction and the PI3K-AKT pathway, these results suggest an association with PI3 K inhibitors.

Genomic characteristics in subtypes separated by FGFR1 and HSPA5 gene expression patterns

We analyzed the relationship between the expression patterns of HSPA5 and FGFR1 and the genomic data (Figure 5(A)). As a result, in the Mesenchymal subtype, we observed prominent mutations of the NF1 gene in the co_up group, while in the Proneural subtype, Mutations of the IDH1 and TP53 were detected in the co_down group. Furthermore, the CGGA data showed the co_up group comprised 53% (153/286) of IDH Wild Type cases (chi-square test: p-value = 5.253e-09 and Supplementary Figure 3B). The results were consistent with our previous findings.

Figure 5. Co-regulation of HSPA5 and FGFR1 at the genomic level. (A) A oncogrid of the relationship between the HSPA5 and FGFR1 expression pattern, transcription-based subtype, clinical information, and genomic mutations. (B) Utilizing the mutation data within the genome, we analyzed the mutation signature resulting from the differences in HSPA5 and FGFR1 expression in the mesenchymal and proneural subtypes.

Additionally, we analyzed mutation signatures to examine the characteristics of mutation signatures between the transcriptome-based subtype analysis and the expression patterns of HSPA5 and FGFR1. In GBM overall, we identified several mutation signatures in the mutation signature analysis, including single base substitution SBS2 associated with APOBEC enzymes, SBS18 related to reactive oxygen species damage, SBS30 linked to base excision repair, SBS32 associated with azathioprine-induced immunosuppression, and SBS44 related to DNA mismatch repair. We discovered SBS19 and SBS40, whose functions are not yet known. Furthermore, we identified a common signature including SBS19 and its unknown function and base excision repair related SBS30 in the expression patterns of HSPA5 and FGFR1 within both the mesenchymal and proneural subtypes. Specifically, within the mesenchymal subtype, only SBS3 related to HRD signature was found in the co_up group.

Therefore, to explore the relationship between the newly discovered group (M_coup) and germline mutations in homologous recombination-related genes, we conducted chi-square tests on 58 HRD-related genes. Approximately 12% (6/52) of cases in the newly discovered group (M_coup) possessed ERCC4 germline mutations. (p-value = 0.055) (Supplementary Table 4). However, we could not find differentially expressed genes related to HRD (Supplementary Table 5).

Discussion

GBM is a type of cancer that shows a very poor prognosis and is very challenging to treat. GBM subtypes are divided based on transcriptional criteria, and treatment strategies tailored to each subtype are necessary. Despite the treatment strategies targeting RTKs proposed in the past, clinical reports have indicated their lack of effectiveness. Therefore, detailed classifications of GBM are necessary in addition to the existing subtypes.

In this study, we discovered poor GBM prognosis due to the interaction between HSPA5 and FGFR1 within the existing subtypes. Especially in the mesenchymal subtype, we suggested a potential direction for treatment strategies by presenting characteristics and mutation signature features according to the co-regulation and expression of HSPA5 and FGFR1.

Temozolomide is currently known as an important chemotherapeutic for GBM treatment, but its limitations are highlighted by drug resistance (Chien et al. 2021; Roh et al. 2023; Teraiya et al. 2023). Other treatment methods have been studied, including clinical trials for targeted therapies, such as bevacizumab, which targets VEGF (Fu et al. 2023) and EGFR (erlotinib, lapatinib, and gefitinib) inhibitors that focus on EGFR, and PI3 K/AKT/mTOR inhibitors aimed at the PI3 K/AKT/mTOR signaling pathway (Pan and Magge 2020). However, the clinical efficacy of treatments targeting RTK exhibits limitations, regardless of their single or combined use (Qin et al. 2021; Shao et al. 2023). While the PI3 K inhibitor has been approved by the FDA as an anticancer drug, its widespread clinical application is limited due to frequent and severe side effects (Lin et al. 2021; Yu et al. 2023).

Despite the classification of GBM subtypes based on existing transcriptional standards, the need for effective treatments persists. Hence, integrated analysis of multi-omics data encompassing transcriptomes, genomes, and proteomes is needed. Through our integrated analysis, we have proposed a new GBM subtype driven by the co-regulation of HSPA5 and FGFR1 in existing GBM subtypes. FGFR1 has shown genomic abnormalities not only in GBM, but in most types of cancer as well, and its increased expression has been reported in GBM (Morrison et al. 1994; Yamaguchi et al. 1994; Jimenez-Pascual and Siebzehnrubl 2019). HSPA5 is recognized for its role in managing oxidative stress protection and cell survival, and its elevation has been observed in a variety of cancers, including GBM (Iglesia et al. 2019).

Phillips et al. (2006), Verhaak et al. (2010), and Wang et al. (2017) have classified existing GBM subtypes through criteria, such as expression signatures, chromosome gain/loss, and mutated genes (Phillips et al. 2006; Verhaak et al. 2010; Wang et al. 2017). With the current ability to correlate RNA and protein expression through integrated multi-omics data analysis, classification criteria should be expanded based on the conclusions obtained from this study. Additionally, the mutation signature presented by Alexandrov et al. (2013) can serve as a method to describe the characteristics of mutation types occurring in specific mutation-inducing processes and to represent the classification features of GBM subtypes.

In this study, we discovered that the interaction between HSPA5 and FGFR1 regulates PI3 K/AKT cell signaling. Specifically, we found co-regulation of RNA and protein in both cell adhesion and RTK, which are highly relevant to the upstream mesenchymal side (Mathew et al. 2014,; Behnan et al. 2019). Therefore, regulating PI3 K/AKT in downstream signal transduction could enable new treatment strategies (Figure 3).

Additionally, mutation signature analysis showed only SBS3 in mesenchymal subtype cases where HSPA5 and FGFR1 were co-upregulated (Figure 5). SBS3 is related to HRD, a major DNA repair mechanism within tumors. Reports have linked HRD components to poor GBM prognosis (Knijnenburg et al. 2018; Lim et al. 2020). In cases of GBM with HRD, treatment strategies have been suggested that either utilize a PARP inhibitor or combine it with temozolomide (Berte et al. 2016; Bisht et al. 2022). Therefore, our study suggests the possibility of combined drug use including PI3 K and PARP inhibitors, along with existing treatments for GBM (Zhang et al. 2021). We wanted to assess the effects of a PARP inhibitor on the M_coup cell line, but there was no available data on the administration of this compound to the M_coup cell line (Supplementary Table 3).

Most drugs showing significant differences in drug sensitivity between the co_up and co_down groups target the PI3 K/AKT or MAPK pathways. We could identify other RTK inhibitors such as R-428 and TELATINIB (Supplementary Table 2 and Supplementary Figure 5B). However, due to limitations in in-silico based analysis, we could not compare which drug was more effective within the same cell.

We anticipate that the PI3 K inhibitor, consistently associated with the co_up group from transcriptomics to phosphoproteomics levels, would be most effective in that group. However, experimental validation is necessary to confirm this (Supplementary Figure 4B).

Therefore, in our subsequent study, we intend to utilize a mouse model to compare the efficacy of candidate drugs. Although the ERCC4 germline mutation was not significant in this study, the high proportion of approximately 12% and the low p-value (0.055) suggest that it may still serve as a promising molecular marker in M_coup group. (Supplementary Table 4). Therefore, we aim to collect samples from M_coup patients to assess the prevalence of ERCC4 germline mutations within these patient groups and investigate whether this mutation can serve as a molecular marker for the newly discovered subtype.

Furthermore, we anticipate identifying candidate genes capable of overcoming RTK drug resistance mechanisms in other cancer types. Consequently, we conducted analyses for 12 cancer types including LUAD, BRCA, COAD, PAAD, KIRC, STAD, ESCA, UCEC, LUSC, LIHC, and BLCA. Among these, significant correlations between RTKs and CSGs were observed in 8 cancer types. Based on these findings, we are planning the next stage of research (Supplementary Table 6).

We identified co-regulation of HSPA5 and FGFR1 in the existing transcription-based GBM classification, and we further refined the existing mesenchymal subtype based on their expression characteristics. In addition, our mutation signature proposes a more detailed classification at the genomic level. Therefore, this refined classification can improve the categorization of clinical treatment strategies and provide information on drugs that can be combined with existing treatments.

Acknowledgments

This work was supported in part by grants from the National Research Foundation of Korea (NRF) grant funded by the Korea government (Ministry of Science and ICT) (NRF-2021M3H9A2097227, NRF-2022R1A2C3008162, and RS-2023-00220840), the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI), funded by the Ministry of Health & Welfare, Republic of Korea (RS-2023-00265923), and the Basic Medical Science Facilitation Program through the Catholic Medical Center of the Catholic University of Korea funded by the Catholic Education Foundation. We thank the Global Science experimental Data hub Center (GSDC) and the Korea Research Environment Open NETwork (KREONET) service for data computing and network provided by the Korea Institute of Science and Technology Information (KISTI).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported in part by grants from the National Research Foundation of Korea (NRF) grant funded by the Korea government (Ministry of Science and ICT) (NRF-2021M3H9A2097227, NRF-2022R1A2C3008162, and RS-2023-00220840), the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI), funded by the Ministry of Health & Welfare, Republic of Korea (RS-2023-00265923), and the Basic Medical Science Facilitation Program through the Catholic Medical Center of the Catholic University of Korea funded by the Catholic Education Foundation. We thank the Global Science experimental Data hub Center (GSDC) and the Korea Research Environment Open NETwork (KREONET) service for data computing and network provided by the Korea Institute of Science and Technology Information (KISTI).