Identification of a Signature for Predicting Prognosis and Immunotherapy Response in Patients with Glioma

Glioma is a deadly tumor that accounts for the vast majority of brain tumors. Thus, it is important to elucidate the molecular pathogenesis and potential diagnostic and prognostic biomarkers of glioma. In the present study, gene expression profiles of GSE2223 were obtained from the Gene Expression Omnibus (GEO) database. Core modules and hub genes related to glioma were identified using weighted gene coexpression network analysis (WGCNA) and protein-protein interaction (PPI) network analysis of differentially expressed genes (DEGs). After a series of database screening tests, we identified 11 modules during glioma progression, followed by six hub genes (RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1) that can predict the prognosis of glioma and were validated in glioma tissues by qRT-PCR. The CIBERSORT algorithm was used to analyze the difference of immune cell infiltration between the glioma and control groups. Finally, Identification VSIG4 for immunotherapy response in patients with glioma demonstrating utility for immunotherapy research.


Introduction
Glioma is a deadly tumor that accounts for about a third of all brain tumors [1]. ey originate from the glial cells in the central nervous system and have a more malignant histological appearance compared to other brain tumors [2]. e World Health Organization classifies glioma into low and high grades. Low-grade glioma (LGG) (astrocytoma, oligodendroglioma, and oligoastrocytoma) is a welldifferentiated tumor with high 5-year survival rates of almost 60% [3][4][5]. Among the high-grade gliomas, glioblastoma (GBM) is the most aggressive tumor with a median survival of only 15 months [6]. About 70% of LGG patients progress to a high grade within 5-10 years [7].
Glioma usually occurs relatively late in life, which makes complete resection difficult [8]. Current treatment for glioma includes surgical removal, followed by chemotherapy and radiation. Aggressive surgical treatment is associated with serious side effects, including the development of resistance against chemotherapy and radiotherapy, leading to treatment failure, tumor recurrence, and ultimately death [9]. erefore, it is necessary to identify new diagnostic and prognostic biomarkers for personalized gene and molecular therapies.
Weighted gene coexpression network analysis (WGCNA) is a systems biology method used to identify the correlations between genes in microarray samples, as well as identify modules of highly correlated genes [10]. WGCNA can be used to discover genes and biological processes of unknown function, candidate diseases, or transcriptional regulatory work. Although WGCNA cannot prove a causal relationship, coexpression networks can identify regulatory genes of different phenotypes. e network approach bridges the gap between individual genes and systemic tumors [11,12].
In this study, the GSE2223 dataset, including 50 glioma samples and four control samples, were used to conduct WGCNA. Gene coexpression networks and gene modules were identified. Overlapping genes were evaluated using the black-blue module and protein-protein interaction (PPI) network of differentially expressed genes (DEGs) using the cytoHubba plugin to identify diagnostic and prognostic related hub genes of glioma. Finally, identification of VSIG4 for immunotherapy response in patients with glioma indicates that the VSIG4 signature may have the ability to predict the effect of immunotherapy in glioma. Figure 1 depicts the procedures used for data preparation, processing, analysis, and validation. Gene expression profiles of GSE2223 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi). e GPL1833 (HG-U133A) Affymetrix Human Genome U133A array was used to extract the expression profile of the GSE2223 dataset. e dataset consisted of 50 glioma and four control samples. e original expression profile was downloaded, followed by background correction and quantile normalization using a robust multiarray averaging (RMA) algorithm. Subsequently, mRNA expression matrices were developed for patients with glioma.

Screening and Analysis of DEGs.
Cuffdiff was used to screen for DEGs between glioma and control samples. e threshold for DEGs was set as follows: adjusted p value <0.02; log2 (fold change) > 1.5 or log2 (fold change) < −1.5. e ggplot2 package in R software (R Foundation for Statistical Computing, Vienna, Austria) was used to display heat and volcano maps.
DEGs with unique biological significance were identified using Gene Ontology (GO) analysis. e Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to search for important pathways. e ClusterProfiler package in R was used for GO annotation and KEGG pathway analysis.

