Comprehensive Analysis of CD163 as a Prognostic Biomarker and Associated with Immune Infiltration in Glioblastoma Multiforme

Background Glioblastoma multiforme (GBM) is the most common and aggressive primary malignancy in adults with high aggression. The prognosis of GBM patients is poor. There is a critical need for novel biomarkers for the prognosis and therapy of GBM. Methods Differentially expressed genes (DEGs) in GBM were screened using TCGA cohort. Univariate and multivariate Cox regression analyses were performed on DEGs to identify the optimal prognosis-related genes. qRT-PCR was performed to verify the result. Results A total of 5216 DEGs, including 2785 upregulated and 2458 downregulated genes, were obtained. Enrichment analysis revealed that these DEGs were mainly involved in the p53 signaling pathway and cell cycle, immune response, and MAPK signaling pathways. Moreover, the top 50 DEGs were associated with drug resistance or drug sensitivity. Prognosis analysis revealed that GBM patients with a high expression of CD163 and CHI3L2 had a poor overall survival, prognosis-free survival, and disease-specific survival. The univariate and multivariate analyses revealed that CD163 and age were independent factors affecting the prognosis of GBM patients. A validation study revealed that CD163 was upregulated in GBM tissues and associated with poor overall survival. Moreover, further analysis revealed that CD163 showed significant correlation with immune cells, immune biomarkers, chemokines, and chemokine receptors. We also identified several CD163-associated kinase, miRNA, and transcription factor targets in GBM, including LCK, miR-483, and ELF1. Conclusions In conclusion, our study suggested CD163 as a prognostic biomarker and associated it with immune infiltration in GBM.


Introduction
Glioblastoma multiforme (GBM) is the most common and aggressive primary malignancy in adults with high aggression [1]. Until now, the etiology and pathogenesis of GBM are still far from clarified [2]. The standard treatment for GBM includes surgical tumor removal followed by ionizing radiation and alkylating chemotherapy [3]. However, there is no prognosis in the standard treatment for glioblastoma in the past two decades [4,5]. The prognosis of glioblastoma patients is poor, with a median survival timeline of about 12 months, with a 5-year survival of about 10% [6]. These sober-ing data illustrate a critical need for novel biomarkers for the prognosis and therapy of GBM.
The immune microenvironment has been chronicled to exert a significant function in biological progress in cancer [7]. Immunotherapy based on immune checkpoint blockade is ever-increasingly suggested as the most promising therapy for GBM in addition to operative treatment, especially for patients with advanced GBM [8]. Though many immunotherapy methods, such as GBM vaccines, oncolytic viral therapies, immune-checkpoint suppressors, and chimeric antigen receptor T cell therapy, have been conducted in clinical trials, none of these have been applied for clinical treatment [7,8]. Similarly, though some prognosis biomarkers, including MLK3 and P4HA1, have been identified in GBM at the genetic level, little of these have been applied for the prediction of the prognosis of patients [9][10][11]. Thus, it is necessary to identify novel biomarkers for the prognosis and therapy of GBM.
In recent years, with the development of sequencing technology and the establishment of various cancer databases, genomic research has become one of the most reliable means to accelerate the clinical and translational research and treatment of cancer. In our study, we aim to identify the prognosis biomarkers and therapy targets for GBM by mining databases. Our result may provide more suitable strategies to improve the anti-immune performance and prognosis prediction of GBM by using a high-throughput sequencing database.

