Immune Infiltrating Cells-Derived Risk Signature Based on Large-scale Analysis Defines Immune Landscape and Predicts Immunotherapy Responses in Glioma Tumor Microenvironment

The following search terms were used as: (((survival OR prognosis OR prognostic OR outcome OR death OR relapse OR recurrence))) AND (Glioma[Title]) OR (Astrocytoma*[Title]) OR (Glioblastoma*[Title]) OR (Ependymoma*[Title]) OR (Oligodendroglioma*[Title]) OR (Gliosarcoma*[Title]) OR (Astroglioma*[Title]) OR (LGG[Title]) OR (HGG[Title]) OR (glial cell tumor[Title]). All datasets were manually examined. Patients lacking survival information were excluded from further evaluation.

Publicly available glioma gene-expression datasets together with clinical annotations were downloaded and examined. 2877 samples from 12 patient cohorts diagnosed with gliomas were included in this study (Table S1). Four external datasets were included for validation: GSE13041↗, GSE16011↗, GSE61335↗, GSE68838↗. The microarray datasets were downloaded from the Gene-Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/↗). Chinese Glioma Genome Atlas (CGGA) datasets were downloaded from the CGGA website (http://www.cgga.org.cn/↗), while The Cancer Genome Atlas (TCGA) datasets were downloaded from UCSC Xena (https://xenabrowser.net/↗).

Raw data from microarray datasets were generated using Affymetrix and Agilent. The RMA algorithm was used to perform quantile normalization and background correction of the raw data from Affymetrix in the Affy software package. The final summarizing of oligonucleotides for each transcript was based on the consensus median polish algorithm in the Affymetrix software. The raw data from Agilent was processed using limma software. RNA-sequencing data were downloaded from the TCGA and CGGA data portals and the fragments per kilobase million (FPKM) values were transformed into transcripts per kilobase million (TPM) values, which were similar to those resulting from microarrays and comparable between samples. The TPM values from TCGA and CGGA had similar signal intensity with the RMA-standardized values from microarray datasets. R package sva was then used to remove the computational batch effect.

The CIBERSORT algorithm was used to predict the presence and quantify the abundance of immune cells in glioma samples (11). The LM22 gene signature was applied since it allowed for sensitive and specific discrimination of 22 human infiltrating immune cell phenotypes. Gene-expression profiles were uploaded to the CIBERSORT web portal (http://cibersort.stanford.edu/↗). This algorithm was run using the LM22 signature and 1000 permutations. Single factor analysis was performed on the 22 immune cells to determine their prognostic values in gliomas. The cellular correlation among the 22 immune cells was performed using Pearson correlation analysis. TIMER algorithm (9), EPIC algorithm (12), MCPcounter algorithm (13), quanTlseq algorithm (14), xCell algorithm (15), and ssGSEA algorithm (16) were also used for estimating the abundance of immune infiltrating cells.

Tumors with qualitatively diverse TME infiltrating patterns were classified using Partitioning Around Medoid (PAM) (17), which identified TME patterns and grouped patients for further analysis. The optimal number of clusters and their stability and reliability in the meta-cohort and TCGA cohort were determined using the ConsensuClusterPlus R package. Infiltration level of stromal cells and immune cells in glioma samples was assessed using the consensus ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression) algorithm (18).

To identify genes associated with TME cell-infiltrating patterns, the patients were grouped into two distinct TME clusters based on the diverse expression of infiltrating immune cells. The enrichment levels of immune signatures were quantified by the xCell algorithm to validate TME clusters (15). The R package limma (19) was used to determine DEGs associated with the two TME cell-infiltrating patterns. The adjusted P-value < 0.01 was used to determine DEGs among the TME subtypes.

The DEGs among the TME clusters were standardized in all samples in the TCGA glioma cohort. Univariate cox regression analysis identified representative DEGs. The unsupervised clustering method (20) was used to classify patients into two TME gene clusters for further analysis. The clusterProfiler R package (21) was used to annotate the TME pattern genes. The consensus clustering algorithm (22) was performed to define the gene clusters. χ² contingency test was used to determine the correlation between the TME gene clusters. The SVM-RFE algorithm was used for dimension reduction and to mitigate the interference effect of redundant genes (23). The top 300 DEGs between two TME gene clusters were selected (24), among which 63 most representative genes were identified with the highest accuracy of separating samples. Principal component analysis (PCA) was performed and the extracted principal component 1 served as the signature score. A method similar to GGI was then applied (25) to define the TME-score of each patient after the prognostic value of gene signature score was obtained:

where i is the signature score of clusters with HR>1, and j represents the expression of genes with HR<1.

All gene sets were downloaded from the MSigDB database. Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) were performed on the TME gene signatures using the clusterProfiler R package and GSVA R package (21). Pathways enriched in TME immune-related gene patterns were identified in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) with the false discovery rate (FDR) < 0.05 and a strict cutoff of P < 0.01.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to infer individual responses to immunotherapy such as immune checkpoint blockade (e.g. anti-PD-1 therapy) (26). The submap analysis was applied to compare differences in response to anti-PD-1 and anti-CTLA-4 therapies. For the melanoma data set (GSE78220↗, N=28), GSE78220↗ expression profiles (FPKM normalized) were transformed into TPM values, which were used to calculate the TME-score (27). T cell-inflamed gene expression profile (GEP) was defined through the expression of the 18 genes (28). Cytotoxic activity (CYT) was determined based on the gene expression value of two cytolytic markers (GZMA and PRF1) (29), and the geometric mean of these two markers was used to perform the calculations. Seven types of immune checkpoints were collected from previously published work (30).

