Overview of m6A transcriptome methylation within lncRNAs in GBM and normal brain tissues
To understand the pattern of the m6A modification profiles of lncRNAs in GBM, 3 human GBM tumor tissues versus 3 normal brain tissues were conducted for transcriptome m6A-seq, which revealed that a considerable fraction of lncRNAs were extensively m6A-modified. And, a total of 2113 m6A peaks were identified in the normal group, representing transcripts of 1538 genes, including 544 lincRNAs (intergenic lncRNAs), 720 antisense lncRNAs and 262 pseudogenes and 12 others; While, in GBM group, 2412 m6A peaks were identified, corresponding with the transcripts of 1508 genes, consist of 445 lincRNAs, 771 antisense lncRNAs and 283 pseudogenes and 9 others (Fig. 1A). Further analysis found that the majority of lncRNAs (74.5% of the m6A-methylated genes in the normal group, while 65.3% of the methylated genes in the GBM group) contained only one peak in both groups, while a relatively small number of genes contained two peaks,and very few lncRNAs contained three or more peaks (Fig. 1B). To further analyze the distribution profiles of m6A peaks within lncRNAs, we found that m6A site is distributed in intron and exon regions in the almost the same proportion in normal group; While, A little bit more inclined to distribute in exon regions in GBM tissues, compared to normal group (Fig. 1C).
To uncover the significance of m6A methylated lncRNAs in GBM, the differences and overlaps in m6A modified lncRNAs between the individuals were analyzed by the Venn diagram. As shown in Fig. 1D, 781 m6A-modified lncRNAs were detected within both groups. Compared with the normal group, the GBM group had 727 new genes, with disappearance of 757 genes, indicating a significant difference in global lncRNA m6A modification patterns between the GBM and normal brain groups. To explore the effect of m6A on lncRNA expression, the differentially expressed lncRNAs between 153 primary GBM and 5 normal brain tissues in TCGA GBM cohort were compared. Compared with normal samples, 4375 lncRNAs were differentially expressed in GBM tissues (LogFC ≥ 1 and padj ≤ 0.05), with 2614 upregulated genes and 1761 downregulated genes (Table S3). And, the global abundance of the m6A peaks between GBM and normal brain tissues were also compared. Further, conjoint analysis of m6A diff-modified lncRNAs and TCGA diff-regulated lncRNAs were conducted, and the distribution of genes with a significant change in both m6A level (|FC| ≥ 1.2, p ≤ 0.05) and overall transcript expression level (|FC| ≥ 2, padj ≤ 0.05) were shown in Fig. 1E, which were mainly divided into four groups, including 208 hypermethylated as well as upregulated genes (‘hyper-up’), 394 hypomethylated as well as downregulated genes (‘hypo-down’), 56 hypermethylated but downregulated genes (‘hyper-down’) and 64 hypomethylated but upregulated genes (‘hypo-up’) in GBM tissues relative to normal brain tissues. And we discovered a positive correlation of differentially methylated m6A peaks and their corresponding gene expression levels (Fig. 1E, pearson correlation r = 0.495, p < 0.0001), and the m6A modification sites of the above four group genes were all distributed in the exon region, reconfirming that m6A modification can regulate the expression of mature lncRNAs. These results showed that there was an obviously distinct m6A modification patterns between GBM and normal brain tissues. And m6A modification could regulate a large amount lncRNAs, may via regulating their stability, degradation or some other functions, having been reported previously.
Identification of immune-stromal-m 6 A-related pseudogene HSPA7 a novel risk prognostic factor in GBM
To investigate the effects of these m6A regulated lncRNAs on TME cell infiltration, we assessed tumorpurity, stromal and immune scores of 153 primary GBM cases in TCGA GBM cohort, using the ESTIMATE algorithm (see Methods, Table S4). Then we performed differential analysis of all RNA-seq data from 153 GBM samples in TCGA database, based on the median cutoff immune/stromal scores, respectively. And the volcano plot of the high/low scores of the stromal/immune scores revealed differential gene expression profiles between the samples, in which 791 up-regulated lncRNAs and 459 down-regulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were obtained based on the difference in immune scores (Fig. S1A, Table S5), simultaneously, 734 up-regulated lncRNAs and 269 down-regulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were obtained based on the differential analysis of stromal scores (Fig. S1B, Table S6), As can be seen from the venn diagram (Fig. S1C-D), there were 4 identical up-regulated genes and 11 identical down-regulated genes, related to immune, stromal and m6A. Then, we performed Kaplan–Meier curve analysis on the obtained 15 differential genes, only HSPA7 (Fig. 2B) and AC011899.9 (Fig. S2B) had prognostic significance (p ≤ 0.05) in 153 TCGA GBM samples based on the median expression. The GEPIA database showed that HSPA7 was significantly overexpressed in GBM tissues, compared with GETx normal brain tissues (Fig. 2A), while AC011899.9 was not (Fig.S2A). Thus, we focused only on HSPA7 in the following work. And the m6A methylation peak distribution and abundances in HSPA7 transcript of GBM and normal brain tissues, detected by m6A-seq, was visualized using IGV software (Fig. 2C). And We found that HSPA7 was highly enriched in the m6A precipitated fraction, and the m6A modification enrichment level could be regulated by methyltransferase like 3 (METTL3), containing catalytic activity domain to catalyze the m6A formation (Fig. 2D). Also, the expression of HSPA7 could be significantly inhibited by knocking down the METTL3 (Fig. 2E) and overexpressing the FTO, a m6A demethylases (Fig. 2F). These results provided evidence that HSPA7 harboured high m6A modification level, and thus could be regulated in an m6A dependent manner.
To explore the association between expression of HSPA7 and clinical characteristics, we first compared the expression in the TCGA GBM cohorts stratified by molecular subtypes, IDH1 status, G-CIMP status, MGMT promoter status, age and gender, separately. As shown in Fig. 2G, the HSPA7 expression in the proneural subtype were signaficantly lower than those in the classical and mesenchymal subtypes, and the expression in the mesenchymal samples was the highest. The expression was lower in IDH mutant samples than those with wild-type IDH1. For G-CIMP status, the expression was lower in patients with G-CIMP than those without G-CIMP. And there were no obvious changes in the characteristics of MGMT promoter methylation, age and gender. Furthermore, univariate cox regression analysis of overall survival of GBM samples within TCGA showed that high HSPA7 expression (HR:1.511, P = 0.056) was the independent risk factors associated with prognosis of GBM. And high HSPA7 expression (HR:1.418, p = 0.034) still remained statistically significance in GBM samples, after taking into account age, gender, IDH status, MGMT promoter methylation status and G-CIMP status in subsequent multivariate cox regression analysis (Fig. 2H). These results indicated that the pseudogene HSPA7 is a novel prognosis biomarker and indicate therapy outcome for GBM patients.
HSPA7 is correlated with immunophenotypes and TME landscapes.
To further gain insight into the exact role of HSPA7 in determining immunophenotype, we analyzed 29 immune-associated gene sets, representing diverse immune cell types, functions and pathways (see Methods). As shown in the heatmap (Fig. 3A), we found that HSPA7 high expression group had significantly greater TME cell infiltration characterization, higher immunescore, stromalscore and lower tumorpurity, confirming that HSPA7 could indeed regulate immune cell infiltration and immune related genes expression. And we then explored the specific difference of 16 immune cells between high and low HSPA7 expression patients. We found tumors with high expression of HSPA7 presented significantly increased infiltration in immunosuppression populations, such as macrophages, neutrophils and Tregs (regulatory T cells), compared to patients with low expression; While, some immunity activating cells also enriched, including activated DCs (aDCs), immature DCs (iDCs), plasmacytoid DCs (pDCs) and TILs (Tumor infiltrating lymphocytes) (Fig. S3A), indicating the complexity of TME, in which GBM cells elicited multiple biological behavior changes through direct or indirect interactions with other TME components. Macrophages, neutrophils and Tregs and immune checkpoints, suppressing cytotoxic T lymphocyte (CTL) function, are emerging as critical regulators of cancer immune escape, leading to an immunosuppressive microenvironment that is permissive to tumor outgrowth[32, 33]. And two recent single-cell mapping of human brain cancer literature[7, 8] synergistically reported that abundant cellular component of the brain TME are microglia and monocyte-derived macrophages, having both tumor-promoting and immunosuppressive functions. Moreover, they also found the elevation of neutrophil infiltration and relative frequencies of Tregs in tumor tissues. They also found that both activation and exhaustion of lymphocytes were prevalent. revealing the diversity of, and the dynamic relationships between cellular components, determining the brain TME complexity and multifaceted functions. And the Univariate and multivariate analysis also showed that the macrophage and neutrophils were risk factors in TCGA GBM patients (Table S7). We found that HSPA7 had a significantly positive correlation with the enrichment scores of the three immunosuppressive immune cells and immune checpoints (Fig. 3B). We then explored the specific difference of markers of suppressor function of myeloid lineages (CD33, NOS2, CD163, CD68, SPP1, CD14, CD206), immune inhibitory checkpoints (LAG3, CTLA4, HAVCR2, PDCD1, PDCD1LG1, CD27, TIGIT, CD274 and TNFRSF9) and major neutrophil-recruiting chemokines and their receptors (ITGA3, CD177, MET and CXCL8) between high and low HSPA7 expression patients. The results showed that tumors with high expression of HSPA7 presented significantly increased expression compared to patients with low expression (Fig. 3C-E), and there were also positive correlations between HSPA7 and these molecules (Fig. S4A-C), indicating that HSPA7 could restrain the intra-tumoral antitumor immune response.
To further understand the exact role of HSPA7 in determining the TME profile and immunophenotype, we performed unsupervised consensus clustering based on the TME cell populations and immune related function genesets, developing a TME pattern with two clusters, TME immune low phenotype (Immune-L) and high (Immune-H) (Fig. S3B). And compared to the Immune-L group (102 patients), patients in Immune-H (51 patients) experienced the outcome of poorer prognosis (Fig. 3F) and higher HSPA7 expression (Fig. 3G). The PN (proneural), CL (classical) and MES (mesenchymal) subtypes of GBM have been most consistently described in the literature, with PN relating to a more favorable outcome and MES relating to poorer survival. And the association between MES gene expression signature, characterized by NF1 mutant, with reduced tumor purity, elevated invasion, migration capacity and infiltration of immunosuppressive cells (macrophages or microglia, mesenchymal stem cells or others), has been identified as a common theme across cancer[33, 34]. Our results showed that majority of MES GBM samples were high expression of HSPA7, as well as high immune-infiltrating phenotypes (Fig. 3H-I). Furthermore, we found tumors with high expression of HSPA7 presented significantly increased expression in MES phenotype signature, compared to patients with low expression (Fig. S4D). Overall, HSPA7 may function as a robust indicator of immunophenotype, significantly correlated with a poorer immune response in GBM.
HSPA7 is correlated with stromal and carcinogenic activation pathways.
To explore the biological behaviors among these distinct HSPA7 expression samples, we used GSVA algorithm to estimate pathway enrichment scores for each sample (see Methods). Compared to low HSPA7 expression group, high expression group were more markedly enriched in stromal activation pathways (angiogenesis, epithelial mesenchymal transition, VEGF and TGF beta signaling pathways), carcinogenic signaling pathways (hypoxia, apoptosis, PI3K Akt MTOR signaling, glycolysis, KRAS signaling and some others) and immune responses (IL6 Jak stat3 signaling, playing immunosuppressive effects on T cell function and mediating ICB resistance in cancers, inflammatory response, complement and some others). Whereas, they were less enriched in DNA replication, DNA damage response and mismatch repair (MMR) related functions (Fig. 4A). Previous studies demonstrated that the activation of stroma cells in TME can induce the T-cell suppression and oncogenic signaling pathway via a complicated cellular and biological reconfiguration mechanisms[32, 36], attenuating tumor response to PD-L1 blockade. As shown in Fig. 4B, stromal activation pathways were positively correlated with carcinogenic signaling pathways and immune suppression pathways, while negatively correlated with DNA damage response and repair pathways in TCGA GBM samples. The further analyses for the activity of stroma-related pathways indicated high HSPA7 expression were significantly associated with higher stromal activation signatures (Fig. 4C), constructed by Mariathasan et al., positively correlated with fibroblasts enrichment (Fig. 4D), calculated by MCP counter, and positively correlated with the typical stromal cell activation related pathways (Fig. 4E). Furthermore, TGF-β is a pleiotropic cytokine associated with poor prognosis in GBM, playing a pro-tumorigenic role by promoting immunosuppression, angiogenesis, metastasis, mesenchymal transition, fibroblast activation[27, 38–40] and has been proved to promote the progression of GBM in an autocrine signaling loop. In our data, compared to low HSPA7 expression group, TGF-β ligand (TGFB1), and a TGF-β receptor (TGFBR2), two key regulators in stromal activation and EMT pathways, and other genes encoding ECM and matricellular proteins (COLA4A1, COLA4A2, ECM1 and FN1), all showed increased expression in high HSPA7 group (Fig. 4F). Previous studies report that ECM can be shaped by monocyte differentiated macrophages via over-expression of the cathepsin proteases CTSB and CTSW, which also showed high expression in HSPA7 high group (Fig. 4F). These results suggested that immune cells and stromal cells in the tumor microenvironment can work together to regulate the immunosuppressive microenvironment synergistically, thereby promoting tumor escape.
Then, we explored the genes, positively correlated with HSPA7 (pearson r ≥ 0.3, p ≤ 0.05, Table S8), enriched pathways via metascape database (see Methods). And, they were significantly enriched in immune response, including cytokine-mediated signaling pathway, leukocyte migration, differentiation and some other immune related pathways. Additionally, they can also enriched in the pathways involving the stromal activation, such as extracellular structure organization and response to wounding. These enriched terms interacted with each other to form a protein–protein interaction network (Fig. 5A). Further, enrichment analysis in PaGenBase showed that the genes were almost specifically expressed in spleen, blood, bone marrow and some other peripheral immune cell aggregation tissues (Fig. 5B), suggesting that HSPA7 could indeed promote the infiltration of immune cells into tumor tissues. And the enrichment analysis in TRRUST, a TF-target interaction database, showed that the genes were substantially transcriptionally regulated by NFκB1 and RELA, two important transcription factors involving in the NF-κB signalling pathways (Fig. 5C). Further GSEA analysis,using GO and KEGG database, also showed that HSPA7 positively correlated with immune response and extracellular structure organization (Fig. 5D). This demonstrated again that the HSPA7 was of great significance in shaping different TME landscapes.
The functions of HSPA7 was verified in two CGGA cohorts
To further validate the function of HSAP7 in GBM, we explored the expression pattern in two CGGA RNA-seq cohorts. We found that the expression of HSAP7 was highest in GBM (WHO IV), among the three different WHO grade glioma specimens (Fig. S5A, Fig. S6A). Kaplan-Meier survival curve analyses showed that GBM patients with high HSPA7 expression had poor survival outcome in both cohorts (Fig. S5B, Fig. S6B). Further GSVA enrichment analysis also showed that HSPA7 expression correlated with immunophenotypes, stromal activation and oncogenic pathways in GBM samples (Fig. S5C, Fig. S6C). These results confirmed that HSPA7 could be an independent risk prognosis biomarker for GBM patients. Moreover, it can regulate the TME immune response and stromal activation, promoting in the malignant progression of GBM tumors.
WTAP potentially regulate recruitment of the m 6 A methyltransferase complex to HSAP7 directly
To explore the potential molecular mechanism in regulating HSPA7 RNA metabolism, we first screen HSPA7-interacting proteins via NCBI database (Table S9). We verified WTAP, a regulatory subunit of the m6A methyltransferase, required for their localization into nuclear speckles, and enhancing the RNA-binding capability of METTL3 to target RNA[31, 45], as putative HSPA7-binding regulator. To better map the association between WTAP and HSPA7, we explored the expression pattern in different grade of glioma in TCGA, CGGA, Gravendeel and Rembrandt datasets. As shown in Fig. S7A, the expression of WTAP is highest in WHO IV (GBM) group, among the three WHO grades glioma and normal brain samples. Further, the overexpression of WTAP correlated with the poor overall survival of GBM patients in all four datasets (Fig. S7B). What’s more, the expression of WTAP was postively correlated with the expression of HSPA7 in TCGA GBM dataset (calculated by GEPIA database) and two CGGA GBM datasets (Fig. S7C). Further GSVA enrichment analysis showed that compared to WTAP low expression group, WTAP high group also showed higher TME immune cells infiltration, stromal activation and carcinogenic pathways activation, with higher immunescore, stromalscore and lower tumor purity in TCGA GBM cohort (Fig. S7D), CGGA GBM cohort 1 (Fig. S8A) and CGGA GBM cohort 2 (Fig. S8B), same as HSPA7. From above, we could speculate that m6A methylation modification HSPA7 was regulated by WTAP directly, recruiting the m6A methyltransferase complex. The conclusion and functional mechanisms need to be further demonstrated in the future work.
The following pathway and process enrichment analysis of HSPA7-interacting proteins showed that they can invovled in many biological behaviors, such as mitotic cell cycle phase transition, Neddylation, proteasome-mediated ubiquitin-dependent protein catabolic, HIF1α pathway and some others (Fig. S9A, Table S9). Then, we further expolred the enriched terms across the RBPs (RNA binding proteins), detected by CLIP-seq data dopoisted in stabase database. We found that they were mainly enriched in mRNA processing, nucleic acid transport, RNA catabolic process and so on (Fig. S9B, Table S10). And we can see component annotation in NCBI database that HSPA7 can be expressed in various regions of the cell (cytoplasm, nucleus, plasma membrane and nucleus), indicating that it may be involved in multiple cellular processes, and the exact mechenism of HSPA7 regulating the GBM TME need us spend a lot of effort and experimentation to clarify in the future.
HSPA7 holds promise in predicting therapeutic response to PD-L1 blockade
Previous studies indicated that GBM patients is resistant to immune checkpoint inhibitors (ICIs), due to its low mutation rate, PTEN-deficient immunosuppressive microenvironment, infiltration by myeloid-derived suppressor cells, and activation of tumor stromal cells [46, 47]. Then, we analyzed the distribution differences of somatic mutation between low and high HSPA7 expression in TCGA-GBM cohort using maftools package, and we found that the mutant rate of PTEN was increased signaficantly in HSPA7 high expression group, compared to the low expression (Fig. 6A, Low: 20%; High: 36%). Additionally, the NF1 mutant rate was also elevated obviously (Fig. 6A, Low: 8%; High: 11%). Chen et al.find that PTEN deficiency in GBM increases macrophage infiltration via upregulating lysyl oxidase (LOX) expression, in turn secreting SPP1 to support GBM survival. And NF1 mutant is a marker of mesentromal GBM subtype, which also showed that high macrophage infiltration and activation. And then we analysised the association between HSPA7 and macrophage assciated genes ( LOX, SPP1, MACRO, the maker of SPP1 macrophages, and two important chemokines, CCL2 and CCR2). We found that HSPA7 had a postive correlation with these genes significantly (Fig. S10). What‘s more, Robert M Samstein et. found that TMB predicted survival after immunotherapy across multiple cancer types genomic (including GBM), and we downloaded the mutantion data of 82 GBM samples, and then analyzed the distribution differences of somatic mutation between bottom 30 low TMB and top 30 high TMB samples. We found that the mutant rate of PTEN and NF1 increased signaficantly in low TMB samples, compared to the high group, same results as changes regulated by HSPA7 (Fig. S11A, PTEN: Low, 53%, High: 30%; NF1: Low, 20%, High: 7%). We further explored the association between HSPA7 expression and TMB. As shown in Fig. S11B, HSPA7 was statically negative with TMB. And Touat M et. recently found that mismatch repair (MMR)-deficient gliomas were characterized by a lack of prominent T cell infiltrates, extensive intratumoral heterogeneity, poor patient survival and a low rate of response to PD-1 blockade. In our data, we also found that high HSPA7 expression group showed the reduced MMR associated pathways (Fig. 4A), which nagatively correlated with the immunesupression, stromal activation and oncogene pathways (Fig. 4B). Some other reserches also showed the same results[27, 51]. Furthermore, HSPA7 is nagetively correlated with MMR genes, such as MLH1, MSH2, MSH6 and PMS2 (Fig. S11C). These results promoted us that HSPA7 maybe a potential marker for ICB thrrapy.
To further investigate the predicapacity of immune checkpoint therapy of the HSPA7 in GBM. However, there was few published immunotherapy datasets received ICB thrapy. Then we used the melanoma cohort that received anti-CTLA4 therapy and urothelial cancer cohorts that received anti-PD-1 therapy (IMvigor210) to test the immunotherapy response of HSPA7 as a complement. We found that patients with low HSPA7 expression exhibited significantly clinical benefits in anti-CTLA4 cohort (Fig. 6B, anti-CTLA cohort, HR 1.927(0.8997–4.129)). Further, compared to HSPA7 high expression group, we found the significant therapeutic benefit and clinical response to anti-CTLA-4 immunotherapy in those with low HSPA7 group (Fig. 6C-D). Same result were also found in the anti-PD-L1 cohort, Kaplan–Meier curve analysis showed that the bottom 50 low HSPA7 expression patients exhibited markedly prolonged survival, compared to those in the top 50 HSPA7 expression group (Fig. 6E, HR 1.670 (1.022–2.792)). Also, there was a significant better therapeutic and clinical response to anti-PD-L1 immunotherapy in HSPA7 low expression group, compared to HSPA7 high group (Fig. 6F-G).
Further GSVA enrichment analysis showed that compared to bottom low HSPA7 expression samples, HSPA7 high group also showed higher TME immune cells infiltration, with higher immunescore, stromalscore and lower tumor purity (Fig. S12A), and higher stromal activation, higher carcinogenic pathways activation, as well as lower MRR pathways (Fig. S12B), consistent exactly with those enrienched functions in GBM (Fig. 3A, Fig. 4A, Fig. S6C and Fig. S7C). The PD-L1 expression of immune cells (IC) and tumor cells (TC) was also detected in the IMvigor210 cohort, we then explored the the difference of HSPA7 expression among different expression level groups. As we can see in Fig. 6H-I, patients with higher PD-L1 expression levels, both IC and TC, were observed with higher HSPA7 expression ( anova suumary, IC: p = 0.0001; TC: p = 0.0009), indicating that HSPA7 could upregulated the PD-L1 expression, suppressing the immune activation. Moreover, previous studies indicated that F-box and WD repeat domain containing 7 (FBXW7) is a vital tumor suppressor in varius cancers, controling proteasome-mediated degradation of oncoproteins such as cyclin E, c-Myc, Mcl-1, mTOR, Jun, Notch and so on, and its loss-of-function mutation promotes resistance to anti-PD-1 through the downregulation of viral sensing pathways. We found that compared to the FBXW7 wildtype group, HSPA7 experssion showed a obvious decrease in mutant group (Fig. 6J). Fibroblast growth factor receptor 1 (FGFR1) is frequently mutanted in various tumors, whose inhibitors have shown promising therapeutic value in several preclinical models. Palakurthi S et. found that the combination of FGFR inhibition and PD-1 blockade can promote tumor-intrinsic induction of antitumor immunity. And we also found that HSPA7 experssion is siganificantly higher in mutant group, compared to wildtype group (Fig. 6K). Taken together, our work strongly indicated that HSPA7 expression would contribute to predicting the response to immune checkpiont therapy.
We investigated whether HSPA7 expression was correlated with prognosis in pan-cancer patients. The expression of HSPA7 in tumor tissues and GETx normal brain tissues were calculated via GEPIA database. As shown in Fig. S13A, the expression of HSPA7 were statistically lower in ACC, COAD, DLBC, LAML, LUAD, LUSC, READ, THCA and THYM, while signaficantly higher in GBM, KIRC, KIRP and PAAD, compared to the corresponding normal brain tissues. There are also genetic alteration in many other tumors, although not significant. And the COX survival rates and Kaplan–Meier curve analysis were evaluated using the sangerbox tool. The relationships between HSPA7 expression and prognosis of different tumors were shown in Fig. S14-16. Statically, HSPA7 was a risk factor in most cancers (such as ACC, LDBC, CESC and so on). Moreover, we found that it exhibited relatively favourable factor in several cancer types, such as SKCM, KIRP, SARC and so on. We next evoluted the immunescore and stromalscore across 33 cancers from TCGA database. We found that HSPA7 had significantly postive correlations with the immunescore (Fig. S17) and stromalscore (Fig. S18) among almost all 33 cancer types. Moreover, HSPA7 had significantly postive correlations with immune checkpoints expression (Fig. S19A) and immune response (Fig. S19B) in most cancer types, such as GBM, LGG, OV, LUAD and some others. However, there was no relevance in some cancers, such as SKCM, UVM, TGCT, PAAD and so on, indicating that HSPA7 may functioned in other manner other than regulating TME in these cancers, and the molecular mechanisms of this gene needs to be studied in the specific tumor in the future. Together these data highlight HSPA7 may be a key factor, assisting various cancers acquired various immunophenotype, and have potentials to be as a consideration in the development of more-effective immunotherapies for GBM and other immunotherapy-resistant tumors. Because different immune cells may play different roles in different tumours, functions of HSPA7 under specific tumor are necessary to explore in the future.