Database and Gene Expression. Level 3 gene expression profiles (level 3 data) for GBM patients were isolated from
The Cancer Genome Atlas (TCGA) data portal (https:// tcga-data.nci.nih.gov/tcga/) and Chinese Glioma Genome Atlas (CGGA) data portal (http://www.cgga.org.cn/). The limma package (version: 3.40.2) of R software was used to explore the differential expression genes. The adjusted p value was analyzed to correct for false positive results in TCGA. "Adjusted p < 0:05 and log ðfold changeÞ > 2 or log ðfold changeÞ < −2" were defined as the thresholds for the screening of differential expression genes (DEGs).

Heat Maps and Volcano Plots.
Heat maps and volcano plots about DEGs were obtained using an R Project.

Enrichment Analysis.
To further confirm the underlying function of potential targets, the data were analyzed by functional enrichment. Gene Ontology (GO) is a widely used tool for annotating genes with functions, especially molecular function (MF), biological pathways (BP), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis is a practical resource for analytical study of gene functions and associated high-level genome functional information. To better understand the carcinogenesis of mRNA, the clusterProfiler package in R was employed to analyze the GO function of potential targets and enrich the KEGG pathway.

Drug Sensitivity Analysis.
To analyze the correlation of DEGs and drug sensitivity, we collected 265 small molecules from the Genomics of Drug Sensitivity in Cancer (GDSC). We downloaded the area under the dose-response curve (AUC) values for drugs and gene expression profiles for all cancer cell lines. The Pearson correlation analysis was utilized to explore the correlation between DEGs and small molecules or drugs.

Survival Analysis.
After separating the high/low expression group of GBM with the medium expression of DEGs, we used the Kaplan-Meier survival analysis to analyze the survival difference between these two groups. Log-rank tests were used to calculate p values and hazard ratio (HR) with 95% confidence interval (CI). Univariate and multivariate cox regression analyses were conducted to detect the proper terms to build the nomogram. The forest was used to show the p value, HR, and 95% CI of each variable through the "forestplot" R package. A nomogram was developed based on the results of a multivariate Cox proportional hazard analysis to predict the 1-year, 2-year, and 3-year overall survival. The nomogram provided a graphical representation of the factors, which can be used to calculate the risk of survival for an individual patient by the points associated with each risk factor through the "rms" R package. p < 0:05 was considered as statistically significant.
2.5. Genetic Mutation Landscape. The genetic mutation data was obtained from the TCGA database. The "maftools" package in R software was applied to analyze and visualize the genetic mutation landscape. A horizontal histogram showed that the genes have the higher mutation frequency in GBM patients.
2.7. LinkedOmics. LinkedOmics (http://www.linkedomics .org/) is a bioinformatics platform for genomic analysis based on the TCGA dataset [15]. The "interpreter module" of Lin-kedOmics performs pathway and network analyses of CD163. Data from the LinkFinder results were signed and ranked, and GSEA was used to perform analyses of GO (CC, BP, and MF), KEGG pathways, kinase-target enrichment, miRNA-target enrichment, and transcription factortarget enrichment. The rank criterion was a p value < 0.05, and 500 simulations were performed.
2.8. GeneMANIA. GeneMANIA (http://www.genemania.org/) is a bioinformatics tool developed for protein-protein interaction (PPI) network analysis and for promoting understanding of the functional association data of target genes [16]. In order to better understand the function behind the CD163associated kinase_LCK network, miR-483 network, and transcription factor ELF1 network, we submitted these genes to GeneMANIA to construct a PPI network.
2.9. Validation of the Expression and Prognosis Value. The immunohistochemistry staining of target genes in GBM tissues and normal tissues was obtained from The Human Protein Atlas (https://www.proteinatlas.org/), a bioinformatics tool aimed at mapping all the human proteins in cells, 2 BioMed Research International     BioMed Research International tissues, and organs using an integration of various omics technologies [17].
GBM and normal brain tissues (n = 52) were obtained from patients from the Shengjing Hospital of China Medical University. All patients provided informed consent. Each patient did not receive any treatment before operation. Total RNA of human tissues was extracted with a TRIzol reagent (Vazyme, Nanjing, China). The synthesis of cDNAs corresponding to the mRNAs of interest depended on PrimeScript RTpolymerase (Vazyme) and SYBR-Green Premix (Vazyme) with specific PCR primers (Sangon Biotech Co., Ltd., Shanghai, China). Glyceraldehyde-3-phosphate dehydrogenase was used as an internal control. The 2 −ΔΔCt method was used to calculate fold changes. Primer sequences were as follows: GAPDH, forward: GCACCGTCAAGGCTGAGAAC, reverse: TGGTGA AGACGCCAGTGGA, and CD163, forward: CTACGAACT TCAGCCAGAGTGCACCTCAT, reverse: GTCATAATGAA CTTCAGCTATTGCACAC. The differences of the CD163 expression and the prognosis of CD163 in GBM were evaluated with Student's t-test and Kaplan-Meier analysis in Graph-Pad Prism 7 software (GraphPad, Inc., La Jolla, CA, USA).

Identification of DEGs and Associated Functions in GBM.
The DEGs between GBM tissues and brain tissues were explored using TCGA GBM dataset. As a result, we obtained 5216 DEGs including 2785 upregulated and 2458 downregulated genes (Figure 1(a)). Figure 1(b) shows the top 50 upregulated and downregulated genes. In order to explore the potential functions of DEGs in GBM, we performed enrichment analysis, including GO analysis and KEGG pathway analyses. As shown in Figure 1(c), these upregulated genes were mainly involved in the p53 signaling pathway and  Figure 3: Drug sensitivity analysis of differentially expressed genes in GBM. The positive correlation means that the gene high expression is resistant to the drug and vice versa. 6 BioMed Research International ribosome, proteoglycans in cancer, focal adhesion, DNA replication, and cell cycle in KEGG pathways. Moreover, GO analysis revealed that upregulated genes were mainly involved in vital transcription, vital gene expression, translational initiation, neutrophil activation involved in immune response, T cell activation, and DNA replication (Figure 1(c)). As for the result of downregulated genes, KEGG pathway analysis suggested that these genes were mainly involved in cAMP signaling pathways, the synaptic vesicle cycle, oxytocin signaling pathways, MAPK signaling pathways, and calcium signaling pathways (Figure 1(d)). Furthermore, GO analysis suggested that these downregulated genes were mainly involved in vesiclemediated transport in synapse, synaptic vesicle exocytosis, signal release from synapse, regulation of neurotransmitter levels, potassium ion transport, cognition, and axonogenesis ( Figure 1(d)).

Somatic Mutations in the GBM.
To identify the somatic mutations of the patients with GBM in the TCGA database, we downloaded mutation data from TCGA and visualized using the "maftools" package in R software. A horizontal histogram showed that the genes have a higher mutation frequency in GBM patients, including TTN (25%), TP53 (30%), PTEN (30%), EGFR (24%), and MUC16 (14%) (Figures 2(a) and 2(b)). Missense mutations and nonsense mutations were the two most common types of mutation in GBM patients (Figures 2(a) and 2(b)). Scanning the variant types of mutations in GBM, single nucleotide polymorphism (SNP) was the most common type (Figure 2(b)). Moreover, C>T was the predominant mutation type in GBM (Figure 2(b)).

Drug Sensitivity Analysis of Top 50 DEGs in GBM.
To develop cancer pharmacotherapy, a crucial way is to assess the link between DEGs and existing drug targets. In our study, we selected the top 50 DEGs (Supplementary Table 1) for further analysis and performed a drug sensitivity analysis. To explore the correlation of DEGs and drug sensitivity, the Pearson correlation coefficients of transcript levels and AUCs were used and normalized based on Fisher's Z transformation based on 265 small molecules from the Genomics of Drug Sensitivity in Cancer (GDSC), which was used before [18]. We observed that most of the top 50 DEGs show drug resistance (positive correlation) or drug sensitivity (negative correlation) ( Figure 3).

Prognosis Value of Top 50 DEGs in GBM.
We then performed prognosis value of top 50 DEGs in GBM, and the genes that were statistically significant in the overall survival, prognosis-free survival, and disease-specific survival are shown in Table 1. In overall survival, a total of 14 genes were   9 BioMed Research International associated with the prognosis of GBM patients. In prognosis-free survival, a total of 13 genes were associated with the prognosis of GBM patients. In disease-specific survival, a total of 11 genes were associated with the prognosis of GBM patients. Interestingly, the result revealed that GBM patients with a high expression of CD163 and CHI3L2 had a worse overall survival, prognosis-free survival, and disease-specific survival (Table 1, all p < 0:05).

Building a Predictive Nomogram.
We then resorted to a nomogram to construct a predictive model, considering clinicopathologic features and potential prognostic biomarkers,

12
BioMed Research International to construct a clinically applicable method that could predict the survival probability of the GBM patient. The univariate and multivariate analyses revealed that CD163 and age were independent factors affecting the prognosis of GBM patients (Figures 6(c) and 6(d)). We generated a nomogram to predict the 1-year, 2-year, and 3-year OS rates in the discovery group using the Cox regression algorithm (Figure 6(e)). The calibration plots for the 1-year and 3-year OS rates were predicted relatively well compared with an ideal model in the entire cohort ( Figure 6(f)).

Validation of the Expression and Prognostic
Value of CD163 in GBM. CD163 was selected for further study, and we performed validation research then. The immunohistochemistry staining from The Human Protein Atlas revealed that the immunohistochemistry staining of CD163 in GBM tissues and normal tissues was medium and not detected (Supplementary Figure 1A). Due to a similar role of CD8 and CD163 in immune infiltration, we also detected CD8 expression in GBM. The immunohistochemistry staining of CD8 in GBM tissues and normal tissues was medium and not detected (Supplementary Figure 1B). Moreover, qRT-PCR was performed to verify the expression and prognostic value of CD163 in GBM. As expected, CD163 expression was increased in GBM tissues (Supplementary Figure 2A Figure 2C). These data further verified our result obtained above.

CD163-Associated Functions in GBM.
In order to clarify the CD163-associated functions in GBM, we performed enrichment analysis using GSEA. The items in GO analysis are shown in Figure 8(a), revealing that CD163 were mainly involved in adaptive neutrophil-mediated immunity, immune response, leukocyte cell-cell adhesion, cytokine receptor binding, cytokine binding, and immunoglobulin binding. Furthermore, CD163 were mainly involved in cytokine-cytokine receptor interaction, NOD-like receptor signaling pathway, NF-kappa B signaling pathway, and TNF signaling pathway in KEGG pathway analysis (Figure 8(b)).
3.9. CD163-Associated Kinase, miRNA, or Transcription Factor Targets in GBM. In order to further clarify the underlining mechanisms about how CD163 affected the tumorigenesis  and progression of GBM, we finally explore CD163-associated kinase, miRNA, or transcription factor targets in GBM using GSEA in LinkedOmics. As a result, the result indicated that the top 5 most significant CD163-associated kinase targets in GBM were LCK, LYN, SYK, HCK, and ATR (Table 3, All p < 0:05). The PPI network based on the correlated genes of kinase LCK constructed by GeneMANIA indicated that kinase LCK was mainly related to the antigen receptor-mediated signaling pathway, immune response, and T cell receptor signaling pathway and activation ( Figure 9). Moreover, the top 5 most significant CD163-associated miRNA targets in GBM were miR-483 (AGGAGTG), miR-485-5P (CAGC CTC), miR-197 (GTGGTGA), miR-499 (AGTCTTA), and miR-331 (CCAGGGG) ( Table 3, All p < 0:05). The PPI network based on the correlated genes of miR-483 constructed by GeneMANIA indicated that miR-483 were mainly related to the regulation of lymphocyte activation, regulation of cell activation, and immune response (Supplementary Figure 4). The top 5 most significant CD163-associated transcription factor targets in GBM were V$ELF1_Q6, V$PEA3_Q6, V$E2F1_Q6, V$BACH2_01, and V$TEL2_Q6 (Table 3, All p < 0:05). The PPI network based on the correlated genes of ELF1 constructed by GeneMANIA indicated that ELF1 were mainly related to the regulation of transcription initiation from RNA polymerase II promoter, mediator complex, and nuclear hormone receptor binding (Supplementary Figure 5).

Discussion
The oncogenesis and progression of GBM is a complex multistep process, involved in a variety of gene expression patterns and other factors. Considering the heterogeneity and complex mechanism of GBM, clarifying the molecular mechanism of GBM is of significant importance for the ther-apy of GBM patients [23]. Moreover, the prognosis of GBM patients was poor. Though multidisciplinary comprehensive treatment, including surgery and chemo-and radiation therapy, had been used for GBM patients, the median survival time of GBM patients is only about 15 months [24]. And the five-year survival rate of GBM is about 0.05% to 4.7% [25]. Thus, it is quite necessary to explore new therapeutic targets and prognostic markers of GBM. In order to explore new therapeutic targets and prognostic markers of GBM, we first identify the DEGs by comparing GBM tissues with normal tissues in the TCGA cohort. As a result, a total of 5216 DEGs including 2785 upregulated and 2458 downregulated genes were obtained. We then selected the top 50 DEGs for further analysis. In order to explore the potential of the 50 DEGs as the therapy targets of GBM patients, we detect the relation between DEGs and existing drug targets. Interestingly, we observed that most of the top 50 DEGs show drug resistance (positive correlation) or drug sensitivity (negative correlation). Therefore, these 50 DEGs had potential as the therapy targets of GBM patients, and further study should be performed to verify this result.
We explore the potential of the 50 DEGs as the prognostic biomarkers of GBM patients by performing prognosis analysis. And the data indicated that CD163 and CHI3L2 may serve as prognostic biomarkers in GBM and GBM patients with a high expression of CD163 and CHI3L2 which predicted a poor overall survival, prognosis-free survival, and disease-specific survival. These were consistent with a previous result, which found that CD163 predicts poor prognosis in glioma patients [26]. Actually, these CD163 and CHI3L2 have been suggested as prognostic biomarkers in other types of cancers. In oral squamous cell carcinoma, CD163 was a prognostic biomarker and associated with poor survival

20
BioMed Research International [27]. Another study suggested that high CD163 expression indicated a poor prognosis of patients with urothelial cell carcinoma [28]. In breast cancer, CD163 was related to postoperative radiotherapy and poor prognosis, indicating CD163 as a prognostic marker in breast cancer [29]. CHI3L2 was suggested as a prognostic biomarker for renal cell carcinoma, predicting high risk for postoperative progression [30]. Moreover, univariate and multivariate analyses were performed and demonstrated that CD163 and age were independent factors affecting the prognosis of GBM patients. And we select CD163 for further analysis.
Another important finding of our study is that CD163 was positively correlated with immune cells, immune bio-markers, chemokines, and chemokine receptors. We found that CD163 expression was associated with the abundance of CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells. Moreover, CD163 expression was positively correlated with most gene markers of these tumorinfiltrating immune cells in both the TCGA and CGGA cohorts, including CD8A, CD8B, STAT6, STAT5A, and PDCD1. Actually, increasing evidences revealed that these immune cells and immune biomarkers exerted vital functions in tumor immune infiltration or served as a therapy target in GBM. The CD4+ T cell was linked to tumor angiogenesis and tumor progression in glioma patients [31]. Another study suggested that neutrophil-induced ferroptosis   [31]. Moreover, monocytes could serve as a promising predictor for therapy responses of glioma patients [32]. Hung et al. suggested that PDCD1 and TIGIT dual checkpoint blockade enhances antitumor immunity and survival in GBM [33]. Moreover, STAT5A was a prognosis marker for GBM and involved in immune infiltration in GBM [34]. These evidences indicated that the possible association between CD163 and immune infiltrates in GBM and CD163 may serve as an immunotherapy target of GBM patients.
There is no doubt that our study had some limitations. Firstly, most analysis was performed at the mRNA level but not the protein level, and double immunohistochemistry staining should be performed to verify the protein expression of CD163 in GBM. Furthermore, it would be better to validate our results by performing in vivo and in vitro experiments.
In conclusion, our study suggested CD163 as a prognostic biomarker and associated it with immune infiltration in GBM.

Data Availability
The analyzed datasets generated during the study are available from the corresponding author on reasonable request.