For each glioma patient of the 48 samples, major exclusion criteria were incomplete follow-up data, poor quality of samples, and missing baseline clinicopathological features. RNAstore-fixed tumor tissues were then collected for sequencing. Briefly, 1 μg RNA per sample was used as input material for RNA sample preparations. RNA was extracted and sheared followed by sequencing library preparation using NEBNext^® UltraTM RNA Library Prep Kit. Subsequently, PCR was performed with Phusion High-Fidelity RNA polymerase, Universal PCR primers and the Index (X) Primer. After PCR primer removal, biotin-labeled probe was used for capturing target regions. The captured libraries were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated. Raw data (raw reads) of fastq format were first processed through in-house perlscripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, ploy-N, and low-quality reads from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All downstream analyses were based on clean data with high quality. Reference genome and gene model annotation files were downloaded from the genome website directly. The reference genome index was built using Hisat2 v2.0.5 and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5 and Hisat2 was selected as the mapping tool. FeatureCounts v1.5.0-p3 was then used to count the reads numbers mapped to each gene. TPM of each gene was calculated based on the gene length and reads count mapped to this gene.

Immunohistochemistry (IHC) staining was conducted as previously described (31, 32). Paraffin-embedded tissues of 40 glioma samples with the corresponding sequencing data from the Xiangya Neurosurgey (XYNS) cohort were used for performing IHC. The paraffin-embedded glioma sections were incubated with CD68, CD163, CD117 (Rabbit, 1:500, Proteintech, China), and CD66b (Rabbit, 1:200, Abcam). The IHC marker was detected with microscope.

The normality of variables was tested using the Shapiro-Wilk normality test (33). For normally distributed variables, unpaired Student’s t-test was used to compare differences between two groups, while the Wilcoxon test was used to compare nonnormally distributed variables. For multiple groups, one-way analysis of variance was used as a parametric method to compare mean values between groups while Kruskal–Wallis tests were used as a nonparametric method.

