The N6-Methyladenosine-Modified Pseudogene HSPA7 Correlates With the Tumor Microenvironment and Predicts the Response to Immune Checkpoint Therapy in Glioblastoma

Background Glioblastoma (GBM), one of the most aggressive tumors of the brain, has no effective or sufficient therapies. Identifying robust biomarkers for the response to immune checkpoint blockade (ICB) therapy, a promising treatment option for GBM patients, is urgently needed. Methods We comprehensively evaluated lncRNA m6A modification patterns in m6A-sequencing (m6A-seq) data for GBM tissues and systematically investigated the immune and stromal regulators of these m6A-regulated lncRNAs. We used the single-sample gene-set enrichment analysis (ssGSEA) algorithm to investigate the difference in enriched tumor microenvironment (TME) infiltrating cells and the functional annotation of HSPA7 in individual GBM samples. Further, we validated that HSPA7 promoted the recruitment of macrophages into GBM TME in vitro, as well as in our GBM tissue section. We also explored its impact on the efficacy of ICB therapy using the patient-derived glioblastoma organoid (GBO) model. Results Here, we depicted the first transcriptome-wide m6A methylation profile of lncRNAs in GBM, revealing highly distinct lncRNA m6A modification patterns compared to those in normal brain tissues. We identified the m6A-modified pseudogene HSPA7 as a novel prognostic risk factor in GBM patients, with crucial roles in immunophenotype determination, stromal activation, and carcinogenic pathway activation. We confirmed that HSPA7 promoted macrophage infiltration and SPP1 expression via upregulating the YAP1 and LOX expression of glioblastoma stem cells (GSCs) in vitro and in our clinical GBM tumor samples. We also confirmed that knockdown of HSPA7 might increase the efficiency of anti-PD1 therapy utilizing the GBO model, highlighting its potential as a novel target for immunotherapy. Conclusions Our results indicated that HSPA7 could be a novel immunotherapy target for GBM patients.


Methods:
We comprehensively evaluated lncRNA m 6 A modification patterns in m 6 Asequencing (m 6 A-seq) data for GBM tissues and systematically investigated the immune and stromal regulators of these m 6 A-regulated lncRNAs. We used the single-sample gene-set enrichment analysis (ssGSEA) algorithm to investigate the difference in enriched tumor microenvironment (TME) infiltrating cells and the functional annotation of HSPA7 in individual GBM samples. Further, we validated that HSPA7 promoted the recruitment of macrophages into GBM TME in vitro, as well as in our GBM tissue section. We also explored its impact on the efficacy of ICB therapy using the patient-derived glioblastoma organoid (GBO) model.

Results:
Here, we depicted the first transcriptome-wide m 6 A methylation profile of lncRNAs in GBM, revealing highly distinct lncRNA m 6 A modification patterns compared to those in normal brain tissues. We identified the m 6 A-modified pseudogene HSPA7 as a novel prognostic risk factor in GBM patients, with crucial roles in immunophenotype determination, stromal activation, and carcinogenic pathway activation. We confirmed that HSPA7 promoted macrophage infiltration and SPP1 expression via upregulating the YAP1 and LOX expression of glioblastoma stem cells (GSCs) in vitro and in our clinical GBM tumor samples. We also confirmed that knockdown of HSPA7 might increase the efficiency of anti-PD1 therapy utilizing the GBO model, highlighting its potential as a novel target for immunotherapy.

