Analysis of Gene Co-Expression Network to Identify the Role of CD8 + T Cell Infiltration-Related Biomarkers in High-Grade Glioma

Introduction

Glioblastoma (GBM) is a cancer condition that can occur at any age, but especially frequently among adult gliomas (approximately 45% of all gliomas).^1,2 However, the prognosis is extremely poor with a median survival rate of 12–15 months.³ A recent study showed that the inhibition of the surveillance of immune response can contribute to the growth of brain cancer tissue. The interaction between glioblastoma cells and extracellular matrix (ECM) has been found to be another key point in the occurrence of glioma. In addition, exposure to radiation causes GBM; long-term exposure to higher doses of non-ionizing radiation is related to the formation of GBM. Finally, viral infections such as cytomegalovirus (CMV)-induced mutations are considered to be a potential cause of GBM formation.⁴ Often, surgery is the first choice of treating patients with GBM, where the tumor tissue is removed with great precision, fully protecting the brain function of patients, combined with radiotherapy, chemotherapy, and electric field treatment with integrated adjuvant treatment measures.⁵ However, GBM, being highly invasive, is difficult to eradicate with surgery, and therefore the patients’ prognosis may be poor and the mortality rate remains high. Hence, new treatment methods need to be designed to control tumor progression and improve patient prognosis.

Table 1 Univariate Cox Proportional Hazards Analyses of Overall Survival of GBM Patients

Table 2 Multivariate Cox Analysis Showing the Hazard Ratios (HR) of Different Factors

In recent years, immune checkpoint inhibitors have played a key role in providing anti-tumor effect by intervening in the interaction of the checkpoint molecular receptors with their ligands, thereby regulating tumor immune response. Further, a lot has been discovered about the molecular mechanisms associated with GBM genesis and progression, the immune system, and the central nervous system (CNS).^6–8 With this background, immune checkpoint inhibitors (ICIs) are now being considered for first-line therapeutics for malignant GBM. Nonetheless, molecular markers for malignant GBM are unavailable at present. As a result, it is important to explore molecular markers related to immunity for GBM.

The strong changes induced in tumor microenvironment (TME) features affect immunotherapy response and play an important role.^9–12 Among the diverse tumor immune cells, CD8+ T cells are known to be associated with the adaptive immunity of tumors occupying the greatest fraction and function. In many solid tumors, a high degree of CD8 + T cell infiltration was found to be the main factor for tumor immunotherapy.^13–15 Immune biomarkers for high-grade gliomas have been analyzed extensively, yet it is not possible to use these observed candidate factors in clinical practice directly. Therefore, the identification of biomarkers related to the immune infiltration mechanism of CD8+ T cells may contribute to monitoring the response to GBM immunotherapy.

A recent study showed that a higher CD8/CD4 ratio is a predictor of overall survival in glioblastoma patients, and also the proliferation of CD8+ T cells could enhance survival, and measures may be taken to promote the dominance of CD8+ T cells in GBM to interfere with the immunoregulatory effects of IFNγ in the tumor microenvironment.¹⁶

In recent years, technology involving bioinformatics has rapidly developed, and numerous approaches have been established for identifying biomarkers.¹⁷ Weighted gene co-expression network analysis (WGCNA) has been developed as an efficient approach for mining gene correlations in identifying pivotal anticancer genes and related modules.²⁰ This algorithm approach has been extensively used for identifying RNA transcripts.^18,19 Typically, a cell type is identified by estimating relative subsets of RNA transcripts. Similarly, CIBERSORT is a bioinformatics approach for the analysis of gene expression profiles that utilizes the deconvolution algorithm for quantifying immunocytes²¹ and estimates the degree of immune cell infiltration within each malignant tumor type, like prostate cancer (PCa)²² or kidney cancer.²³

To explore the functional role of TME and identify a suitable biomarker for GBM, this study adopted WGCNA based gene expression profiles in GBM and CIBERSORT algorithm for calculating the composition of T cells in samples. Thereafter, key CD8+ T cell infiltration-related genes and modules were identified, and database analysis was conducted to verify the therapeutic and immune features of the selected genes. Through this study, we also discovered and validated candidate prognostic biomarkers.

Materials and Methods

Sources of Gene Expression Profiles