Pearson analysis and distance correlation analyses were used to calculate correlation coefficients. Contingency tables were analyzed by χ² contingency test. The OS and TME-score were calculated using the R package survival and cutoff values determined. Based on the dichotomized TME-score, patients were grouped into high or low TME-score in each data set, and the computational batch effect was reduced by the R package sva. Data were mainly visualized using the R package ggplot2. R package, igraph was used to visualize the cellular interactions within the TME. For the differential gene expression analysis, we used the Benjamini–Hochberg method that converts the P values to FDRs to identify significant genes (34). The package pROC (35) was utilized to establish receiver operating characteristic (ROC) curves and calculate the area under the curve (AUC). OncoPrint was used to delineate the mutation landscape of TCGA via the maftools R package (36). The Kaplan–Meier method was applied to generate and visualize survival curves for the subgroups, and the differences between data sets were compared using the log-rank test. The univariate and multivariate Cox proportional hazards regression models were utilized to calculate the hazard ratios in univariate and multivariate analyses and to determine independent prognostic factors using the R package survival. All survivorship curves were generated via R package survminer. All heatmaps were generated based on pheatmap. All statistical analyses were performed in R version 3.6.1 (https://www.r-project.org/↗). All tests were two-sided and P values < 0.05 were considered to be statistically significant.

The TME cell infiltration patterns and gene signatures were evaluated (Figure 1). The ConsensusClusterPlus package was used to assess clustering stability to determine the optimal cluster number (Figure S1A), which supported two robust subtypes of gliomas in a meta-cohort. The optimal cluster number was also identified in TCGA (Figure S1B). We built an integrated TME cell network that comprehensively depicted a landscape of TME cell lineages, tumor immune cell correlations, and their prognostic values on the OS of glioma patients (Figure 2A and Tables S3, S4). A similar TME cell network was constructed in the TCGA cohort (Figure S2A and Tables S3, S4). Partitioning Around Medoid (PAM) was performed for the 2877 tumors with the corresponding TME cell expression profiles of the 12 included glioma cohorts (Figure 2B and Table S2). PAM was subsequently performed in the TCGA cohort (1027 patient samples), and two phenotypes were separated by different clinical factors (Figure S2B and Table S6). Two TME phenotypes were identified by TME cell infiltration and exhibited significant differences in the OS in the meta-cohort and TCGA alone (log-rank test, p < 0.001; Figures 2C and S2C). Among the two clusters, PCA distribution was separated in both the meta-cohort and TCGA alone (Figures 2D and S2D). Figure 2E shows that the two TME clusters exhibited significant differences in TME cell infiltration patterns, and these differences were reported in multiple immune suppressive cell types. Additionally, PAM was performed in another two cohorts for validation, and there was a significant correlation between identified TME phenotypes and TME infiltration cell patterns in the CGGA and GSE108474↗ datasets, respectively (Figures S2E, F and Table S6). We further analyzed immune cells by the xCell algorithm (15). The two TME clusters identified by xCell algorithm were significantly different in terms of the survival probability (Figure S3E). Additionally, the TME clusters identified by xCell algorithm were consistent with those identified by the CIBERSORT algorithm (Figure S3F).

To elucidate the correlation between the immune infiltrating environment and TME clusters, 20 of both immune-related and DNA regulation-related signaling pathways in GO analysis were identified in the meta-cohort (). We found that TME cluster1 was associated with immunosuppressive pathways (). Additionally, TME cluster 1 was associated with pathways regulating tumor cell proliferation (). Similar results were observed in the TCGA database (), showing differences in 20 signaling pathways in the two TME clusters (). 2 2 2 2 2

We acquired a total of 1312 DEGs (Table S5) using the limma package to classify the patients into genomic subtypes and to investigate the potential biological characteristics of each TME infiltration cell pattern. The analysis was significantly consistent with the clustering results of the TME phenotype groups (χ² contingency tests, p =1.95 × 10^-12). The TCGA glioma cohort population was grouped into two TME gene clusters 1 and 2 (Figure S4A). The survival analysis of the two patient clusters indicated that gene cluster 1 correlated with worse survival outcomes than cluster 2 (Figure S4B). The GO and KEGG enrichment analyses showed that gene clusters 1 and 2 were enriched in distinct biological processes. In GO enrichment analysis, gene cluster 1 was involved in tumor proliferation (Figure S4C). Overexpression of genes involved in immune activation pathways was enriched in gene cluster 2 (Figure S4D). Additionally, the KEGG enrichment analysis showed that gene cluster 1 was associated with tumor proliferation and was a prognostic marker for poor survival outcomes (Figure S4E). Gene cluster 2 was associated with immune activation and was a prognostic marker for better survival outcomes (Figure S4F).

The SVM-RFE algorithm was used in dimension reduction to extract phenotype signatures with high classification accuracy and further explore the role of TME phenotypes. Sixty-three most representative DEGs were identified (and); the chromosomal distribution and expression of these genes are displayed in. Almost all of the 63 genes were significantly differentially expressed between glioma molecular subtypes, isocitric dehydrogenase (IDH) mut glioma and IDH wt glioma (). The regulatory networks identified by the clusterProfiler R package suggested that immune activation and tumor proliferation pathways were involved and exhibited significant overlaps with other pathways (). 2 1 2 1 2

The PCA algorithm was used to define the TME-scores of the 12 cohorts (Table S8). Based on the 63 DEGs, PCA distribution was separated among the two TME gene clusters in TCGA (Figure S5F). The interconnections among TME clusters, TME gene clusters, TME-scores, patient survival, and tumor grade are presented in Figure S5D. Contingency table revealed the significant consistency between TME clusters and TME-scores, which TME score could be considered a collection of the features of the two TME clusters (Figure S5E). The distribution of TME-scores in TME clusters in CGGA and GSE108474↗ datasets are shown in Figure S5G. Consistent with the findings in TME gene cluster 1, a high TME-score was a prognostic marker for poor clinical outcomes in TCGA and meta-cohort (Figures 3A, B). Given that gliomas consist of various types and grades of glial tumors, the differences in TME landscapes among different types of glial tumors were explored. The prognostic values of TME-scores were verified in TCGA LGG (Figure S6A) and TCGA GBM (Figure S6B). The TME-score was also a prognostic marker for the IDH status in LGG (Figure S6C), GBM (Figure S6D), and pan-gliomas (Figure S6E). A high TME-score was associated with metastatic and immunosuppressive signatures (Figure 3C). The correlations between TME-scores and these known signatures are shown in Figure S5H.

The associations between TME-scores and the immune infiltrating environment was further examined. High TME-scores were correlated with Estimate Scores, Immune Scores, and Stromal Scores (Figure 3D) and also associated with the infiltration of M2 macrophages, mast cells, and neutrophils. Thus, high TME-scores were an indication of the immunosuppressive environment and poor survival outcomes while low TME-scores were prognostic for activated immune environments (Figures 3E, S5I). Further, TME-scores were significantly correlated with CD8 T cell, NK cell, regulatory T cells (Tregs), macrophages, fibroblasts, Th1 cells, and dendritic cells based on TIMER algorithm, EPIC algorithm, MCPcounter algorithm, quanTlseq algorithm, and xCell algorithm (Figure 3F). The immunocyte infiltrating characteristics of TME-scores were verified in LGG (Figure S6F) and GBM (Figure S6G). Gliomas with a high TME-score expressed more immune checkpoints, such as LAG3, CD40, and PDCD1LG2 (Figure 3G). Figure S7B displays the expression differences of TME-scores in relation to several clinical factors. Gliomas with unmethylated MGMT, IDH WT, 1p19q non-codeletion, higher grade, and mesenchymal gliomas with poor survival outcomes had high TME-scores.

We performed a copy number variant (CNV) (Figure 4A) and somatic mutation analysis (Figures 4B, C) of the TCGA dataset to determine the associations between TME-score and glioma genomic profiles. Various frequently amplified and deleted genomic regions were observed in high and low TME-score samples (Figure 4A). The somatic mutation analysis showed mutation frequencies of various genes as follows: TP53 (45%), IDH1 (40%), ATRX (27%), TTN (21%), and EGFR (20%) in the high TME-score (Figure 4B), while IDH1 (62%), TP53 (31%), and CIC (27%) in the low TME-score cluster (Figure 4C). A CNV comparison between high and low TME-score clusters identified significantly different CNV regions (Table S11). The mutated genes were compared between high and low TME-score clusters, identifying different mutated genes (Table S12).

High TME-scores demonstrated significant enrichment of PD-1 signaling, T cell signaling, Hypoxia signaling, exosome signaling, immunosuppressive cells including Tregs, myeloid-derived suppressor cells (MDSCs), tumor-associated macrophages (TAMs), and cancer-associated fibroblasts (CAFs) (Figure 5A). m6A signatures that were associated with antitumor immunity were also enriched in high TME-scores (Figure 5B). In term of antigen presentation capacity, high TME-scores presented higher antigen processing and presenting machinery (APM) score (Figure 5C). Cancer testis antigen (CTA) and neoantigens were vital source of tumor-specific antigens, and they were significantly higher in high TME-scores (Figures 5D, E). A series of factors associated with tumor immunogenicity was then assessed. High TME-scores presented higher level of nonsilent mutation rate, number of segments, and aneuploidy score, all of which were significant indicators for genome alteration (Figures 5F–H). Stroma signatures including TGF-beta response, leukocyte fraction, and stromal fraction were higher in high TME-scores (Figures 5I–K). Intratumor heterogeneity and tumor mutation burden (TMB) predicting better immunotherapy responses were also higher in high TME-scores (Figures 6A, B).

We assessed the TME-scores in glioma cohorts because of their prognostic significance associated with poor outcomes in glioma datasets (Figure 4D). Survival analysis in the 10 included cohorts indicated an association of the high TME-score with poor survival outcomes in all datasets (Figure S7A). Univariate and multivariate Cox regression models in both TCGA and CGGA cohorts showed that the TME-score model was an independent prognostic factor (Figure S8A). TME-scores were also validated in several external datasets and high TME-scores were found to be prognostic markers for poor survival outcomes (Figure S9A).

The ability of the TME-score to predict patients’ response to immune-checkpoint therapy was explored by assigning the GSE35640↗ cohort patients (melanoma dataset) to different TME-score groups. Patients with high TME-scores exhibited better immunotherapeutic responses (Figure 6E). The ROC analyses confirmed that TME-score was a predictive biomarker in patients with melanoma (Figure 6F). In another melanoma dataset, GSE78220↗, patients with high TME-scores exhibited significantly longer OS compared to patients with low TME-scores (Figure 6G). High TME-scores also correlated with complete anti-PD-1 responses (Figures 6H, I). The expression patterns of TME-scores in 27 melanoma patients with complete anti-PD-1 and partial anti-PD-1 responses and progressive disease are displayed in Figure S8C. The ROC analyses confirmed that TME-score was a predictive biomarker in melanoma patients (Figure 6J). TME-score was also found to be significantly correlated with two classical immune checkpoint molecules, LAG3 and PDCD1LG2 (Figure S8F).

To further elucidate the correlation between TME-score and immunotherapy, the potential response to immunotherapy in TCGA based on the TIDE algorithm was evaluated. Patients with high TME-scores exhibited better immunotherapy response compared to those with low TME-scores (Figure 6C). Subsequently, responses to anti-PD-1 and anti-CTLA-4 therapies were analyzed. The results showed different responses between high and low TME-score groups to both immunotherapies, which a high TME-score exhibited a significant response to anti-PD-1 immunotherapy while TME-score predicted no response for anti-CTLA-4 immunotherapy in TCGA (Figure 6D). Further, TME-scores were significantly associated with CYT (Figures 6K, L). A high TME-score indicated increased expression of GEP (Figures 6M, N). TME gene cluster 1 and TME phenotype cluster 1 also showed a high expression of CYT and GEP (Figures S8D, E).

The potential associations between TME-scores and signaling pathways in GO and KEGG pathways based on GSVA were analyzed in TCGA. GO results showed that high TME-scores were significantly associated with immune-related pathways (). KEGG analysis showed that a high TME-score was associated with pathways in cancer, apoptosis, and VEGF signaling pathway (). These results denoted the complexity of TME and also showed that activated T cells were major components of immune infiltrating cells. GSEA indicated that negative regulation of the immune response and T cell activation were enriched in high TME-scores (). Pathways in cancer were associated with high TME-scores as shown in KEGG analysis (). A high TME-score indicated the presence of higher overall variants (). Additionally, a high TME-score was associated with lower arm SCNA levels and higher chromosome SCNA levels in this study (). 2 2 2 2 2 2

We validated TME-scores in our samples with sequencing data from Xiangya hospital (Table S14). High TME-scores served as prognostic markers for poor clinical outcome in the XYNS cohort (Figure 7A). Figure S10A shows that the high TME-score exhibited significant higher TME cell infiltration level based on ssGSEA algorithm, and these differences were reported in multiple immune suppressive cell types, including macrophages, mast cells, MDSCs, and Tregs. Gliomas with a high TME-score expressed higher levels of immune checkpoint molecules including CD274, CD276, PDCD1LG2, LAG3, PDCD1, TIGIT, IDO1, CTLA4, and TGFB1 (Figure S10B). As shown in Figure 7B, TME-scores were also found to be positively correlated with four classical immune checkpoints (CD274, PDCD1LG2, LAG3, and PDCD1). Given the critical role of macrophages in the tumor microenvironment, CD68 and CD163, specific M2 macrophage markers, were used for the identification of M2 macrophages in the glioma microenvironment (37). IHC staining showed that high TME-scores exhibited a relatively higher expression of CD68 and CD163 compared to low TME-scores (Figure 7C), implying a high infiltration of M2 macrophages in the tumor microenvironment with high TME-scores. Besides, IHC staining was also performed for CD66b (marker of neutrophil) (38) and CD117 (marker of mast cell) (39) (Figure S11). Likewise, high TME-scores exhibited a relatively higher expression of CD66b and CD117 compared to low TME-scores, implying a high infiltration of neutrophils and mast cells in the tumor microenvironment with high TME-scores.