key: cord-0075120-b7lem3c1
authors: Feng, Peng; Li, Yuchen; Tian, Zhijie; Qian, Yuan; Miao, Xingyu; Zhang, Yuelin
title: Analysis of Gene Co-Expression Network to Identify the Role of CD8 + T Cell Infiltration-Related Biomarkers in High-Grade Glioma
date: 2022-02-22
journal: Int J Gen Med
DOI: 10.2147/ijgm.s348470
sha: a3bf1e85b0a7aaf978d69ee7e4f5ef9b24b45449
doc_id: 75120
cord_uid: b7lem3c1

BACKGROUND: High-grade glioma is a type of heterogeneous lethal brain tumor most common in adults. At present, immune checkpoint inhibitors (ICIs) are 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. MATERIALS AND METHODS: The present study adopted a deconvolution algorithm for quantifying immunocyte composition and measuring gene expression, and used weighted gene co-expression network analysis (WGCNA) to analyze GBM expression data obtained from Gene Expression Omnibus (GEO), Chinese Glioma Genome Atlas (CGGA), and the Cancer Genome Atlas Glioblastoma Multiforme (TCGA-GBM) databases. 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. RESULTS: From this study, CD8+ T cell-related modules were identified. By using consistent clustering analysis, two panels of genes (red and green) with the highest correlation with CD8+ T cells infiltration were used to construct high-, low-expression groups, silent and/or mixed group of T cell infiltrations. In the high and low CD8+ T cell infiltration groups, a total of 535 differential genes were obtained, of which ten genes (RPS5, RPS6, FAU, RPS19, RPS23, RPS15A, RPS29, RPS14, RPS16, RPS27A) were identified through protein–protein interactions and co-expression network analysis. Post Cox regression and Kaplan–Meier (K-M) survival analysis, RPS5, RPS6, and RPS16 were selected as candidate prognostic biomarkers related to CD8+ T cells. CONCLUSION: The three associated genes RPS5, RPS6, and RPS16 were markedly related to degree of T cell infiltration and immune-related activated. We identified their potential biomarkers and therapeutic targets associated with the extent of CD8+ T cell infiltration in GBM.

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. 3 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. 4 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. 5 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.

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] [7] [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] [10] [11] [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] [14] [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. 16 In recent years, technology involving bioinformatics has rapidly developed, and numerous approaches have been established for identifying biomarkers. 17 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. 20 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 21 and estimates the degree of immune cell infiltration within each malignant tumor type, like prostate cancer (PCa) 22 or kidney cancer. 23 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.

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.

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. 

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 ).

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-infiltrationrelated 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. 

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.

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 ). 

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 ).

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 ).

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.

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://cytos cape.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) .

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. 24 And the PD-1 pathway could regulate immunological homeostasis and influence the mechanism of autoimmunity protection. 1884 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. 25 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. 28 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. 29 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. 30 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. 31 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.

GBM, glioblastoma multiforme, glioblastoma; CGGA, Chinese Glioma Genome Atlas; TCGA-GBM, The Cancer Genome Atlas Glioblastoma Multiforme; GEO, Gene Expression Omnibus; WGCNA, weighted gene co-expression network analysis; K-M, Kaplan-Meier; ECM, extracellular matrix; CMV, cytomegalovirus; CNS, central nervous system; ICIs, immune checkpoint inhibitors; TME, tumor microenvironment.

Data supporting our findings are already included in the manuscript.

All methods of this study were carried out in accordance with relevant guidelines and regulations. The data of this study are from TCGA, CGGA and GEO public databases, and did not involve animal experiments or human specimens, no ethics-related issues. Ethical approval had been obtained for the patient data in the databases. Users can download relevant data for free for research and publish relevant articles. Our study was based on open source data, so there are no ethical issues.

Every author involved with the manuscript agreed to the publication of the manuscript.

Neuro-oncology in 2015: progress in glioma diagnosis, classification and treatment

Prevalence estimates for primary brain tumors in the United States by age, gender, behavior, and histology

Cytomegalovirus promotes murine glioblastoma growth via pericyte recruitment and angiogenesis

Glioblastoma and other malignant gliomas: a clinical review

PD-1 blockade with nivolumab in relapsed or refractory Hodgkin's lymphoma

Pembrolizumab for the treatment of non-small-cell lung cancer

Recurrent glioma clinical trial, CheckMate-143: the game is not over yet

Microenvironmental regulation of tumor progression and metastasis

Overcoming resistance to checkpoint blockade therapy by targeting PI3Kγ in myeloid cells

The tumor microenvironment underlies acquired resistance to CSF-1R inhibition in gliomas

PI3Kγ is a molecular switch that controls immune suppression

Type, density, and location of immune cells within human colorectal tumors predict clinical outcome

CD8 tumor-infiltrating lymphocytes are predictive of survival in muscle-invasive urothelial carcinoma

Tumor-infiltrating CD8+ lymphocytes predict clinical outcome in breast cancer

Proliferating CD8+ T cell infiltrates are associated with improved survival in glioblastoma

A new risk score based on twelve hepatocellular carcinoma-specific gene expression can predict the patients' prognosis

Characterization of long non-coding RNA and messenger RNA profiles in laryngeal cancer by weighted gene co-expression network analysis

Molecular signature of subtypes of non-small-cell lung cancer by large-scale transcriptional profiling: identification of key modules and genes by weighted gene co-expression network analysis (WGCNA)

WGCNA: an R package for weighted correlation network analysis

Robust enumeration of cell subsets from tissue expression profiles

The immune landscape of prostate cancer and nomination of PD-L2 as a potential therapeutic target

Immune infiltration in renal cell carcinoma

Anti-PD-1 blockade and stereotactic radiation produce long-term survival in mice with intracranial gliomas

Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma

Multicellular gene network analysis identifies a macrophage-related gene signature predictive of therapeutic response and prognosis of gliomas

Single-cell transcriptome-based multilayer network biomarker for predicting prognosis and therapeutic response of gliomas

Inferring microenvironmental regulation of gene expression from single-cell RNA sequencing data using scMLnet with an application to COVID-19

The potential role of ribosomal protein S5 on cell cycle arrest and initiation of murine erythroleukemia cell differentiation

Phosphorylation of ribosomal protein S6 attenuates DNA damage and tumor suppression during development of pancreatic cancer

USP1-dependent RPS16 protein stability drives growth and metastasis of human hepatocellular carcinoma cells

The International Journal of General Medicine is an international, peer-reviewed open-access journal that focuses on general and internal medicine, pathogenesis, epidemiology, diagnosis, monitoring and treatment protocols. The journal is characterized by the rapid reporting of reviews, original research and clinical studies across all disease areas. The manuscript management system is completely online and includes a very quick and fair peer-review system

We acknowledge TCGA, CGGA and GEO databases for providing their platforms and contributors for uploading their meaningful datasets.

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

This study was supported by a grant (No. 2021SF-355) from the Shaanxi Province Key R&D Program.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.