The mRNA expression profiles were obtained from 1092 GBM samples in Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/), CGGA, and TCGA databases. A total of 1092 patients’ mRNA expression with detailed survival information was download from TCGA-GBM from TCGA database (https://portal.gdc.cancer.gov/), CGGA dataset and GEO datasets including GSE4412 and GSE7696. After merging of multiple data-sets across platforms, the normalization of data was required. Expression values were log-transformed and meanwhile we used the “ComBat” algorithm to reduce the possibility of batch effects due to non-biotech bias between different data-sets.

Evaluation of Tumor-Infiltrating Immune Cells

CIBERSORT, an analytical algorithm used to analyze RNA expression profiles for assessing diverse cell type proportions from every sample, showed 22 TIIC scores. And then, we used the 22 TIIC scores as the trait data of WGCNA.

Co-Expression Network Construction

The R package “WGCNA” was used to build the co-expression network considering the levels of 4187 genes. Firstly, a single transcript expression profile was transformed to a similarity matrix according to Pearson’s correlation coefficients of matched genes. Subsequently, the similarity matrix was transformed into an adjacent matrix, like amn = |cmn|β (amn = paired gene adjacency; cmn = paired gene Pearson’s correlation). Typically, β enhanced the potent gene-gene correlation while reducing poor correlations. To construct a scale-free network, this study selected β = 7 to be a soft thresholding power (Figure 1A and B). At this β, the adjacency matrix was converted into a topological overlap matrix. Thereafter, the hierarchical clustering tree was built by dynamic hybrid cutting, where one leaf stood for one gene. To classify genes whose expression profiles were close in diverse modules, a dynamic hybrid cutting approach was adopted for which the bottom-up algorithm was used, and a minimal module size threshold was set to be 30. The focus of the study was to show red and green modules associated with CD8+ T cells and to identify the two central modules (Figure 2A).

Figure 1 Selection of appropriate beta value to construct hierarchical clusters by (A) analysis of the unscaled fitting index of 1–20 soft threshold power (β). (B) Analysis of the average connectivity of 1–20 soft threshold powers. (C) Classification of genes into different modules through hierarchical clustering (different colors are used to represent different modules).

Figure 2 (A) The heat map showing a correlation between characteristic genes of the module and infiltration of T cells. Enrichment analysis of genes in (B) green pane. (C) red panel.

Genomic Variation Analysis (GSVA) and Functional Annotation

To explore differences in BPs in the previously mentioned four groups, the R package “GSVA” was adopted. Based on MSigDB database, this study acquired “c2.cp.kegg.v6.2.symbols” for GSVA. Adjusted P<0.05 value stood for statistical significance. Group-related genes were later annotated by R package cluster Profiler upon the threshold of FDR<0.05.

We included the use of the gray and brown plates as T CD8+ high-infiltration-related genes and low-infiltration-related genes, respectively, and used consensus ClusterPlus (parameters: rep¼100, pItem¼0.8, pFeature¼1) for clustering. Euclidean distance and D2 were used to measure distance and clustering algorithm separately.

Results

Tumor-Infiltrating Immune Cells (TIIC) Assessment

CIBERSORT, an analytical algorithm used to analyze RNA expression profiles for assessing diverse cell type proportions from every sample, showed 22 TIIC scores. And then, we used the 22 TIIC scores as the trait data of WGCNA.

GBM Gene Co-Expression Network

The co-expression network was built based on 4187 gene expression profiles by R package “WGCNA”. Thereafter, the mean Pearson’s correlation and linkage coefficients were determined for clustering 1092 samples, while genes showing close expression profiles were clustered for forming a tree branch, which stood for one gene module. To build a scale-free network, we picked β = 8 (scale free R2 = 0.87676) as the soft-thresholding power (Figure 1A and B). There were eight modules formed (Figure 1C).

Identification and Enrichment Analysis of Hub Modules

Among these eight modules, we selected CD8+ T cells as candidate markers and therefore focused on showing red and green modules associated with CD8 + T cells (Figure 2A). The green module showed a strong negative correlation with CD8+ T cells (R2 = 0.41, P = 1e-44), and the red module showed a strong positive correlation with CD8 + T cells (R2 =−0.35, P = 1e−33). Among them, the red section showed a total of 81 genes, and the green section showed a total of 128 genes. Many genes observed from the red section were found to be related to pathways that regulate immunity, like the pathway modulating stem cell pluripotency, the ErbB pathway, and endocrine resistance (Figure 2C). The genes in the green section were enriched with pathways that negatively regulate immunity, such as the MAPK signaling pathway and the Hippo signaling pathway (Figure 2B).

Immunity Section Grouping

The median expression level of co-expressed gray panel genes and the median expression level of brown panel genes were used to assign groups that showed CD8+ T cells with high infiltration group (gray panel gene>0, brown panel gene <0) and CD8+ T cells with low infiltration group (gray panel gene <0, brown panel gene>0), mixed group (gray panel gene>0, brown panel gene>0) and silent group (gray panel gene<0, brown panel gene<0) (Figure 3A–C). The CD8+ T cells were found to belong to the high infiltration group, which was further verified for acceptable infiltration levels (Figure 3D).

Figure 3 (A) The heat map depicting a consistent clustering solution (k = 3) of related genes represented in red and green blocks in the GBM sample (n = 1092). (B) Scatter plot showing the median expression level of a red plate (x-axis) and green plate (y-axis) genes co-expressed in each GBM sample with subgroups assigned based on the relative expression levels of related genes in red and green panels. (C) The heat map depicting the expression levels of genes co-expressed in red and green panels in each subgroup. (D) The distribution of immune cells in each subgroup.

Grouped Enrichment Analysis

The biological behavior between different groups by GSVA (Figure 4) showed the green plate to be significantly enriched with oxidative phosphorylation and ribosome. The red plate showed the enrichment pathways related to complete immune suppression, including the TGF beta, ERBB, and MTOR pathways.

Figure 4 (A–F) The biology function analysis between different subtypes.

Differential Gene Analysis and Prognostic Analysis of High and Low Infiltration Groups

Analysis of differential genes by Limma, between high and low infiltration groups to identify key genes that affect CD8+ T cell infiltration, showed altogether 535 differential genes (Figure 5A). The search approach (STRING; http://string-db.org) to identify interacting genes for construction of PPI network was based on co-regulated pivot genes and evaluated significant differences between these DEG-encoded proteins. Interaction and Cytoscape (v3.0, https://cytoscape.org/) software for identifying and visualizing interaction between these proteins along with cytoHubba software were used to screen the top 10 core genes (Figure 5B). Following that, we extracted these 10 core genes and tested them for survival and prognosis analysis, wherein we found that a total of 3 genes showed survival differences, including RPS5, RPS6, and RPS16 (p<0.5, Figure 5C–E). In COX analysis, these genes were shown to become independent prognostic factors (Tables 1 and Tables 2).

Figure 5 (A) Differential gene analysis of the red and green panel group. (B) Screening of core genes from differential genes. (C–E) Survival curves of RPS5, RPS6, and RPS16.

Discussion

The immune system of a human body is equivalent to a “scavenger” in its functioning. During the normal physical condition, it can identify and remove cancer cells from the TME. However, the malignant tumor cells bypass this effect and can develop into a variety of types under the supervision of the same immune system, thereby escaping the attack strategies of the immune system, which include the various stages of the anti-tumor immune response.

GBM is a frequently occurring aggressive adult cancer condition which, like other tumors, genetically alters the host causing a diseased state due to the loss of control over growth and inadaptability of the human body (including the effects of the immune system).

The current treatment for GBM includes radiotherapy and chemotherapy, surgical treatment, and immunotherapy. At present, immunotherapy is a promising treatment method which can effectively improve the prognosis of patients with GBM. The efficacy of immunotherapy largely depends on the expression of PD-1. The interaction of PD-1 and PD-L1 (especially in the T cell) can promote invasion of glioblastoma multiforme by creating an immunoregulatory axis.²⁴ And the PD-1 pathway could regulate immunological homeostasis and influence the mechanism of autoimmunity protection. The PD-L1 expression of GBM can activate the PD-1 receptor, leading to the negative regulation of T cells, especially CD8+ T cells. A recent study showed that the PD-L1 expression of a tumor can be activated by activation of various receptors such as toll like receptor (TLR), epidermal growth factor receptor (EGFR), interferon alpha receptor (IFNAR), and interferon-gamma receptor (IFNGR). A recent study identified that the PD-L1 expression in T cells was correlated with WHO grading and could be identified as a tumor biomarker.

Notably, at the tumor site, there are several inhibitory signals which frequently lead to dysfunction of tumor-specific T cells. Selecting such signals via immune checkpoint inhibitors (ICI) may rejuvenate T cells and improve clinical efficacy. ICIs have been successfully applied in GBM and have attracted great interest in exploring immune factors for immunotherapy.²⁵ Of these, CD8+ T cells serve as the main mediator for anticancer immunity in modulating their response, which has been the main focus of several cancer immunotherapy strategies. After recognition of specific antigen peptides presented via human leukocyte antigen (HLA, in humans) or major histocompatibility complex (MHC, in vertebrates), CD8+ T cells have been shown to be activated for killing cancer cells. Studies extensively report their important function in immunotherapy. Thus, for tumor patients, it is important to identify the prognostic factors and the target of immunotherapy to improve survival rate. Currently, single cell sequencing and analysis have been widely used in various cancer analyses. Especially in GBM, single cell analysis can be used to identify the key prognostic factors^26,27 and understand the tumor microenvironment.²⁸

In this study, we used the cibersort and WGCNA analysis to identify the related genes of CD8+ T cell infiltration, and green and red modules were marked as the positive and negative genes with CD8+ T cell infiltration. Our study showed that the green module genes were correlated with oxidative phosphorylation and ribosome, the red module genes were enriched with the TGF beta, ERBB, and MTOR pathways. With the two algorithms, we were able to accuratelyidentify the related genes of CD8+ T cell infiltration. In order to confirm the high/low CD8+ infiltration samples, we used the unsupervised clustering analysis to identify different subtypes of CD8+ T cell infiltration which actuated different prognosis. A total of four subtypes were identified, which showed significant differences in survival, indicating their potential to promote the progression of GBM. Based on that, we identified the green subtype was high CD8+ T cell infiltration and activated immune-related pathways and the red was marked as cold immune response. Then we used the GSVA analysis to understand the underlying biology function, and the result showed that the red cluster was significantly enriched with glycolysis and immune pathways such as the B-cell and T-cell receptor pathways and green plate showed the enrichment pathways related to complete immune activation, including the TGF beta, MAPK, and MTOR pathways, which is consistent with our classification. To obtain the key genes, we used the limma packages to confirm 535 differential genes between red and green groups, and 10 core genes (RPS5, RPS6, FAU, RPS19, RPS23, RPS15A, RPS29, RPS14, RPS16, RPS27) were obtained by using a protein interaction network. The expressions of RPS5, RPS6, and RPS16 showed visible differences in survival, and also, the COX analysis of these genes showed that these genes are independent prognostic factors. The deregulation of RPS5 gene expression (constitutive expression) affects RPS5 protein level and delays both the onset of initiation of erythroid maturation and entrance in cell cycle arrest in inducer-treated MEL cells.²⁹ Ribosomal protein S6 (rpS6), a component of the 40S ribosomal subunit, has been reported to be associated with various physiological and pathophysiological functions. A recent study showed that phosphorylation of rpS6 attenuates KRAS-induced DNA damage and p53-mediated tumor suppression.³⁰ The transfection-enforced re-expression of RPS16 restores oncogenic-like activity in USP1-deficient HCC cells. Importantly, the high expression of USP1 and RPS16 in liver tissue is a prognostic factor for poor survival of HCC patients.³¹

Our research showed that the green subtypes had high immune cell infiltration and were enriched with oxidative phosphorylation and ribosome. Between the green subtype and red subtype, we found that RPS5, RPS6, and RPS16 had significant expression in the green subtype and could be prognostic factors of GBM. As the core genes of green subtypes, RPS5, RPS6, and RPS16 are involved in CD8+ T cell infiltration regulation and the activation of oxidative phosphorylation and ribosome pathways. Altogether, the three associated genes RPS5, RPS6, and RPS16 were markedly related to degree of T cell infiltration and immune-related activity. Our study has one limitation. In our study, we investigated these genes which may be potential biomarkers and targets for GBM immunotherapy, with certain notable limitations in the study. The mechanism and validation of these genes still need to be further researched in clinical and molecular biology experiments.