INTRODUCTION
Glioblastoma (GBM), one of the most aggressive brain tumors, currently has no effective and sufficient therapies due to its intratumoral heterogeneity and molecular complexity. Immune checkpoint blockade (ICB) therapy is being actively pursued as a promising treatment option for GBM. However, very few patients respond to this therapy (1)(2)(3)(4)(5) partly because of the contribution of prominent immunosuppressive factors in the brain tumor microenvironment (TME) in GBM, including tumor-associated macrophages (TAMs), neutrophils, and regulatory T cells (Tregs) (6)(7)(8). The development of biomarkers and identification of the definite molecular mechanism underlying resistance to ICB therapy are urgently needed to identify effective therapeutic strategies for GBM. N 6 -Methyladenosine (m 6 A), the most abundant reversible methylation modification of mRNA, critically affects processes in mRNA metabolism, including splicing, export, translation, and decay. Dysregulation of this modification is clearly linked to diverse pathological processes and disease progression (9,10), including GBM tumorigenesis (11)(12)(13)(14)(15). Recent studies have described the role of m 6 A modification in regulating the immune response (16)(17)(18)(19), prompting us to reveal the importance of the spectrum of m 6 A-regulated genes and m 6 A regulatory mechanisms in shaping the TME. Numerous studies have demonstrated that m 6 A is also present in numerous long non-coding RNAs (15,(20)(21)(22) (lncRNAs, transcripts longer than 200 nucleotides but lacking functional coding capacity), and lncRNA m 6 A modification has emerged as a fundamental player in cancer progression and immune regulation, suggesting a potential association between the tumor immune response and m 6 A lncRNA modification. However, the lncRNA m 6 A methylation profile has not been systematically clarified in GBM tumors. Additionally, none of these studies have specifically investigated the roles of m 6 A-modified lncRNAs in the overall TME landscape in GBM. Therefore, it is worthwhile to obtain comprehensive knowledge of the cellular TME infiltration characteristics mediated by m 6 A-modified lncRNAs, as this knowledge would contribute to our understanding of the role of m 6 A modification in immune regulation and guide the development of more effective immunotherapeutic strategies.
Here, we depicted the first transcriptome-wide lncRNA m 6 A methylation profile in GBM and normal brain tissues via m 6 A sequencing (m 6 A-seq) data analysis, revealing the highly distinct lncRNA m 6 A modification patterns between these two groups. Key immune-stromal-related lncRNAs were identified by differentially expressed gene (DEG) analysis in primary GBM cohorts from The Cancer Genome Atlas (TCGA). Integrating the m 6 A-regulated lncRNAs revealed in this study, we identified HSPA7 as a novel prognostic factor in GBM patients. Through detailed bioinformatic analyses of HSPA7, we identified its crucial role in immunophenotype determination, stromal activation and carcinogenic pathway activation and highlighted its robust capacity to predict the ICB response. We confirmed that HSPA7 facilitated macrophage infiltration via the YAP1-LOX axis in vitro. We also confirmed that HSPA7 enhanced the efficiency of anti-PD1 therapy utilizing GBM patient-derived glioblastoma organoids, an ex vivo model. These results demonstrated that HSPA7 could be a novel immunotherapy target for GBM patients.

Patient Specimens and Public Patient Cohorts
Human GBM and normal brain tissues for m 6 A-seq were obtained from patients admitted to Qilu Hospital. All participants provided written informed consent, and the research was approved by the Ethics Committee on Scientific Research of Shandong University Qilu Hospital (approval number: KYLL-2018-324).
The RNA sequencing (RNA-seq) transcriptome, somatic mutation data, and corresponding clinicopathological parameters of the TCGA GBM cohort were obtained from the TCGA database (http://cancergenome.nih.gov/). Two Chinese Glioma Genome Atlas (CGGA) GBM RNA-seq datasets and the corresponding clinicopathological parameters were obtained from the CGGA database (http://www.cgga.org.cn/). For ICB data, genomic and clinical information from the IMvigor210 cohort, complete expression data, and detailed clinical annotations were obtained from http://research-pub.Gene.com/ imvigor210corebiologies based on the Creative Commons Attribution 3.0 license. For data from GBM patients treated with PD-1 inhibitors (pembrolizumab or nivolumab), clinical information was obtained from Supplementary paper data, and the sequencing data were obtained from SRA PRJNA482620 (2). For patients with melanoma treated with anti-CTLA4 therapy, the expression data were downloaded from the cBioPortal database (http://www.cbioportal.org/), and the detailed clinical characteristics of individual patients were obtained from the supplementary data of a previous paper (23). In addition, the somatic mutation data and information in Figure S11A were obtained from the cBioPortal database. The m 6 A-seq sequencing data have been deposited in SRA PRJNA661159 (the data are being processed, submission ID: SUB8069560, released when the paper is published). The processed data are available from the corresponding author upon reasonable request.

Estimation of TME Cell Characterization
We used the single-sample gene set enrichment analysis (ssGSEA) algorithm to calculate the enrichment score of immune cell infiltration into the GBM TME for each sample. Immune cell-related genes were obtained from Bindea et al. (24) and Robert L (25). (Table S1) and included genes related to immune cell types, immune-related pathways and functions. Based on the ssGSEA results, samples from the TCGA GBM cohort were classified into the high immune cell infiltration (immune-H) group or low immune cell infiltration (immune-L) group by using the "hclust" R package.

Functional Annotation and Pathway Enrichment Analysis
To explore the differences in biological behavior among the samples with distinct HSPA7 expression levels, we used selected HALLMARK (26) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (27) from gene sets from the Molecular Signatures Database (MSigDB) and other commonly used gene signatures (28)(29)(30) to estimate pathway enrichment scores for each sample (Table S2) by gene set variation analysis (GSVA) using the "GSVA" R package.

Transwell Assay
Transwell assays were performed in 24-well multiwell insert systems according to the protocol of the manufacturer. THP-1 cells were incubated with 100 ng/ml PMA (Sigma-Aldrich) for 24 h in vitro to induce their differentiation into macrophages and then added to the top chamber in serum-free media. The bottom chamber was filled with 10% FBS 1640 and GSC growth media.
After 24 h of incubation, the top chamber cells were removed using a cotton swab, and the membrane was fixed in 4% paraformaldehyde for 15 min and stained with crystal violet for 15 min. Five fields of adherent cells in each well were photographed randomly.

Statistical Analysis
Kaplan-Meier survival analysis was performed using GraphPad Prism 7.04, and significant differences between two groups were compared by the log-rank (Mantel-Cox) test. The waterfall function in the "maftools" package was used to visualize the mutational landscape in patients in the high-and low-HSPA7 groups or the high-and low-tumor mutation burden (TMB) groups. Student's t-test was used for two-group comparisons. For comparisons among more than two groups, the Wilcoxon test and one-way ANOVA were used for non-parametric and parametric data (32). P > 0.05 was considered nonsignificance (ns), P ≤ 0.05 was considered statisticallysignificant (*P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001). All data processing with R packages was performed using R Studio (version 3.6.3)

Overview of Transcriptome-Wide m 6 A Methylation in lncRNAs in GBM and Normal Brain Tissues
To understand the pattern of the lncRNA m 6 A modification profiles in GBM, three human GBM tumor tissues and three normal brain tissues were subjected to transcriptome m 6 A-seq, which revealed that a considerable proportion of lncRNAs were extensively m 6 A-modified. In addition, 2,113 m 6 A peaks were identified in the normal group, corresponding to the transcripts of 1,538 genes, including 544 long intergenic ncRNAs (lincRNAs), 720 antisense lncRNAs, 262 pseudogenes, and 12 others. In the GBM group, 2,412 m 6 A peaks were identified, corresponding to the transcripts of 1,508 genes, namely, 445 lincRNAs, 771 antisense lncRNAs, 283 pseudogenes, and 9 other genes ( Figure 1A). Further analysis showed that most of the lncRNAs (74.5% of the m 6 A-methylated genes in the normal group but 65.3% of the methylated genes in the GBM group) in both groups contained only one peak, while a relatively small number of lncRNAs contained two peaks, and very few lncRNAs contained three or more peaks ( Figure 1B). Upon further analysis of the distribution profiles of m 6 A peaks within lncRNAs, we found that m 6 A sites were distributed in almost the same proportion in intronic and exonic regions in the normal group but exhibited a slightly increased tendency to be distributed in exonic regions in GBM tissues ( Figure 1C).
To reveal the significance of m 6 A-methylated lncRNAs in GBM, the differences and overlaps in the m 6 A-modified lncRNAs between the GBM and normal brain groups were analyzed by constructing a Venn diagram. As shown in Figure 1D, 781 m 6 A-modified lncRNAs were common to both groups. A total of 727 new genes were expressed, and the expression of 757 genes was lost in the GBM group compared with the normal group, indicating a significant difference in global lncRNA m 6 A modification patterns between the GBM and normal groups. To explore the effect of m 6 A on lncRNA expression, the differentially expressed lncRNAs between 153 primary GBM and five normal brain tissues in the TCGA GBM cohort were compared. Compared with normal samples, GBM tissues exhibited 4,375 differentially expressed lncRNAs (LogFC ≥ 1 and padj ≤ 0.05), with 2,614 upregulated and 1,761 downregulated (Table S3). In addition, the global abundance of m 6 A peaks between GBM and normal brain tissues was also compared. Furthermore, integrated analysis of these differentially m 6 A-modified lncRNAs and differentially regulated lncRNAs from the TCGA dataset was conducted. The distribution of genes with a significant change in both the m 6 A level (|FC| ≥ 1.2, P ≤ 0.05) and the overall transcript expression level (|FC| ≥ 2, padj ≤ 0.05) is shown in Figure 1E. These genes were divided into four main groups: 208 were hypermethylated and upregulated ("hyperup"), 394 were hypomethylated and downregulated ("hypodown"), 56 were hypermethylated but downregulated ("hyperdown") and 64 were hypomethylated but upregulated ("hypoup") in GBM tissues relative to normal brain tissues. We also discovered a positive correlation between differentially methylated m 6 A peaks and the expression levels of their corresponding genes [ Figure 1E, Pearson correlation coefficient (r) = 0.495, P<0.0001], and the m 6 A modification sites in all of the above four groups of genes were distributed in exonic regions, reconfirming that m 6 A modification can regulate the expression of mature lncRNAs. These results revealed obviously distinct m 6 A modification patterns between GBM and normal brain tissues. Moreover, m 6 A modification could regulate numerous lncRNAs, possibly by regulating their stability, degradation or other functions, as reported previously.

Identification of the Immune-Stromal-m 6 A-Related Pseudogene HSPA7 as a Novel Prognostic Risk Factor in GBM
To investigate the effects of these m 6 A-regulated lncRNAs on cell infiltration into the TME, we assessed the tumor purity, stromal score, and immune score of 153 primary GBM cases in the TCGA GBM cohort using the ESTIMATE algorithm (see Materials and Methods, Table S4). Then, we performed differential analysis of all RNA-seq data from these 153 GBM samples in the TCGA database based on the median cutoff immune/stromal scores. The volcano plot of the high/low stromal/immune scores revealed differential gene expression profiles between the samples. A total of 791 upregulated lncRNAs and 459 downregulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were identified based on the difference in immune scores (Supplementary Figure S1A and Table S5). Simultaneously, 734 upregulated lncRNAs and 269 downregulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were identified based on the differential analysis of stromal scores ( Figure S1B and Table S6). As the Venn diagram ( Supplementary Figures S1C, D) indicates, four identical upregulated genes and 11 identical downregulated genes were related to immune activation, stromal activation, and m 6 A modification. Then, we performed Kaplan-Meier analysis on patients stratified by the expression levels of these 15 differentially expressed genes, and only HSPA7 ( Figure 2B) and AC011899.9 (Supplementary Figure S2B) were found to have prognostic significance (P ≤ 0.05) in the 153 TCGA GBM samples based on the median expression levels. The Gene Expression Profiling Interactive Analysis (GEPIA) database was used to show that HSPA7 was significantly overexpressed in GBM tissues compared with normal brain tissues in the Genotype-Tissue Expression (GTEx) database (Figure 2A), while AC011899.9 was not (Supplementary Figure S2A). Thus, we focused only on HSPA7 in the following work.
We then tested the RNA expression level of HSPA7 in GBM cell lines (U87MG, U251MG, A172, and LN229 cells), and discovered that its expression was significantly higher in GBM cell lines than in normal human astrocyte (NHA) cells ( Figure 2C). The m 6 A methylation peak distribution and abundance in HSPA7 transcripts from GBM and normal brain tissues, as detected by m 6 A-seq, were visualized using IGV software ( Figure 2D). We found that HSPA7 was highly enriched in the m 6 A-precipitated fraction, and the m 6 A modification enrichment level could be regulated by methyltransferase-like 3 (METTL3), which contains a catalytic activity domain to catalyze m 6 A formation ( Figure 2E). Additionally, the expression of HSPA7 was significantly inhibited by knocking down METTL3 ( Figure 2F) and overexpressing FTO, an m 6 A demethylase ( Figure 2G). These results provided evidence that HSPA7 harbored high m 6 A modification levels and thus could be regulated in an m 6 A-dependent manner.
To explore the association between HSPA7 expression and clinical characteristics, we first compared HSPA7 expression levels in patients in the TCGA GBM cohort stratified separately by molecular subtype, IDH1 status, CpG island methylator phenotype (G-CIMP) status, MGMT promoter status, age, and sex. As shown in Figure 2H, HSPA7 expression in the PN subtype was significantly lower than the HSPA7 expression in the classical (CL) and MES subtypes and was highest in mesenchymal samples. HSPA7 expression was lower in samples with IDH mutations than in samples with wild-type IDH1. Regarding the G-CIMP status, HSPA7 expression was lower in patients with G-CIMP tumors than in those without G-CIMP tumors. No obvious correlations between HSPA7 expression and MGMT promoter methylation, age, or sex were observed. Furthermore, univariate Cox regression analysis of the overall survival of GBM patients in the TCGA cohort showed that high HSPA7 expression (HR: 1.511, P = 0.056) was an independent risk factor associated with the prognosis of GBM. Moreover, high HSPA7 expression (HR: 1.418, P = 0.034) remained a statistically significant factor in GBM patients after adjustment for age, sex, IDH status, MGMT promoter methylation status, and G-CIMP status in subsequent multivariate Cox regression analysis ( Figure 2I). These results indicated that the pseudogene HSPA7 is a novel risk prognostic biomarker and indicates therapeutic outcomes of GBM patients.

HSPA7 Is Correlated With Immunophenotypes and TME Landscapes
To gain further insight into the exact role of HSPA7 in immunophenotype determination, we analyzed 31 immuneassociated gene sets representing diverse immune cell types, functions, and pathways (see Materials and Methods). As shown in Figures 3A, B, the HSPA7-high expression group had significantly greater cell infiltration into the TME, higher immune and stromal scores, and lower tumor purity than the HSPA7-low expression group, confirming that HSPA7 could indeed regulate immune cell infiltration and immune-related gene expression. We then explored the specific differences in 31 immune cell phenotypes with high and low HSPA7 expression. Compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased infiltration of immunosuppressive cell populations, such as macrophages, neutrophils, and Tregs. However, some immuneactivating cells, including activated dendritic cells (aDCs), immature DCs (iDCs), plasmacytoid DCs (pDCs), and tumor-infiltrating lymphocytes (TILs), were also enriched (Supplementary Figure  S3A), indicating the complexity of the TME, in which GBM cells elicit multiple biological behavioral changes through direct or indirect interactions with other TME components. Two recent publications (7, 8) on single-cell mapping of human brain cancer reported that tumorassociated macrophages (TAMs) are the most abundant cellular components in brain TMEs, which can be subdivided ontogenetically into tissue-resident microglia (MGs) and macrophages of embryonic origin, and bone marrow-derived macrophages (BMDMs), showing tumor-promoting and immunosuppressive functions. To identify specific macrophage cell populations linked to HSPA7 expression in GBM, we examined the TCGA GBM dataset for TAM MG and TAM BMDM using validated gene set signatures (25), which demonstrated that tumors with high HSPA7 expression exhibited significantly increased infiltration of TAM BMDMs, while TAM MG was significantly decreased ( Figure 3A and Supplementary Figure S3A), suggesting that HSPA7 enhanced recruitment of tumor-promoting BMDMs into the GBM TME rather than MG. These studies also showed that both activation and exhaustion of lymphocytes were prevalent, with increased relative frequencies of Tregs in the brain tumor TME.
Moreover, we found elevated neutrophil infiltration in tumor tissues, revealing the complexity and multifaceted functions of the brain TME. We found that HSPA7 had a significant positive correlation with the enrichment scores of the three immunosuppressive immune cells and immune checkpoints ( Figure 3C). We next investigated chemokines and immune modulators associated with immune suppression states. We then explored the specific differences in markers of myeloid lineages with suppressive functions (CD33, NOS2, CD163, CD68, SPP1, CD14, CD206), immune inhibitory checkpoints (LAG3, CTLA4, HAVCR2, PDCD1, PDCD1LG2, CD27, TIGIT, CD274, and TNFRSF9), and major neutrophil-recruiting chemokines and their receptors (ITGA3, CD177, MET, and CXCL8) between patients with high and low HSPA7 expression. Tumors with high HSPA7 expression exhibited significantly increased expression levels of these markers compared to tumors with low HSPA7 expression ( Figures 3D-F), and positive correlations were found between HSPA7 expression and these molecules (Supplementary Figures  S4A-C). Furthermore, to explore the direct involvement of the HSPA7 in the biological pathway causing immune suppression, we then analyzed myeloid cell-derived macrophage-restricted chemokines, which were the main factors that cause immunosuppression in GBM. We found that compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased myeloid cell-derived macrophagerestricted chemokines (8,25) (Supplementary Figures S4D, E), including CCL17, CXCL2, CXCL3, and CXCL16 (involved in wound healing); immunosuppressive cytokines IL-10, CTSB, and CTSW (participating in multiple tumor-promoting processes, including invasion and metastasis); CCL2 and CCR2 (involved in macrophage chemotaxis); and other chemokine contributions of myeloid cell populations to the inflammatory TME milieu, indicating that HSPA7 could promote immunosuppressive phenotypes and suppress intratumoral antitumor immune responses.
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 functional gene sets identified by Bindea et al. (24). This analysis identified a TME pattern with two clusters corresponding to a TME immune-low (immune-L) and a TME immune-high (immune-H) phenotype (Supplementary Figure S3B). Patients in the immune-H group (n = 51) exhibited poorer prognosis ( Figure 3G) and higher HSPA7 expression levels ( Figure 3H) than patients in the immune-L group (n = 102). The PN, CL, and MES subtypes of GBM have been most consistently described in the literature; the PN subtype is related to a more favorable outcome, and the MES subtype is related to poorer survival (23). The association between the MES gene expression signature, characterized by NF1 mutation, with reduced tumor purity, elevated invasion, enhanced migration capacity, and infiltration of immunosuppressive cells (macrophages, microglia, mesenchymal stem cells, or other cells) has been identified as a common theme across cancers (33,34). Our results showed that most MES GBM samples had high expression of HSPA7 as well as the immune-H phenotype (Figures 3I, J). Moreover, we observed the same results as all GBM samples by analyzing MES glioblastoma alone, in which tumors with high HSPA7 expression exhibited  Figures S5A, B). As MES glioblastomas are strongly associated with higher immunosuppressive cell infiltration, we hypothesized that HSPA7 activated the immune microenvironment by promoting the phenotypic transformation of the GBM MES subtype. As shown in Supplementary Figure S5C, compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased expression levels of genes in the MES phenotype signature genes (33). Furthermore, HSPA7 positively correlated with CD44, a marker of the MES subtype ( Figure 3K). GSEA also showed that GBM samples with high HSPA7 expression were enriched in the MES subtype compared to GBM samples with low HSPA7 expression ( Figure 3L). In addition, we found that HSPA7 expression was higher in MES GSCs than in PN GSCs ( Figure 3M). Moreover, Western blot assays revealed that knockdown of HSPA7 reduced the expression of CD44 ( Figures 3N, O). Overall, these results indicated that in GBM, HSPA7 may be a robust indicator of the immunophenotype and be significantly correlated with a poorer immune response.

HSPA7 Is Correlated With Stromal and Carcinogenic Activation Pathways
To explore the differences in biological behaviors among these samples with distinct HSPA7 expression levels, we used the GSVA algorithm to estimate pathway enrichment scores for each sample (see Materials and Methods). Compared to the low HSPA7 expression group, the high HSPA7 expression group exhibited marked enrichment of stromal activation pathways [angiogenesis, epithelial-mesenchymal transition (EMT), VEGF, and TGF-b signaling pathways], oncogenic signaling pathways (hypoxia, apoptosis, PI3K Akt MTOR signaling, glycolysis, KRAS signaling, and other pathways), and immune responses (IL6 Jak stat3 signaling, which exerts immunosuppressive effects on T cell function and mediates ICB resistance in cancers (35), inflammatory response; interferon response; complement; and other pathways). However, the high HSPA7 expression group exhibited lower enrichment in pathways related to DNA replication-and DNA damage responserelated functions ( Figure 4A). Previous studies demonstrated that activation of stromal cells in the TME can induce T cell suppression and oncogenic signaling pathway activation via complex cellular and biological reconfiguration mechanisms (36,37), attenuating the tumor response to PD-L1 blockade (28). As shown in Figure 4B, stromal activation pathways were positively correlated with carcinogenic signaling pathways and immune suppression pathways but negatively correlated with DNA damage response and repair pathways in TCGA GBM samples. Further analyses of the activity of stroma-related pathways indicated that high HSPA7 expression was significantly associated with higher stromal activation signatures ( Figure 4C), as constructed by Mariathasan et al. (28), positively correlated with fibroblast enrichment (Figure 4D), as calculated by the MCP-counter method (38), and positively correlated with typical stromal cell activation-related pathways ( Figure 4E). Furthermore, TGF-b is a pleiotropic cytokine associated with poor prognosis in GBM, playing a protumorigenic role by promoting immunosuppression, angiogenesis, metastasis, mesenchymal transition, and fibroblast activation (28,(39)(40)(41), and has been proven to promote the progression of GBM via an autocrine signaling loop (42). In our analysis, a TGF-b ligand (TGFB1) and a TGF-b receptor (TGFBR2), two key regulators in stromal activation and EMT pathways (28), and other genes encoding ECM and matricellular proteins (COLA4A1, COLA4A2, ECM1, and FN1) (8), which form a barrier to lymphocyte infiltration, showed increased expression in the high-HSPA7 group compared to the low-HSPA7 group ( Figure 4F). These results suggested that immune cells and stromal cells in the TME can cooperate to synergistically regulate the immunosuppressive microenvironment, thereby promoting tumor immune escape. Then, we explored the pathways enriched with genes positively correlated with HSPA7 (Pearson's r ≥ 0.3, P ≤ 0.05, Table S7) via the Metascape database (see Materials and Methods). These genes were significantly enriched in pathways related to the immune response, including cytokine-mediated signaling pathways, leukocyte migration, differentiation, and other immune-related pathways. Additionally, they were enriched in pathways involving stromal activation such as extracellular structure organization and response to wounding. These enriched pathways interacted with each other to form a protein-protein interaction network ( Figure 5A). Furthermore, enrichment analysis in the PaGenBase (43) showed that these genes were almost completely specifically expressed in spleen, blood, bone marrow, and some other tissues with peripheral immune cell aggregation ( Figure 5B), suggesting that HSPA7 can indeed promote the infiltration of immune cells into tumor tissues. In addition, enrichment analysis in the TF-target interaction database Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) (44) showed that the genes were substantially transcriptionally regulated by NF-kB1 and RELA, two important transcription factors involved in NF-kB signaling pathways ( Figure 5C). Further GSEA using the Gene Ontology (GO) and KEGG databases also showed that HSPA7 expression was positively correlated with the immune response and extracellular structure organization ( Figure 5D).

Verification of the Functions of HSPA7 in Two CGGA Cohorts
To further validate the function of HSPA7 in GBM, we explored its expression pattern in two CGGA RNA-seq cohorts. The expression of HSPA7 was highest in GBM (WHO IV) among the glioma specimens with three different WHO grades (Supplementary Figures S6A, S7A). Kaplan-Meier survival analysis showed that GBM patients in both cohorts with high HSPA7 expression had poor survival outcomes ( Supplementary  Figures S6B, S7B). Further GSVA enrichment analysis also showed that HSPA7 expression was correlated with the immunophenotype, stromal activation, and oncogenic pathway activation in GBM samples (Supplementary Figures S6C, S7C). These results confirmed that HSPA7 can be used as an independent prognostic risk biomarker in GBM patients. Moreover, HSPA7 can regulate the TME immune response and stromal activation, promoting malignant progression of GBM tumors.

M 6 A-Modified HSPA7 Is Potentially Regulated by the Methyltransferase WTAP
To explore the potential molecular mechanism regulating HSPA7 RNA metabolism, we first screened HSPA7-interacting proteins via the NCBI database (Table S8). We verified that WTAP (45), a regulatory subunit of the m 6 A methyltransferase (46,47), is required for the localization of the m6A methyltransferase into nuclear speckles as a putative HSPA7binding regulator. To better map the association between WTAP and HSPA7, we explored the expression pattern in different grades of glioma in the TCGA, CGGA, Gravendeel and Rembrandt datasets. As shown in Supplementary Figure S8A, the expression of WTAP was highest in the WHO IV (GBM) group among glioma samples of three WHO grades and normal brain samples. Furthermore, overexpression of WTAP correlated with poor overall survival of GBM patients in all four datasets (Supplementary Figure S8B). Moreover, the expression of WTAP was positively correlated with the expression of HSPA7 in the TCGA GBM dataset (as calculated via the GEPIA database) and two CGGA GBM datasets (Supplementary Figure S8C). Further GSVA enrichment analysis showed that compared to the WTAP-low group, the WTAP-high group showed increased immune cell infiltration into the TME, stromal activation and oncogenic pathway activation, accompanied by higher immune and stromal scores and lower tumor purity, in the TCGA GBM cohort (Supplementary Figure  S8D), CGGA GBM cohort 1 (Supplementary Figure S9A) and CGGA GBM cohort 2 (Supplementary Figure S9B). These patterns were the same as the patterns observed for HSPA7 expression. The above findings indicate that m 6 A methylation of HSPA7 is regulated by WTAP, although this conclusion and the functional mechanisms need further demonstration. Subsequent pathway and process enrichment analyses of HSPA7-interacting proteins showed that they can mediate many biological behaviors, such as mitotic cell cycle phase transition, neddylation, proteasome-mediated ubiquitin-dependent protein catabolism, HIF1a pathway activity, and other functions (Supplementary Figure S10A, Table S8). Then, we further explored the enriched terms across the RNA-binding proteins (RBPs) detected by crosslinking immunoprecipitation coupled with sequencing (CLIP-seq) from data deposited in the starBase database. These RBPs were enriched mainly in terms such as mRNA processing, nucleic acid transport, and RNA catabolic process (Supplementary Figure S10B and Table S9). However, the exact mechanism by which HSPA7 regulates the GBM TME requires further clarification.

HSPA7 Holds Promise for Predicting the Therapeutic Response to ICB Therapy
We desired to further investigate the capacity of HSPA7 to predict the response to immune checkpoint therapy in GBM. However, few published immunotherapy datasets include GBM patients who received ICB therapy. Thus, we used a cohort of melanoma patients who received anti-CTLA4 therapy (23) and urothelial cancer cohorts of patients who received anti-PD-1 therapy (IMvigor210) to perform a complementary evaluation of the ability of HSPA7 to predict the immunotherapy response. Patients in the anti-CTLA4 cohort with low HSPA7 expression exhibited significant clinical benefits [ Figure 6A, anti-CTLA cohort, HR 1.927 (0.8997-4.129)]. Furthermore, compared to the HSPA7-high group, the HSPA7-low group exhibited a significant therapeutic benefit and clinical response to antiCTLA-4 immunotherapy (Figures 6B, C). Similar results were found in the anti-PD-L1 cohort. Kaplan-Meier analysis showed that the patients with the 50 lowest HSPA7 expression levels exhibited markedly prolonged survival compared to the patients with the 50 highest HSPA7 expression levels [ Figure 6D, HR 1.670 (1.022-2.792)]. Additionally, the HSPA7-low group exhibited significantly better therapeutic and clinical responses to anti-PD-L1 immunotherapy than the HSPA7-high group ( Figures 6E, F).
Further GSVA enrichment analysis showed that compared to the HSPA7-low group, the HSPA7-high group exhibited higher TME immune cell infiltration with higher immune and stromal scores, lower tumor purity (Supplementary Figure S11A), higher stromal activation, higher oncogenic pathway activation, and lower MMR pathway activation (Supplementary Figure  S11B), completely consistent with the functional enrichment patterns in GBM ( Figures 3A, 4A and Supplementary Figures  S5C, S6C). PD-L1 expression in immune cells (ICs) and tumor cells (TCs) was also assessed in the IMvigor210 cohort, and we examined the difference in HSPA7 expression among groups with different PD-L1 expression levels. As indicated in Figures 6G, H, patients with higher PD-L1 expression levels in either immune cells or tumor cells exhibited higher HSPA7 expression levels (ANOVA summary IC: P = 0.0001, TC: P = 0.0009), indicating that HSPA7 can upregulate PD-L1 expression, suppressing immune activation. Moreover, previous studies indicated that F-box and WD repeat domain containing 7 (FBXW7) is a vital tumor suppressor in various cancers, controlling proteasome-mediated degradation of oncoproteins such as cyclin E, c-Myc, Mcl-1, mTOR, Jun, and Notch (48) and that its loss-of-function mutation promotes resistance to anti-PD-1 therapy through downregulation of vital sensing pathways (49). Compared to the FBXW7 wildtype group, HSPA7 expression was obviously decreased in the FBXW7 mutant group ( Figure 6I). Fibroblast growth factor receptor 1 (FGFR1) is frequently mutated in various tumors, and inhibitors of FGFR1 have shown promising therapeutic value in several preclinical models (50). Palakurthi S. et al. (51) found that the combination of FGFR inhibition and PD-1 blockade can promote tumor-intrinsic induction of antitumor immunity. We also found that HSPA7 expression was significantly higher in the FGFR mutant group than in the FGFR wild-type group ( Figure 6J). Collectively, our work strongly indicates that HSPA7 expression contributes to predicting the response to immune checkpoint therapy.

HSPA7 Facilitated Macrophage Infiltration and Could Be a Potential Immunotherapy Target for GBM Patients
We then analyzed the differences in the distribution of somatic mutations between the low and high HSPA7 expression groups in the TCGA GBM cohort using the "maftools" package and found that the PTEN mutation rate was significantly increased in the HSPA7-high group compared to the HSPA7-low group ( Figure 6K; low: 20%, high: 36%). Additionally, the NF1 mutation rate, which often occurs in the MES subtype and drives recruitment and activation of TAMs (33), was also obviously elevated in the HSPA7-high group compared to the HSPA7-low group ( Figure 6K; low: 8%, high: 11%). Chen et al. (52) found that PTEN deficiency in GBM increases macrophage infiltration by activating YAP1 signaling, which directly upregulates lysyl oxidase (LOX) expression, an MES subtype marker and a macrophage chemoattractant (33), in turn inducing SPP1 secretion to support GBM survival. We found that YAP1 signaling was significantly upregulated in the HSPA7high group compared to the HSPA7-low group ( Figure 7A). We also confirmed that LOX and YAP1 were higher in MES GSCs than in PN GSCs ( Figure 7B). Moreover, LOX and YAP1 expression levels were enhanced upon overexpression of HSPA7 in PN subtype GSCs 8-11 ( Figures 7C, D) but decreased significantly upon knockdown of HSPA7 in MES subtype GSCs 267 ( Figure 7D), suggesting that HSPA7 may facilitate tumor promoting macrophage infiltration by enhancing LOX expression, which could promote macrophage migration into the GBM TME and enhance angiogenesis by upregulating TAMderived SPP1 (52). Finally, to further confirm the role of HSPA7 in facilitating macrophage migration by enhancing the YAP1-LOX axis, using the transwell assay, we found that conditioned medium from HSPA7 knockdown GSCs 267 inhibited THP1-differentiated macrophage migration significantly compared to the NC group ( Figure 7E). Then, THP1-differentiated macrophages were cultured in conditioned medium collected from GSCs 267 transfected with HSPA7 ASO or NC. As shown in Figure 7F, conditioned medium from HSPA7-knockdown GSCs inhibited SPP1 expression in macrophages compared to the NC group. The immunofluorescent staining of our local clinical tumor tissue section analyses also showed that the expression of YAP1, LOX, CD68 (human macrophage marker), SPP1, and PD-L1 was enhanced in HSPA7-high GBM tissues compared to HSPA7-low tissues ( Figure 7G). Moreover, Robert M. Samstein et al. (53) found that TMB predicted survival after immunotherapy across multiple cancer types (including GBM). Thus, we downloaded the mutation data of 82 GBM samples and then analyzed the differences in the distribution of somatic mutations between the samples with the 30 lowest TMBs and the 30 highest TMBs. The PTEN and NF1 mutation rates were significantly increased in low-TMB samples compared to the high-TMB samples, and the same results were found for the HSPA7 expression level (Supplementary Figure  S12A). We further explored the association between the HSPA7 expression level and TMB. As shown in Supplementary Figure  S12B, HSPA7 expression was statistically negatively correlated with TMB. In addition, Touat M. et al. (54) recently found that mismatch repair (MMR)-deficient gliomas were characterized by a lack of prominent T cell infiltration, extensive intratumoral heterogeneity, poor patient survival, and a low rate of response to PD-1 blockade therapy. Our data also showed that the high HSPA7 expression group showed reduced activation of MMRassociated pathways ( Figure 4A), which was negatively correlated with immunosuppression, stromal activation, and oncogene pathway activation ( Figure 4B). Other studies have also shown the same results (28,55). Furthermore, HSPA7 expression was negatively correlated with the expression of MMR genes such as MLH1, MSH2, MSH6, and PMS2 (Supplementary Figure S12C). Previous studies indicated that GBM patients exhibit resistance to immune checkpoint inhibitors (ICIs) due to the low mutation rate, the PTENdeficient immunosuppressive microenvironment, infiltration by myeloid-derived suppressor cells (including TAMs, prominent players in brain cancers), and activation of tumor stromal cells (56,57), indicating that knockdown of HSPA7 may enhance the efficacy of ICIs such as PD1 inhibitors in GBM. To illustrate this function, we generated patient-derived glioblastoma organoids (GBOs), which maintain cell-type heterogeneity (including macrophages, T cells, and vascular cells) and molecular signatures of their respective parental tumors, according to experimental protocols already reported by Jacob F. et al. (31). Immunofluorescence assays of our GBO sections also confirmed the presence of macrophages and vascular cells, as detected by CD68 and CD31 markers, respectively (Supplementary Figure  S13A), supporting the use of three-dimensional GBO models for holistic study of the TME of GBM and evaluating the efficacy of PD1 inhibitors. We next applied our GBO model to test treatment responses in vitro. To mimic the postsurgical standard of care treatment, we subjected GBO samples from patient whose GBM tissues expressed high CD68 and SPP1 (Supplementary Figure  S13B) to a single exposure conditioned medium from HSPA7knockdown GSCs or with PD1 inhibitor (5 µM) treatment for 5 days. The therapeutic response was evaluated by quantifying the percentage of cells expressing Ki67 (proliferation index) and CD44 (invasion index). As shown in Figure 7H, knockdown of HSPA7 significantly enhanced the efficacy of the PD1 inhibitor in GBM. In summary, HSPA7 facilitated tumor promoting macrophage migration into the GBM TME by activating the YAP1-LOX axis, and our results indicated that HSPA7 might be a potential immunotherapy target for GBM patients.

Characterization of HSPA7 Across 33 Cancer Types
We investigated whether HSPA7 expression is correlated with prognosis in a pancancer patient cohort. The expression of HSPA7 in tumor tissues and GTEx normal brain tissues was determined in the GEPIA database. As shown in Supplementary Figure S14, the expression of HSPA7 was significantly lower in tumor tissues in the ACC, COAD, DLBC, LAML, LUAD, LUSC, READ, THCA, and THYM cohorts but significantly higher in tumor tissues in the GBM, KIRC, KIRP, and PAAD cohorts than in the corresponding normal tissues. Genetic alterations were also observed in many other tumors, although these alterations were not significant. In addition, Cox regression analysis of survival rates and Kaplan-Meier analysis were conducted using the SangerBox tool. The relationships between HSPA7 expression levels and the prognoses of different tumors are shown in Supplementary Figures S15-S17. Statistically, HSPA7 is a risk factor in most cancers (such as ACC, LDBC, and CESC). However, we found that it is a relatively favorable factor in several cancer types such as SKCM, KIRP, and SARC. We next evaluated the immune score and stromal score across 33 cancers in the TCGA database. HSPA7 expression was significantly positively correlated with the immune score (Supplementary Figure S18 , and quantification histogram represented the relative protein expression of LOX and YAP1 (right); data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001. Means were compared with Student's t-test for two groups and one-way ANOVA for three groups, and the NC group is indicated as the control. (E) Representative transwell migration assays showed that HSPA7 inhibited human THP1-differentiated macrophage migration by exposing them to conditioned medium from GSC 267 cells transfected with NC, HSPA7 ASO1, or ASO2 for 24 h. Original magnification, ×100 (left), and quantification histogram represented relative migration of THP-1 macrophages (right); data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001. Means were compared with one-way ANOVA, and the NC group is indicated as the control. (F) Western blot assays showed that HSPA7 inhibited the expression of SPP1 in THP1-differentiated macrophages by exposing them to conditioned medium from GSC 267 cells transfected with NC, HSPA7 ASO1, or ASO2 for 48 h. (G) IF staining in a human GBM tissue microarray showed that the expression of LOX, YAP1, PD-L1, CD68, and SPP1 was higher in the HSPA7 high group than in the low group. Histogram representing statistical proportion data of positive area; data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001, means were compared with Student's t-test. (H) GBOs at 2 weeks were cocultured with conditioned medium from GSC 267 cells, transfected with NC, HSPA7 ASO1, or ASO2 at 48 h, and PD1 antibody (5 mM) for 5 days as indicated. IF staining for Ki67 and CD44 in GBO sections showed that knockdown of HSPA7 enhanced the effect of anti-PD1 therapy, original magnification, ×630 (scale bars: 20 mm). Histogram representing statistical proportion data of positive area; data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001. Means were compared with one-way ANOVA, and the NC group is indicated as the control.
types such as GBM, LGG, OV, and LUAD. However, no relationship was found in other cancers such as SKCM, UVM, TGCT, and PAAD, indicating that HSPA7 may act via a mechanism other than TME regulation in these cancers. Thus, the molecular mechanisms of this gene need further study in specific tumors. Collectively, these data emphasize that HSPA7 may be a key factor in facilitating the acquisition of various immunophenotypes in various cancers and may be considered in the development of more effective immunotherapies for GBM and other immunotherapy-resistant tumors. Because different immune cells may play different roles in different tumors, the functions of HSPA7 in specific tumors must be explored.

DISCUSSION
GBM, one of the most aggressive brain tumors, currently has no effective and sufficient therapies due to its intratumoral heterogeneity and molecular plasticity. ICB therapy is being actively pursued as a promising treatment option for GBM, but very few patients respond to ICB therapy. Thus, identifying markers regulating the brain TME that are prominent players in ICB therapy for cancer could reveal promising new targets for therapeutic intervention. LncRNAs and m 6 A modifications are emerging as indispensable regulators of the TME. However, the overall TME infiltration characteristics mediated by m 6 Amodified lncRNAs have not been comprehensively recognized. Therefore, it is worth obtaining comprehensive knowledge of the cellular TME infiltration characteristics mediated by m 6 Aregulated lncRNAs.
Here, based on our m 6 A-seq data, we revealed highly distinct lncRNA m 6 A methylation modification patterns between GBM and normal brain tissues. In addition, we identified immune-stromalm 6 A-related HSPA7 as a novel prognostic risk factor in GBM patients; this gene plays a crucial role in immunophenotype determination, stromal activation, and oncogenic pathway activation and has a robust capacity to predict the ICB response. Furthermore, in this study, WTAP, a methyltransferase mediating m 6 A modification, was identified as a potential regulator of HSPA7 expression. Further analysis showed that WTAP, like HSPA7, can also regulate immunophenotype determination and stromal activation-related pathways (Supplementary Figures S8, S9). This finding again demonstrated that m 6 A modification is highly important in shaping the TME landscape. Many studies have reported that stromal cell activation in the TME can suppress immune infiltration or facilitate an immunosuppressive response in the TME, mediating therapeutic resistance to ICB (28,58). Our data also revealed a markedly negative correlation between HSPA7 expression and TMB, a marker of the response to ICB therapy, in a TCGA GBM cohort ( Figure 6B). In addition, another study showed that PTEN mutations were associated with immunosuppressive expression signatures in non-responders to anti-PD-1 immunotherapy in GBM (2). Chen et al. (52) found that PTEN deficiency in GBM increases the infiltration of SPP1 + macrophages, which can interact with fibroblasts and vascular endothelial cells, inducing angiogenesis, EMT, and some other stromal activation-related pathways in colon cancer and are characterized by expression of the pattern recognition receptor MARCO (59) via a YAP1-LOX-b1 integrin-PYK2 axis. Intriguingly, antibodymediated depletion of MARCO inhibits cancer progression and metastasis, enhancing ICB efficacy (60), suggesting that SPP1 + macrophages play a non-negligible role in the immunorepressive response and immunotherapeutic resistance. We then analyzed myeloid cell-derived macrophage-restricted chemokines, which were the main factors that cause immunosuppression in GBM. We found that compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased myeloid cell-derived macrophage-restricted chemokines (8,25) (Supplementary Figures S4D, E); we also confirmed that HSPA7 could facilitate the macrophage infiltration into the GBM TME via YAP1-LOX axis in vitro, as well as into our local GBM tissue sections (Figures 7A-G). Furthermore, we found that HSPA7 can interact with YAP1 (61), which can induce the secretion of CCL2/CSF1 to recruit monocytes (62), suggesting that HSPA7 may be a target that synergistically regulates SPP1 + macrophages, which then induce stromal activation in the TME. Moreover, we found that HSPA7 and YAP are also contained in GBM extracellular exosomes (63), a fundamental regulator of TME cells. Numerous research reports have indicated that YAP is a hub in the network of signals exchanged within the TME (64), thus regulating the immune response and stromal activation. However, the specific mechanism of HSPA7 needs to be proven by a large number of experiments. We also confirmed the predictive value of HSPA7 in a cohort of melanoma patients who received anti-CTLA4 therapy and an IMvigor210 cohort of patients treated with anti-PD-L1 therapy. A statistically significant difference in HSPA7 expression was found between non-responders and responders. MMR deficiency has recently emerged as a beneficial indicator of the response to PD-1 blockade in patients with cancer (65, 66), whereas we found that HSPA7 was negatively correlated with MMR-related pathway activity ( Figure 4A), which in turn was negatively correlated with stromal activation and oncogenic pathway activation ( Figure 4B and Supplementary Figure S11). This pattern suggests that MMR deficiency is an unfavorable factor for the response to ICB therapy. Other studies have shown the same results (28,55). By comprehensive analysis of data for 10,294 samples, a recent study showed that MMR-deficient gliomas are characterized by a lack of prominent T cell infiltration, extensive intratumoral heterogeneity, poor patient survival, and a low rate of response to PD-1 blockade (54), reconfirming the ability of HSPA7 to predict the response to immunotherapy. We also explored the function of HSPA7 in various cancer types. As shown in Figures S13-S19, HSPA7 expression had prognostic significance in most cancers; in addition, it can regulate the expression of immune checkpoint genes and the activity of immune response pathways and was positively correlated with the immune score and stromal score in the majority of tumor types. HSPA7 was found to be a risk factor in the cohort of melanoma patients who received anti-CTLA4 therapy ( Figures 6A-C) but was found to be a beneficial prognostic factor in the TCGA SKCM cohort (23) (Supplementary Figures  S14, S16), exhibiting no obvious correlation with immune checkpoint expression (Supplementary Figure S20A) or the immune response (Supplementary Figure S20B). The mechanisms underlying this discrepancy need future exploration.
In summary, via integrated analysis of our own m 6 A-seq data and public clinical data, we found highly distinct lncRNA m 6 A methylation modification patterns between GBM and normal brain tissues and determined that HSPA7 could be used to evaluate the prognosis of GBM patients, as well as the corresponding cellular TME infiltration and activation characteristics of individual patients, and could predict the clinical efficacy of anti-PD-1/-PD-L1 immunotherapy in patients. More importantly, this study offers several novel insights regarding cancer immunotherapy; for example, targeting m 6 A regulators to change the m 6 A modification patterns of key immune-regulating targets may reverse unfavorable cellular TME infiltration characteristics. This study contributes to future exploration of novel drug combination strategies or novel immunotherapeutic agents.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, SRA PRJNA661159.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethical Committee on Scientific Research of Shandong University Qilu Hospital (approval number: KYLL-2018-324). The patients/participants provided their written informed consent to participate in this study.