香京julia种子在线播放

RNA-seq data and clinical data on LGGs were extracted from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/). After transcripts per million (TPM) conversion, we held genes with expression levels larger than 0.1 TPM for analysis. Matched Normal data: 105 cortex tissues were obtained from the GTEx project (https://commonfund.nih.gov/GTEx/). We downloaded both methylation data (Illumina Infinium HumanMethylation450 BeadChip) and somatic mutation data generated by TCGA from the UCSC Xena browser (https://xenabrowser.net/hub/) (16). Fibroblast and glioma cell lines from The Cancer Cell Line Encyclopedia (CCLE) (https://portals.broadinstitute.org/ccle) (17). Infiltration Estimation for TCGA-LGGs was collected from TIMER2.0 (http://timer.comp-genomics.org/). The STRING (https://www.string-db.org) database produced the protein-protein interaction network and reconstructed them via Cytoscape software. A flowchart of the entire procedure can be shown in Figure 1 .

The Illumina HumanMethylation450 BeadChip array contains probes covering 99% of reference sequence genes and 96% of CpG islands. The raw methylation intensities for each probe were represented as β-values, which were converted into M-values with R package Lumi for statistics analysis (18). 5’-C-phosphate-G-3’ (CpG) methylation data between different groups were compared with R package limma to identify differentially methylated CpG sites. Benjamini-Hochberg (BH) method was used to adjust p-value as a false discovery rate (FDR). The CpG site and gene mapping files were downloaded from illumine (https://www.illumina.com/). CpGs span various gene regions, including 1500 bp and 200 bp upstream of the transcription start sites (TSS1500 and TSS200, respectively). The average β-values for each region were calculated according to all CpG sites at the corresponding region, and the average β-value was converted to an M-value. Average regional methylation data between different groups were compared with R package limma to identify differentially methylated regions (DMRs). The hypermethylated DMRs with a threshold of adjusted p-value< 0.05 combined a delta β-value > 0.2, and the hypomethylated DMRs with a threshold of adjusted p-value< 0.05 combined a delta β-value< −0.2. Genes harboring DMRs in any part of the gene features were differentially methylated genes (DMGs). The equations described above are listed below.

where M is the intensity of the methylated allele, i = 1, 2, 3,…, k, and k is the number of CpG sites in a region.

CAFs abundances were separately estimated via Estimate the Proportion of Immune and Cancer cells (EPIC) algorithm using the R package Immunedeconv (19). TCGA-LGGs were divided into a high-CAFs-score group and a low-CAFs-score group according to the median score. The estimated immune and stromal scores were computed using the R package ESTIMATE.

Differential expression between the high-CAFs-score and low-CAF-score group samples was analyzed with the R package limma. False discovery rate (FDR) as adjusted p-value using the Benjamini-Hochberg (BH) method. The fold change was log2-transformed. Differentially expressed genes (DEGs) were calculated with a difference > 1.5-fold and p< 0.01. Methylation analysis results were carried out for joint analysis. Venn diagram analyses were performed to calculate the intersection of DMGs and DEGs and explored the differentially methylated and expressed genes (DMEGs). DMEGs were grouped according to four expression patterns: HypoUp, HypoDown, HyperUp, and HyperDown.

To functionally annotate DMEGs of this study, Gene Ontology (GO) including biological process (BP), cellular component (CC), and molecular function (MF) analysis was performed in the R package ClusterProfiler (20). A p-value of< 0.05 and an FDR of< 0.05 were used for the cutoff value. The ClueGO Plugin version 2.5.8 in Cytoscape Version 3.8.2 was employed to identify hub genes and functional analysis (21).

We selected HypoDown and HyperUp genes in TSS200 and TSS1500 for univariate cox regression analysis and filtrated the prognostic-related genes. Subsequently, we used the R package glmnet to conduct the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm and develop a potential risk signature. The minimum value of lambda was derived from 1,000 cross-validations (‘1-se’ lambda), and which corresponding partial likelihood deviance value was the smallest for the risk model. Coefficients with regression were confirmed by the “cvfit” function with 1000 repeats. The risk score calculating equation, which contains five risk genes, is:

where Coef _i means the coefficients, x _i is the expression value of each gene.

The predictive power of the prognostic signature was evaluated by the receiver-operating characteristic (ROC) curve. The independent clinical factor was validated by multivariate Cox regression analysis. Finally, a nomogram was constructed according to independent predictors. The calibration of the nomogram was evaluated by the calibration curve to assess the goodness of 1-, 3-, and 5-year overall survival.

In this study, the mRNA expression matrix of LGGs was analyzed using the CIBERSORT R script downloaded from http://cibersort.stanford.edu. Based on deconvolution, we estimated the abundance of immune cell populations. The relationship between each immune cell and survival was measured by Kaplan-Meier (KM) survival analysis. We evaluated a total of 60 immune checkpoints (ICP) genes in two categories (Inhibitory ICP (22) and Stimulatory ICP (23)) from widely recognized literature (24). Then we assessed the expression and survival of these ICP in the TCGA-LGGs cohort for a comprehensive overview of the immunosuppressive landscape.

The Tumor Immunophenotype Profiling (TIP) was performed to quantify the extent of infiltrating immune cells and anticancer immunity (25). Assessment of antitumor immunity was conceptually divided into seven steps, including tumor cell antigen release (step 1), cancer antigen presentation (step 2), priming and activation (step 3), trafficking of immune cells to tumors (step 4), infiltration of immune cells (step 5), T cell recognition of cancer cells (step 6), and killing of cancer cells (step 7).

TIDE (http://tide.dfci.harvard.edu/), an excellent algorithm, was used to explore the prediction of clinical response to immune checkpoint blockade (ICB) therapy (22). The TIDE score was calculated to simulate two mechanisms of tumor immune evasion: the induction of T cell dysfunction with high infiltration of cytotoxic T lymphocytes and the retard of T cell infiltration in tumors with low cytotoxic T lymphocyte infiltration. The TIDE score is a good reflection of the responsiveness of the ICB. The SubMap (https://www.genepattern.org/) was carried out to validate the reliability of the prediction of TIDE.

The R package pRRophetic was used to predict chemotherapeutic response in LGGs patients (26). In addition to temozolomide, which patients with glioma widely used, we included three drugs with the therapeutic potential for glioma in this study: Axitinib, GDC-0941 (PI3K inhibitor), and Bleomycin.

H4, SW1783, and HMF cell lines were purchased from ATCC. According to the manufacturer’s instructions, total RNA was isolated using Trizol reagent (Invitrogen, USA). 2μg of the total RNA was transcribed into cDNA. SYBR Green PCR kit (Takara, Japan) was used for qRT-PCR. We selected the 2−ΔΔCq method to calculate gene transcription level, with β-actin mRNA as control. Data represent the mean ± SD of triplicate real-time PCR. The primers were synthesized by Tsingke Biotechnology (Shanghai, China) and displayed in Supplementary Table S1 . Immunohistochemistry (IHC) analyzed the protein expression levels. GSAP (ab106630), LATS2 (ab111054), SWAP70 (ab228846), and SLC2A10 (ab110528) antibody was purchased from Abcam. EMP3 (sc-81797) antibody was purchased from Santa Cruz Biotechnology. Clinical characteristics of LGGs patients cohort are displayed in Supplementary Table S2 . All the patients and the hospital’s Ethics Committee approved this research. LGGs tissues were formalin-fixed, paraffin-embedded, and sectioned at 4 µm. Immune complexes were detected with the SP Kit (Solarbio, Beijing, China) and DAB Substrate Kit (Solarbio, Beijing, China). Signals were detected using an Olympus BX41 microscope. Quantification of Immunohistochemistry (IHC) staining was performed in a blinded fashion.

All the data were analyzed using the R software version 4.1.0. The overall survival (OS) between different groups was analyzed using Kaplan-Meier curves. Kruskal–Wallis tests were applied to compare gene expression in two groups. The Fisher test assessed different groups’ responses to ICB treatment. Somatic mutation data sorted in the form of Mutation Annotation Format (MAF) was analyzed using the R package maftools. ImageJ and ImageJ plugin IHC profiler was applied to quantify IHC staining analysis. IHC scoring data and qRT-PCR data were analyzed using GraphPad/Prism 9.0. In addition, tumor mutational burden (TMB) and mutation counts were computed from somatic mutation frequencies. p< 0.05 was marked as ‘*’, p< 0.01 was marked as ‘**’, p< 0.001 was marked as ‘***’. and p< 0.0001 was marked as ‘****’.

To identify DMEGs in LGGs, we first extracted the gene expression and DNA methylation data of TCGA-LGGs and performed a comparative analysis. Samples were divided into high- and low-CAFs-score groups according to the median CAFs score. From the summary estimate, 2393 statistically significant Differentially Expressed Genes (DEGs) were identified (difference > 1.5-fold, p-value< 0.01), including 263 upregulated and 2131 downregulated genes ( Figure 2A , Table S3 ). Promoter regions (TSS200 and TSS1500) were enrolled in the primary study. As shown in volcano plots, 1276 DMGs were identified from two regions, including 896 DMGs in the TSS200 region ( Figure 2B , Table S4 ) and 583 DMGs in the TSS1500 region ( Figure 2C , Table S5 ).

We analyzed the relationship between methylation and gene expression by integrating DMGs and DEGs in two promoter regions (TSS1500 and TSS200). A total of 77 DMEGs were identified ( Figure 2D ), and then we performed a principal component analysis (PCA) of the DMEGs in normal tissue and LGGs ( Figure 2E ). The PCA profile revealed a clear separation of normal samples from LGGs. Afterwards, GO enrichment analysis was performed on these DMEGs, and the result indicated that DMEGs were significantly enriched in positive regulation of oligodendrocyte differentiation, negative regulation of canonical Wnt signaling pathway, and regulation of epithelial to mesenchymal transition ( Figure 2F ). These functions were closely related to the oncogenesis and progression of glioma.

To investigate the differences in DMEGs within each region, we classified these DMEGs into four groups (HypoUp, HyperUp, HyperDown, and HypoDown) for TSS200 and TSS1500, respectively ( Figures 3A, B ). The HyperDown group was the most prevalent in TSS200 and TSS1500 regions, and the HypoUp group had the second-highest proportion in the two regions ( Figure 3C ). After extracting the HyperDown and HypoUp DMEGs in TSS200 and TSS1500 regions separately, we used STRING to construct PPI networks. These genes were analyzed by GO enrichment using the “ClueGO” plugin for Cytoscape software (p< 0.01, Kappa score = 0.5). Functional enrichment analysis revealed that these DMEGs participated in critical biological processes. In the TSS200 region, DMEGs are mainly associated with response to mechanical stimulus, regulation of cell-substrate adhesion, and positive regulation of macrophage migration ( Figure 3D ). As shown in Figure 3D , MMP14 is critically involved in various immune-related functions. Meanwhile, DMEGs in the TSS1500 region are mainly related to regulating epithelial to mesenchymal transition, programmed necrotic cell death, and positive regulation of gliogenesis ( Figure 3E ). These simultaneous differential DEGs and DMGs pooled to DMEGs may be the main factor causing the altered biological function of LGGs.

Among the DMEGs in HyperDown and HypoUp groups, Univariate Cox regression analysis screened 203 DMEGs with prognostic values. Then a Lasso‐penalized Cox analysis was performed to shrink further the scope of DMEGs screening ( Figure 4A ) and lambda. Min was regarded as the optimal value in the cross-validation process ( Figure 4B ). Five DMEGs (GSAP, EMP3, LATS2, SWAP70, and SLC2A10) and corresponding coefficients ( Table S6 ) were identified. We used TCGA-LGGs as the train set CGGA-325 and CGGA-693 as the test sets, and samples were split into high- and low-risk groups by the median value of the risk score. KM survival curves depicted that LGG patients with increased risk scores had worse clinical outcomes in both train set and test sets ( Figures 4C, D , and Figure S1A ). Statistical analysis was performed using a log-rank test (train set p< 0.001, test set p< 0.001). After that, we established ROC curves of the risk score model with 1-year, 3-year, and 5-year. The results revealed that the risk score could effectively distinguish LGGs patients with different survival statuses ( Figure 4E 1-yer AUC = 0.86, 3-year AUC = 0.83, and 5-year AUC = 0.80). The results were similar and slightly lower in the test set ( Figure 4F 1-yer AUC = 0.73, 3-year AUC = 0.79, and 5-year AUC = 0.77). Multivariate Cox regression analysis showed the independent prognostic value of this risk score ( Figure 4G , p< 0.001, HR = 3.844). Then we examined the mRNA expression of five risk genes in the glioma cell lines and fibroblast cell lines with CCLE ( Figure 4H ). Although glioma lines are not representative of LGGs, we found the expression of GSAP was consistent with that of fibroblasts, and the expression of the other four risk genes was lower than that of fibroblasts in glioma cell lines. The mortality of patients increased with the increase in the risk score ( Figure 4I ). Finally, we established the Nomogram model, which contained risk score and age to assess the survival prediction in LGGs patients ( Figure 4J ) and calibration curves for nomogram predicted 1-, 3- and 5-year overall survival performed well with the risk model ( Figure 4K ).

Clinical variables were introduced into the risk score system to analyze the relationship between risk scores and clinical characteristics. We initially compared the CAFs’ scores concerning the risk score. As shown in Figure 5A , the risk score was mildly positively associated with the CAFs score (r = 0.42, p< 0.01). In contrast, tumor purity was negatively correlated with the risk score ( Figure 5B r = -0.56, p < 0.01). Other clinical features were then introduced. WHO grade III LGGs have a higher risk score than WHO grade II LGGs, which is consistent with their poor prognosis ( Figures 5C, D , Figure S1B ). The risk score in MGMT promoter methylated LGGs were significantly lower than that in MGMT promoter unmethylated LGGs samples ( Figures 5E, F , Figure S1C ). It is widely accepted that glioblastoma patients with MGMT promoter methylated are sensitive to temozolomide and suitable for TMZ chemotherapy. For another crucial molecular marker 1p19q codeletion status, the risk score was significantly higher in 1p19q-non-codeletion samples compared to codeletion samples in both datasets ( Figures 5G, H , Figure S1D ). Correspondingly, IDH1/2 wild-type cases showed a valid increased risk score compared to IDH1/2 mutant cases ( Figures 5I, J , Figure S1E ). The high-risk scores were seen in cases aged > 45 years. ( Figures 5K, L , Figure S1F ). Overall, some important molecular markers and clinical features of LGGs responded well to this risk score system.

To investigate the relationship between the prognostic signature and immune cell infiltration in LGGs. We evaluated immune scores, immune cell infiltration, and immune checkpoints separately. CAFs are the most prominent tumor stroma cell type in the TME. Comparing the stromal score and immune score of the LGGs datasets, we found a significantly positive correlation between the risk score and stromal score, and the same was true of the immune score ( Figures 6A, B , stromal score p< 0.001 and Figures 6C, D immune score p< 0.001). CIBERSORT algorithm showed the proportions of distinct immune cell subpopulations in different risk groups ( Figure 6E ). The relative expression is shown in bar diagrams. The proportions of Macrophages M0 and Macrophages M2 in the high-risk group were significantly higher than in the low-risk group ( Figure 6F ). By contrast, the proportions of Monocytes, activated NK-cell and activated Mast-cells were higher in the low-risk group. Subsequently, we performed KM survival analysis to evaluate the OS with differing immunocytes infiltration samples ( Figures 6G–I ). High proportions of activated NK-cell and activated Mast-cell had a better OS ( Figures 6G, I ). Conversely, the level of macrophage M0 expression is inversely related to OS ( Figure 6H ). We next examined whether the risk score is associated with the expression of inhibitory and stimulatory immune checkpoint (ICP) molecules. Fortunately, the expression of many immune checkpoint molecules showed significant differences between high and low-risk groups ( Figure 6J Inhibitory ICP and Figure 6K stimulatory ICP). The KM curves for survival analysis of these immune checkpoints were presented in Figure S2 . In summary, the risk score was associated with the expression of these immune checkpoint molecules.

To explore the risk score stratification and the associated characteristics of the antitumor immune response, we introduced the Tracking Tumor Immunophenotype (TIP) system, which analyzed the status of anticancer immunity and the proportion of tumor-infiltrating immune cells across seven-step Cancer-Immunity Cycle. As shown in Figure 7A , the scores of the high-risk group were significantly higher in tumor antigen release (step 1), immune cell recruitment (step 4), and immune cell infiltration into the tumor (step 5). Conversely, the low-risk group scored significantly higher in T cell priming and activation (step 3), T cell recognition of cancer cells (step 6), and killing of cancer cells (step 7). Then, we used the CGGA dataset to validate these results ( Figure S3A ). Next, we introduce the TIDE algorithm to assess the efficacy of risk score in predicting ICB responsiveness. We found there was a significant difference in response to ICB treatment between the two groups (p< 0.001), and the response to ICB treatment was more sensitive in the low-risk group ( Figures 7B, C ). SubMap was used to compare the prediction response to anti-PD1 and anti-CTLA4 therapy results with another dataset containing 47 patients with melanoma that responded to immunotherapies. Using this tool, we found that the high-risk group in the train and test sets showed comparable performance in predicting the LGGs’ response to anti-PD1 therapy ( Figures 7D, E p< 0.05). Anti-CTLA4 therapy also showed a partial response in the TCGA train set ( Figure 7D Bonferroni corrected p = 0.32).

We tried to analyze the response of two risk groups to chemo-drug efficacy. Here, in addition to TMZ, we selected three other chemo drugs from the literature that may have therapeutic potential for glioma. As expected, the sensitivity in the low-risk group was slightly better than that in the high-risk group ( Figure 7F ). That may be related to the level of MGMT promoter methylation in the low-risk group ( Figure 5E ). For the other three drugs, the estimated IC50 was significantly better in the low-risk group ( Figures 7G–I ). These results were validated using the CGGA325 dataset ( Figure S3B–E ). Although these drugs have not been used in clinical treatment on a large-scale, differences in sensitivity suggest that they have potential as novel therapeutic agents.

Tumor genomic alterations have profound effects on immunity and drug therapy. We investigated the mutation frequencies of different risk groups and showed the top ten most frequently mutated genes ( Figures 8A, B ). As we expected, IDH1 had the highest mutation frequency, predominantly missense. Remarkably, some genes were associated with the immune microenvironment and immunotherapy. These genes were more prominent in the high-risk group, such as ATRX, EGFR, and PTEN. We targeted IDH1, the most frequently mutated of the two groups in our study and analyzed the relationship between the expression of five risk genes and IDH1 mutations. The results were shown in Figures 8C–G . The expression of all risk genes was higher in samples with wild-type IDH1 than in mutant IDH1. These risk genes and corresponding risk scores are consistent with the findings summarized in the preceding text ( Figure 5I ). In addition, we analyzed correlations between tumor mutational burden (TMB) and risk score. Like previous research, LGGs patients in the high TMB group have a poorer prognosis. Hence, we introduced the risk score and analyzed it jointly with TMB. Our results show that the risk score had a low positive correlation with TMB values ( Figure 8H ). Meanwhile, the high-risk group with a poorer prognosis corresponded to high TMB values ( Figure 8I ). Combining the two elements in the analysis, we found that LGGs samples with high-risk scores and high TMB values had the worst prognosis ( Figure 8J ). Finally, we compared the TIDE, dysfunction, and exclusion scores between the different risk groups. The TIDE score in patients with low-risk scores was significantly higher than those with high-risk scores ( Figure 8K p< 0.01). In parallel, patients in the high-risk group had higher dysfunction scores constructed using dysregulated immune genes ( Figure 8L p< 0.001). No difference was found between the two groups for exclusion scores constructed using immune rejection genes ( Figure 8M ). The mechanism of immune escape in these high-risk LGGs samples may be immune dysfunction rather than immune exclusion. These epigenetic alterations may affect the prognostic model for the therapeutic assessment of LGGs patients.

After obtaining the above five risk genes, we identified them at the cellular and protein levels. Given the rarity of fibroblasts in brain tissue, we chose a stable human mammary fibroblast (HMF) cell line as a control group. T98G and U251 cell lines are commonly used glioma cell lines for experiments. GSAP and SWAP70 expressed similar levels in different cell lines, and the expression of the remaining three genes was lower in glioma cell lines ( Figures 9A–E ). For protein expression, three patients were analyzed with IHC. LGGs tissues were obtained from the tumor center (TC) and tumor periphery (TP). We found the expression of five proteins was highly expressed in tumor periphery ( Figure 9F ). After a 4-step grading system was quantified, except for LATS2, other proteins showed high expression in the TP group ( Figures 9G–K ). This founding may be related to the research that tumor cells can influence the recruitment of CAFs precursors and induce the activation of normal fibroblasts into CAFs.