Coexpression Network Constructed Using WGCNA.
e WGCNA R package was used to construct the coexpression network of all genes in glioma and control samples. e algorithm screened the top 25% of genes for further analysis. e WGCNA analysis was performed on 50 glioma and four control samples. e samples were used to calculate the Pearson correlation matrix. Weighted adjacency matrix was established using the formula amn� |cmn|β (where amn� adjacency between genes m and n; cmn� Pearson's correlation; β� soft-power threshold). Furthermore, the weighted adjacency matrix was transformed into a topological overlap metric matrix (TOM) to estimate its connectivity in the network. e average linkage hierarchical clustering method was used to construct a clustering dendrogram of the TOM matrix. e minimum gene module size was set at 30 to obtain suitable modules, and the threshold for merging similar modules was set at 0.25.

Hub Gene Identification Using WGCNA and PPI.
e WGCNA modules that were most significantly correlated with the clinical phenotypes were identified. en, the most relevant MEyellow and MEbrown modules were selected for subsequent analysis. To construct a PPI network, 311 DEGs were uploaded to the STRING database (https://string-db. org/) [13]. e confidence score was set at 0.9. e core modules of the PPI network were analyzed using Cytoscape software and plugin (cytoHubba) [14], and the degree scores of the first 100 genes were screened. Gene overlapping between cytoHubba and the yellow-brown modules were identified as hub genes.

Validation of Hub
Genes. To validate the abnormal expression level of hub genes, the online database GEPIA [15] was used to analyze the hub gene expression profiles of glioma tissue and normal brain samples from the Cancer Genome Atlas (TCGA) database.
Results of the immunohistochemical staining were collected from the Human Protein Atlas to verify the protein levels of hub genes in glioma tissue and normal tissue (HPA, https://www.proteinatlas.org/). HPA is a Sweden-based program initiated in 2003, which maps human proteins in cells, tissues, and organs. We used the GeneMANIA online platform to analyze the hub genes and their networks of coexpressed genes (https://genemania.org/) [16].

Survival Analysis.
e Kaplan-Meier curve of GEPIA [15] was used to analyze the overall survival (OS) and disease-free survival (DFS) for hub genes in TCGA glioma patients.
e methods used were in accordance with the publisher's instructions. Only genes with p values <0.05 were considered potential prognostic genes.

Reverse-Transcription and Quantitative
Real-Time PCR (qRT-PCR). Glioma (n� 12) and normal brain (from traumatic decompression patients, n� 6) tissues were collected from the Neurosurgical Department of Tongji Hospital after obtaining written consent and approval from the Research Ethics Committee of Tongji Hospital (no. TJ-IBR20181111). Tissue total RNA was isolated using the Trizol reagent (Takara Bio Inc., Shiga Japan) according to the manufacturer's instructions. e RT Premix kit (Takara Bio Inc.) was used for reverse transcription. e reaction was conducted at 37°C for 15 min and at 85°C for 5 s. Quantitative real-time PCR (qRT-PCR) was performed using the real-time PCR system (Bio-Rad Laboratories, Hercules, CA, USA) and SYBR Green PCR Master Mix (Toyobo Co., Ltd., Osaka, Japan). Table 1 lists the primer sequences.

Analysis of Immune Cell Infiltration.
e CIBERSORT algorithm was used to evaluate the percentage of 22 immune cell types in each sample. e fraction of 22 immune cells was compared between the glioma and control groups, and the violin plot was drawn by the "vioplot" R package. e correlation coefficient between immune cells was calculated using the "corrplot" R package. Spearman correlation analysis was also performed to investigate the correlation of    Figures 1A and 1B).
To explore the relationship between DEGs, they were annotated using GO and KEGG. e GO analysis showed that the most important GO terms were "modulation of chemical synaptic transmission" (body: biological process), "synaptic vesicle membrane" (body: cell component), and "calmodulin binding" (body: molecular function) (Supplementary Figure 2A). In addition, KEGG pathway analysis revealed that the most significant enrichment pathway was "insulin secretion" (Supplementary Figure 2B).

Construction of the Coexpression Network.
Raw data were normalized using the RMA method in the limma package. Genes with a false discovery rate <0.05 and log2 fold change ≥0.5 were included in the WGCNA analysis. First, genes and samples were examined with the missing values, and all of them met the threshold. Next, the samples were clustered to identify any significant outliers. e height cut-off value was set at 30, and all samples were included in the analysis (Figure 2(a)).
To construct the WGCNA network, the soft threshold power β was calculated, and the coexpression similarity was proposed to calculate the adjacency. e pick soft-threshold function in WGCNA was used to analyze the network topology. In subsequent analyses, the soft threshold power β was set at 11 due to the scale independence of 0.9 and relatively high average connectivity (Figure 2(b)). e gene network was constructed and modules were identified using the one-step network construction function of the WGCNA R package. For cluster splitting, the soft threshold power was set at 11, minimum module size at 30, and deepSplit at 2 (which correlated with medium sensitivity). Finally, 11 gene coexpression modules were constructed (Figure 2(c)).

Construction of Coexpression Modules and Identification of Key Modules.
Relationships between the identified modules were mapped, and the connectivity of eigengenes was analyzed. Eigengenes provide pairings between gene coexpression modules.
e results showed that the 11 modules could be clustered into two clusters ( Figure 2(d)). e associations between modules and clinical characteristics were evaluated to identify the most significant association.
e results of this analysis showed that modules (yellow and brown modules) were most significantly correlated with glioma (correlation coefficients: 0.61 and −0.6, respectively). erefore, the two modules were selected for further analyses (Figure 2(e)). A scatterplot of gene significance and module membership was plotted in the yellow and brown modules (Supplementary Figure 3).

Identification of Hub Genes.
To investigate the functions of DEGs, a PPI network was constructed using the STRING database to provide a visual annotation network for identifying the structural and functional properties of proteins ( Figure 3(a)). en, cytoHubba was used to detect the key genes of the PPI network. e top 100 gene networks are shown in Figure 3(b). Using the Venn package in R, the overlapping genes were screened in the yellow-brown module and cytoHubba ( Figure 3(c)). We screened 60 overlapping genes; the differences in hub gene expression in the top nine hub genes in glioma are shown in Table 2.

Validation of Hub Gene Expression Levels.
e expression levels of the top nine genes in glioma were validated using GEPIA. DEGs between LGG and GBM were considered hub genes. We identified six hub genes: RAB3A, SYP, CAMK2A, and GABRA1 had a lower expression level, while TYROBP and VSIG4 had a higher expression level in 163 GBM and 518 LGG tissues, compared to 207 normal samples ( Figure 4). Meanwhile, the results of immunohistochemical staining showed that the expression level of corresponding proteins in glioma tissues was consistent with the transcription level of hub genes selected from the GEO database (Supplementary Figure 4). To explore the relationship between hub genes in glioma, the interaction network of hub genes and their coexpression genes were analyzed using the GeneMANIA online platform (Figure 3(d)).

Expression Level of Hub Genes in the GSE2223 Dataset.
In the GSE2223 dataset, the expression levels of six hub genes (RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1) were consistent with those from the GEPIA database (Figures 5(a)-5(f )).
e areas under the curve (AUC) for all six hub genes were ≥0.89, while the receiver operating characteristic (ROC) curve was p ≤ 0.03 (Figures 5(g)-5(l), Table 3). erefore, the six hub genes had a good diagnostic value for glioma.

Prognostic Value of Hub Genes in Glioma.
Using the GEPIA database, the survival curves were produced to explore the prognostic value of hub genes. e higher level of expression of RAB3A, SYP, CAMK2A, and GABRA1 in glioma patients was associated with improved OS and DFS (Figures 6(a), 6(a)-6(d), 6(f)-6(g), 6(i)-6(j), and 6(l)), whereas 4 Journal of Oncology     Journal of Oncology 7 the higher level of expression of TYROBP and VSIG4 was associated with reduced OS and DFS (Figures 6(b), 6(e), 6(h), and 6(k)). erefore, the six hub genes had a significant prognostic value in glioma patients and predicted the OS rate and progression-free interval (PFI).

Verification of Hub Gene Expression Levels in Glioma
Tissues.
e levels of mRNA expression of the hub genes were measured using quantitative PCR in glioma and normal tissues. Compared to the control tissues, glioma tissues had a decreased expression of RAB3A, SYP, CAMK2A, and GABRA1 but an increased expression of TYROBP and VSIG4 (Figures 7(a)-7(f )). ese results are in agreement with the findings from the GEO microarray, GEPIA database, and immunohistochemical staining in the Human Protein Atlas database. e AUC and ROC values for the six hub genes were ≥0.83 and p ≤ 0.03, respectively. erefore, all hub genes had a good diagnostic value for glioma. e data from the hub genes were consistent with those from the GSE2223 dataset (Figures 7(g)-7(l), Table 4).  e CIBERSORT algorithm was used to analyze the difference of immune cell infiltration between the glioma and control groups in 22 subpopulations of immune cells. e total value of all immune cells in each sample was set at 100%, and the proportion of each immune cell in these samples is presented in Figure 8(a). e violin plot showed marked differences in the distribution of 13 out of 22 immune cells (Figure 8(b)). Taken together, these results suggest that the heterogeneity of infiltrating immune cells in glioma is evident and may play a role in the pathogenesis of glioma.
To further investigate the correlation of RAB3A, TYROBP, SYP, CAMK2A, VSIG4, GABRA1, and infiltrating immune cells, Spearman correlation was performed and plotted in a lollipop chart (Figures 8(c)-8(h)). ese results indicate that the core gene RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1 is closely related to the level of immune cell infiltration and plays a crucial role in the immune microenvironment of glioma.

Potential of VSIG4 as an Indicator of Response to
Immunotherapy.
e TIDE score reflects the sensitivity to immune checkpoint. We evaluated the correlation between the RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1 signature and the TIDE score. We only found that the TIDE score was significantly negatively correlated with VSIG4 gene expression (Figure 8(i)). We also observed similar outcomes in VSIG4 high and low subgroups. e high expression levels of the VSIG4 group benefited more from immunotherapy (Figure 8(j)). ese findings suggested that patients of the high-expression VSIG4 group may be more sensitive to immunological treatment.

Discussion
In the present study, we analyzed 50 glioma samples and four normal samples from GSE2223 and constructed a gene coexpression network based on WGCNA. Overlapping genes were confirmed in the yellow-brown module and PPI network of DEGs using the cytoHubba plugin. We identified six hub genes (RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1) that predicted the prognosis of glioma. ese findings were supported by the RT-qPCR test, which is used in the clinical setting. Finally, the identification of VSIG4 for immunotherapy response in patients with glioma demonstrates utility for immunotherapy research.
After screening a series of databases (GEO, TCGA, HPA, STRING, GEPIA, and GeneMANIA) in GSE2223, 11 modules were identified for glioma progression, followed by six hub genes for prognosis of glioma. RAB3A regulates calcium-dependent lysosome exocytosis and plasma membrane repair through interaction with two effectors: SYTL4 and myosin-9/MYH9 [17]. RAB3A positively regulates acrosome content secretion in sperm cells by interacting with RIMS1 [18,19]. TYROBP is tyrosine-phosphorylated in the ITAM domain following ligand binding by the associated receptors, which leads to the activation of additional tyrosine kinases and subsequent cell activation [20]. TYROBP stabilizes the TREM2 C-terminal fragment produced by TREM2 ectodomain shedding, which suppresses the release of proinflammatory cytokines [21]. SYP is possibly involved in organizing membrane components and targeting vesicles to the plasma membrane, as well as short-term and longterm synaptic plasticity [22]. CAMK2A regulates dendritic     spine development [23] and migration of developing neurons [24]. VSIG4 phagocytic receptor negatively regulates T-cell proliferation and IL-2 production [25]. e RAB18-VSIG4 interaction was involved in reducing glioma proliferation and increasing apoptosis, as well as reducing TMZ sensitivity [26]. Let-7g-5p could inhibit epithelialmesenchymal transition of glioblastoma by targeting VSIG4, which was consistent with the reduction of the glioma stem cell (GSC) phenotype [27]. High expression of VSIG4 was associated with poor prognosis of OS and PFS in high-grade glioma patients [28]. e GABRA1 ligand-gated chloride channel, which is a component of the heteropentameric receptor for GABA, is the major inhibitory neurotransmitter in the brain [29][30][31]. GABRA1 plays an important role in the formation of functional inhibitory GABAergic synapses and synaptic inhibition of GABAgated ion channels [29,30]. is trend was confirmed by immunohistochemical staining of hub genes in Human Protein Atlas. e survival curves were produced using the GEPIA database to explore the prognostic value of hub genes. Higher expression levels of RAB3A, SYP, CAMK2A, and GABRA1, as well as lower expression levels of TYROBP and VSIG4, in glioma patients were associated with improved OS and DFS. is suggests that all hub genes have a significant prognostic value in glioma and can predict the OS and PFI event.

B cells naive B cells memory Plasma cells T cells CD8 T cells CD4 naive T cells CD4 memory resting T cells CD4 memory activated T cells follicular helper T cells regulatory (Tregs) T cells gamma delta NK cells resting NK cells activated Monocytes
e mRNA expression levels for the hub genes in tissue samples were similar to those obtained from the GEO microarray, GEPIA database, and Human Protein Atlas database. AUC values for the hub genes demonstrated a good diagnostic value for glioma, which suggests that they may be useful as diagnostic biomarkers and therapeutic targets in glioma.
In addition, immune infiltration analysis in this study demonstrated that the changes of infiltrating immune cells in glioma are evident. Interestingly, RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1 were also found to be closely related to the level of immune cell infiltration in the current study. erefore, it could be concluded that RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1 may play a critical role in glioma by regulating immune cells. We also observed the TIDE score was significantly negatively correlated with VSIG4 gene expression, indicating that the VSIG4 high expression group may be more sensitive to immunological treatment. Different from the classical complement receptors CR3 and CR4, the unique function of VSIG4 suggests unique functions in the regulation of innate and acquired immunity [32]. Soluble Vsig4-IG could attenuate the induction of T cell responses in vivo and inhibit cell-dependent responses [25]. VSIG4 signaling inhibited the proliferation of CD4(+) and CD8(+) T cells and the production of IL-2 and IFN-c in coculture in vitro [33]. Besides, after coculture with DCs transfected with hVSIG4 recombinant adenovirus, T cell proliferation potential, cytokine production, and activation marker expression were suppressed [34]. All the above suggests that VSIG4 may contribute to the development of new immunotherapy strategies.
ere were some limitations to our study. First, this study included a small sample size. Second, we did not have data for the survival curve analysis (OS and DFS) of glioma patients. ird, the corresponding protein levels for the hub genes were not measured. Fourth, the molecular mechanisms underlying the relationships between hub genes and glioma diagnosis and prognosis were not studied. erefore, further studies are needed to validate the new therapeutic targets.

Conclusions
After WGCNA and PPI network analysis of DEGs in glioma samples, a series of multiple databases were screened. Six hub genes related to glioma diagnosis and prognosis were identified, including RAB3A, TYROBP, SYP, CAMK2A, VSIG4, and GABRA1. ese genes may help to identify potential therapeutic targets and diagnostic and prognosis biomarkers. We identified VSIG4 for immunotherapy response in patients with glioma, demonstrating utility for immunotherapy research.

Data Availability
Publicly available datasets were analyzed in this study. ese data can be found in the following link: https://www.ncbi. nlm.nih.gov/sra/?term�SRP049695.

Conflicts of Interest
e authors declare no conflicts of interest.