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

Abstract

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. 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 sobering 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–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.

2. Materials and Methods

2.1. 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 \text{fold}\text{change}>2$ or $\log \text{fold}\text{change}<-2$” were defined as the thresholds for the screening of differential expression genes (DEGs).

2.1.1. Heat Maps and Volcano Plots

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

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

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

2.4. 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.6. Immune Infiltration Analysis

Spearman correlation analysis was performed to explore the relation between gene expression and immune cell infiltration and the expression of immune biomarkers as well as immune checkpoints in TIMER (https://cistrome.shinyapps.io/timer/) [12] and CIBERSORT (https://cibersortx.stanford.edu/) [13]. These gene markers of tumor-infiltrating immune cells included markers of CD8+ T cells, T cells (general), B cells, monocytes, TAMs, M1 macrophages, M2 macrophages, neutrophils, natural killer (NK) cells, dendritic cells (DCs), T-helper 1 (Th1) cells, T-helper 2 (Th2) cells, follicular helper T (Tfh) cells, T-helper 17 (Th17) cells, Tregs, and exhausted T cells [14]. A $p$ value of less than 0.05 was considered statistically significant.

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 LinkedOmics 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 factor-target 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 CD163-associated 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, 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 RT-polymerase (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: TGGTGAAGACGCCAGTGGA, and CD163, forward: CTACGAACTTCAGCCAGAGTGCACCTCAT, reverse: GTCATAATGAACTTCAGCTATTGCACAC. 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 GraphPad Prism 7 software (GraphPad, Inc., La Jolla, CA, USA).

3. Results

3.1. 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 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 vesicle-mediated transport in synapse, synaptic vesicle exocytosis, signal release from synapse, regulation of neurotransmitter levels, potassium ion transport, cognition, and axonogenesis (Figure 1(d)).

[figures omitted; refer to PDF]

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

[figures omitted; refer to PDF]

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

[figure omitted; refer to PDF]

3.4. 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 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$). To be more specific, GBM patients with a high expression of CD163 had a poor overall survival (Figure 4(a), $p=0.0437$, $\text{HR}95\%\text{CI}=1.451.01-2.09$), prognosis-free survival (Figure 4(b), $p=0.012$, $\text{HR}95\%\text{CI}=1.601.11-2.30$), and disease-specific survival (Figure 4(c), $p=0.048$, $\text{HR}95\%\text{CI}=1.481.00-2.17$) with a 3-year AUC of 0.744 (Figure 4(a)), 0.63 (Figure 4(b)), and 0.749 (Figure 4(c)), respectively. And the risk score of each patient is shown in Figure 4. Moreover, GBM patients with a high expression of CHI3L2 had a poor overall survival (Figure 5(a), $p=0.0077$, $\text{HR}95\%\text{CI}=1.611.12-2.33$), prognosis-free survival (Figure 5(b), $p=0.0063$, $\text{HR}95\%\text{CI}=1.641.14-2.37$), and disease-specific survival (Figure 5(c), $p=0.022$, $\text{HR}95\%\text{CI}=1.581.07-2.35$) with a 3-year AUC of 0.664 (Figure 5(a)), 0.678 (Figure 5(b)), and 0.698 (Figure 5(c)), respectively. And the risk score of each patient is shown in Figure 5. In order to further verify our result, we then submitted CD163 and CHI3L2 to the CGGA cohort and performed a prognosis analysis. As expected, GBM patients with a high expression of CD163 (Figure 6(a), $p<0.0001$) and CHI3L2 (Figure 6(b), $p<0.0001$) had a poor prognosis in the CGGA cohort. These data demonstrated that CD163 and CHI3L2 might serve as prognostic biomarkers in GBM.

Table 1

Prognosis analysis top 50 differentially expressed genes in GBM.

Genes	Overall survival		Progression-free survival		Disease-specific survival
	$p$ value	HR (95% CI)	$p$ value	HR (95% CI)	$p$ value	HR (95% CI)
CD163	0.0437	1.45 (1.01-2.09)	0.012	1.60 (1.11-2.30)	0.048	1.48 (1.00-2.17)
CHI3L2	0.0077	1.61 (1.12-2.33)	0.0063	1.64 (1.14-2.37)	0.022	1.58 (1.07-2.35)
TOP2A	—	—	0.006	0.60 (0.42-0.86))	—	—
MMP9	—	—	0.044	1.45 (1.01-2.09)	—	—
NUSAP1	—	—	0.029	0.67 (0.47-0.96)	—	—
PIMREG	—	—	0.007	0.60 (0.42-0.87)	—	—
F2R	—	—	0.038	0.68 (0.47-0.98)	—	—
AURKB	—	—	0.026	0.66 (0.46-0.95)	—	—
PDPN	—	—	0.002	1.79 (1.24-2.59)	—	—
ANXA1	—	—	0.022	1.53 (1.06-2.21)	—	—
KIFC1	—	—	0.038	0.68 (0.48-0.98)	—	—
PI3	—	—	0.033	1.48 (1.03-2.13)	—	—
SAA1	0.033	1.49 (1.03-2.15)	0.042	1.46 (1.01-2.11)	—	—
TGFBI	0.019	1.55 (1.07-2.24)	—	—	0.016	1.62 (1.09-2.40)
PACSIN1	0.027	1.51 (1.05-2.18)	—	—	0.035	1.52 (1.03-2.26)
SULT4A1	0.023	1.52 (1.06-2.19)	—	—	0.025	1.56 (1.06-2.30)
VSNL1	0.012	1.59 (1.11-2.27)	—	—	0.039	1.50 (1.02-2.20)
SV2B	0.006	1.66 (1.15-2.39)	—	—	0.018	1.60 (1.08-2.37)
CHGA	0.016	1.57 (1.09-2.27)	—	—	0.033	1.54 (1.04-2.28)
SNCB	0.010	1.61 (1.12-2.32)	—	—	0.039	1.51 (1.02-2.22)
SLC12A5	0.005	1.70 (1.18-2.46)	—	—	0.005	1.76 (1.19-2.62)
NRGN	0.035	1.48 (1.03-2.14)	—	—	0.099	1.39 (0.94-2.06)
CPLX2	0.039	1.48 (1.02-2.14)	—	—	—	—
CREG2	0.031	1.49 (1.04-2.15)	—	—	—	—

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

3.5. Building a Predictive Nomogram

We then resorted to a nomogram to construct a predictive model, considering clinicopathologic features and potential prognostic biomarkers, 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)).

3.6. 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, $p=1.2\ast {10}^{-6}$). Further analysis suggested and GBM patients with a high CD163 level had a poor overall survival (Supplementary Figure 2B, $p=0.026$) with an AUC of 0.723 in the ROC curve (Supplementary Figure 2C). These data further verified our result obtained above.

3.7. CD163 Were Associated with Immune Infiltration in GBM

Increasing evidences revealed that immune infiltration is an independent predictor of sentinel lymph node status and survival in cancers [14, 19, 20]. In order to explore the role of CD163 in GBM, we then detect the association between CD163 expression and immune infiltration in GBM. The result suggested that the CD163 expression was associated with the abundance of CD8+ T cells ($\text{Cor}=-0.252$, $p=1.80E-7$), CD4+ T cells ($\text{Cor}=0.132$, $p=6.89E-3$), macrophage ($\text{Cor}=0.174$, $p=3.39E-4$), neutrophils ($\text{Cor}=0.148$, $p=2.45E-3$), and dendritic cells ($\text{Cor}=0.498$, $p=1.42E-27$) (Figure 7(a)). We then verify this result using the CIBERSORT dataset, which revealed that the CD163 expression was associated with the abundance of neutrophils ($\text{Cor}=0.62$, $p=1.38E-17$), macrophages ($\text{Cor}=0.64$, $p=1.76E-18$), dendritic cells ($\text{Cor}=0.54$, $p=7.41E-54$), and CD4+ T cells ($\text{Cor}=0.22$, $p=0.006$) (Supplementary Figure 3A).

[figures omitted; refer to PDF]

Moreover, we also explore the association between CD163 expression and gene markers of tumor-infiltrating immune cells. As expected, the CD163 expression was positively correlated with most of gene markers of these tumor-infiltrating immune cells in both the TCGA and CGGA cohorts, including CD8A, CD8B, CD3D, CD3E, CD2, CD79A, CD86, CSF1R, CCL2, CD68, IL10, IRF5, COX2, VSIG4, MS4A4A, ITGAM, CCR7, KIR2DL4, HLA-DPB1, HLA-DQB1, HLA-DRA, HLA-DPA1, CD1C, NRP1, ITGAX, GATA3, STAT6, STAT5A, BCL6, STAT3, FOXP3, PDCD1, CTLA4, HAVCR2, and GZMB (Table 2). Immune checkpoints also play a vital role in immune infiltration of cancer [21]. In the current study, we found that CD163 expression increased as the expression of SIGLEC15, TIGFT, CD247, HAVCR2, PDCD1, CTLA4, and PDCD1LG2 increased (Supplementary Figure 3B, all $p<0.05$). Chemokines and their receptors modulate immune surveillance in oncogenesis, metastasis, and response to immunotherapy [22]. Interestingly, the result demonstrated a strong correlation between CD163 and chemokines as well as chemokine receptors (Figures 7(b) and 7(c)). These evidences indicated the possible association between CD163 and immune infiltrates in GBM.

Table 2

Correlation analysis between CD163 and gene biomarkers of immune cells in GBM.

Immune cells	Gene markers	GBM
		TCGA		CGGA
		Cor	$p$ value	Cor	$p$ value
CD8+ T cell	CD8A	0.202	$$\ast$$	0.441	$$\ast \ast \ast$$
	CD8B	0.241	$$\ast \ast$$	0.516	$$\ast \ast \ast$$
T cell (general)	CD3D	0.378	$$\ast \ast \ast$$	0.55	$$\ast \ast \ast$$
	CD3E	0.392	$$\ast \ast \ast$$	0.628	$$\ast \ast \ast$$
	CD2	0.422	$$\ast \ast \ast$$	0.624	$$\ast \ast \ast$$
B cell	CD19	0.097	0.231	0.321	$$\ast \ast \ast$$
	CD79A	0.187	$$\ast$$	0.32	$$\ast \ast \ast$$
Monocyte	CD86	0.58	$$\ast \ast \ast$$	0.673	$$\ast \ast \ast$$
	CD115 (CSF1R)	0.617	$$\ast \ast \ast$$	0.462	$$\ast \ast \ast$$
TAM	CCL2	0.567	$$\ast \ast \ast$$	0.652	$$\ast \ast \ast$$
	CD68	0.677	$$\ast \ast \ast$$	0.746	$$\ast \ast \ast$$
	IL10	0.593	$$\ast \ast \ast$$	0.662	$$\ast \ast \ast$$
M1 macrophage	INOS (NOS2)	-0.037	0.649	0.206	$$\ast \ast$$
	IRF5	0.317	$$\ast \ast \ast$$	0.498	$$\ast \ast \ast$$
	COX2 (PTGS2)	0.558	$$\ast \ast \ast$$	0.447	$$\ast \ast \ast$$
M2 macrophage	VSIG4	0.798	$$\ast \ast \ast$$	0.735	$$\ast \ast \ast$$
	MS4A4A	0.835	$$\ast \ast \ast$$	0.851	$$\ast \ast \ast$$
Neutrophils	CD66b (CEACAM8)	0.002	0.977	—	—
	CD11b (ITGAM)	0.607	$$\ast \ast \ast$$	0.613	$$\ast \ast \ast$$
	CCR7	0.475	$$\ast \ast \ast$$	0.497	$$\ast \ast \ast$$
Natural killer cell	KIR2DL1	0.13	0.11	—	—
	KIR2DL3	0.009	0.911	—	—
	KIR2DL4	0.236	$$\ast \ast$$	0.346	$$\ast \ast \ast$$
	KIR3DL1	0.021	0.799	—	—
	KIR3DL2	0.002	0.985	—	—
	KIR3DL3	0.105	0.195	—	—
	KIR2DS4	0.184	0.0228	—	—
Dendritic cell	HLA-DPB1	0.57	$$\ast \ast \ast$$	0.666	$$\ast \ast \ast$$
	HLA-DQB1	0.328	$$\ast \ast \ast$$	0.424	$$\ast \ast \ast$$
	HLA-DRA	0.573	$$\ast \ast \ast$$	0.736	$$\ast \ast \ast$$
	HLA-DPA1	0.467	$$\ast \ast \ast$$	0.687	$$\ast \ast \ast$$
	BDCA-1 (CD1C)	0.349	$$\ast \ast \ast$$	0.413	$$\ast \ast \ast$$
	BDCA-4 (NRP1)	0.585	$$\ast \ast \ast$$	0.763	$$\ast \ast \ast$$
	CD11c (ITGAX)	0.171	$$\ast$$	0.385	$$\ast \ast \ast$$
Th1	T-bet (TBX21)	0.031	0.707	0.233	$$\ast \ast \ast$$
	STAT4	0.318	$$\ast \ast \ast$$	0.017	0.794
	STAT1	-0.022	0.783	0.45	$$\ast \ast \ast$$
	IFN-g (IFNG)	0.134	0.0977	—	—
	TNF-a (TNF)	0.135	0.0967	0.126	0.0595
Th2	GATA3	0.277	$$\ast \ast \ast$$	0.273	$$\ast \ast \ast$$
	STAT6	0.442	$$\ast \ast \ast$$	0.651	$$\ast \ast \ast$$
	STAT5A	0.399	$$\ast \ast \ast$$	0.623	$$\ast \ast \ast$$
	IL13	-0.169	$$\ast$$	0.09	0.18
Tfh	BCL6	0.218	$$\ast \ast$$	0.295	$$\ast \ast \ast$$
	IL21	-0.009	0.91	—	—
Th17	STAT3	0.248	$$\ast \ast$$	0.681	$$\ast \ast \ast$$
	IL17A	-0.019	0.82	—	—
Treg	FOXP3	0.296	$$\ast \ast \ast$$	0.263	$$\ast \ast \ast$$
	CCR8	0.393	$$\ast \ast \ast$$	—	—
	STAT5B	-0.103	0.204	-0.079	0.238
	TGFb (TGFB1)	0.427	$$\ast \ast \ast$$	0.62	$$\ast \ast \ast$$
T cell exhaustion	PD-1 (PDCD1)	0.354	$$\ast \ast \ast$$	0.516	$$\ast \ast \ast$$
	CTLA4	0.429	$$\ast \ast \ast$$	0.362	$$\ast \ast \ast$$
	LAG3	0.023	0.777	0.37	$$\ast \ast \ast$$
	TIM-3 (HAVCR2)	0.442	$$\ast \ast \ast$$	0.622	$$\ast \ast \ast$$
	GZMB	0.332	$$\ast \ast \ast$$	0.565	$$\ast \ast \ast$$

${}^{\ast }p<0.05$, ${}^{\ast \ast }p<0.01$, ${}^{\ast \ast \ast }p<0.001$. —: no data.

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

[figures omitted; refer to PDF]

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 (CAGCCTC), 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).

Table 3

The kinase and transcription factor-target networks of CD163 in GBM.

Enriched category	Gene set	LeadingEdgeNum	$p$ value
Kinase target	Kinase_LCK	26	0
	Kinase_LYN	23	0
	Kinase_SYK	18	0
	Kinase_HCK	9	0
	Kinase_ATR	34	0
miRNA target	AGGAGTG, miR-483	24	0
	CAGCCTC, miR-485-5P	54	0
	GTGGTGA, miR-197	27	0
	AGTCTTA, miR-499	17	0
	CCAGGGG, miR-331	38	0
Transcription factor target	V$ELF1_Q6	20	0
	V$PEA3_Q6	10	0
	V$E2F1_Q6	98	0
	V$BACH2_01	73	0
	V$TEL2_Q6	62	0

[figure omitted; refer to PDF]

4. 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 therapy 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 [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 biomarkers, 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 tumor-infiltrating 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 promotes tumor necrosis in the progression of GBM [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.

Authors’ Contributions

Hao Li and Zhen Li were responsible for the design of the study and the writing the manuscript. Di Wang, Bolong Yi, Heng Cai, Zhuo Xi, and Xin Lou were responsible for the data analysis work. All authors read and approved the final manuscript.

References

[1] H. G. Wirsching, E. Galanis, M. Weller, "Glioblastoma," Handbook of Clinical Neurology, vol. 134, pp. 381-397, 2016.

[2] J. Yang, Q. Yang, "Identification of core genes and screening of potential targets in glioblastoma multiforme by integrated bioinformatic analysis," Frontiers in Oncology, vol. 10, 2020.

[3] D. J. Voce, G. M. Bernal, K. E. Cahill, L. Wu, N. Mansour, C. D. Crawley, P. A. Campbell, A. Arina, R. R. Weichselbaum, B. Yamini, "CDK1 is up-regulated by temozolomide in an NF- κ B dependent manner in glioblastoma," Scientific Reports, vol. 11, 2021.

[4] R. Stupp, M. E. Hegi, W. P. Mason, M. J. Van Den Bent, M. J. Taphoorn, R. C. Janzer, S. K. Ludwin, A. Allgeier, B. Fisher, K. Belanger, P. Hau, "Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial," The Lancet Oncology, vol. 10 no. 5, pp. 459-466, DOI: 10.1016/S1470-2045(09)70025-7, 2009.

[5] B. Campos, L. R. Olsen, T. Urup, H. S. Poulsen, "A comprehensive profile of recurrent glioblastoma," Oncogene, vol. 35, pp. 5819-5825, 2016.

[6] A. O. Sasmita, Y. P. Wong, A. P. K. Ling, "Biomarkers and therapeutic advances in glioblastoma multiforme," Asia-Pacific Journal of Clinical Oncology, vol. 14, pp. 40-51, 2018.

[7] L. Hu, Z. Han, X. Cheng, S. Wang, Y. Feng, Z. Lin, "Expression profile analysis identifies a novel seven immune-related gene signature to improve prognosis prediction of glioblastoma," Frontiers in Genetics, vol. 12, 2021.

[8] R. Batash, N. Asna, P. Schaffer, N. Francis, M. Schaffer, "Glioblastoma multiforme, diagnosis and treatment; Recent Literature Review," Current medicinal chemistry, vol. 24 no. 27, pp. 3002-3009, DOI: 10.2174/0929867324666170516123206, 2017.

[9] F. Saadeh, S. El Iskandarani, M. Najjar, H. I. Assi, "Prognosis and management of gliosarcoma patients: a review of literature," Clinical Neurology and Neurosurgery, vol. 182, pp. 98-103, 2019.

[10] Y. Zhu, J. M. Sun, Z. C. Sun, F. J. Chen, Y. P. Wu, X. Y. Hou, "MLK3 is associated with poor prognosis in patients with glioblastomas and actin cytoskeleton remodeling in glioblastoma cells," Frontiers in Oncology, vol. 10, 2020.

[11] X. Zhu, S. Liu, X. Yang, W. Wang, W. Shao, T. Ji, "P4HA1 as an unfavorable prognostic marker promotes cell migration and invasion of glioblastoma via inducing EMT process under hypoxia microenvironment," American Journal of Cancer Research, vol. 11, pp. 590-617, 2021.

[12] T. Li, J. Fan, B. Wang, N. Traugh, Q. Chen, J. S. Liu, B. Li, X. S. Liu, "TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells," Cancer Research, vol. 77 no. 21, pp. e108-e110, DOI: 10.1158/0008-5472.CAN-17-0307, 2017.

[13] B. Chen, M. S. Khodadoust, C. L. Liu, A. M. Newman, A. A. Alizadeh, "Profiling tumor infiltrating immune cells with CIBERSORT," Methods in Molecular Biology, vol. 1711, pp. 243-259, 2018.

[14] F. Liang, H. Liang, Z. Li, P. Huang, "JAK3 is a potential biomarker and associated with immune infiltration in kidney renal clear cell carcinoma," International Immunopharmacology, vol. 86, 2020.

[15] S. V. Vasaikar, P. Straub, J. Wang, B. Zhang, "LinkedOmics: analyzing multi-omics data within and across 32 cancer types," Nucleic Acids Research, vol. 46, pp. D956-D963, 2017.

[16] D. Warde-Farley, S. L. Donaldson, O. Comes, K. Zuberi, R. Badrawi, P. Chao, M. Franz, C. Grouios, F. Kazi, C. T. Lopes, A. Maitland, "The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function," Nucleic Acids Research, vol. 38, pp. W214-W220, DOI: 10.1093/nar/gkq537, 2010.

[17] M. Uhlén, L. Fagerberg, B. M. Hallström, C. Lindskog, P. Oksvold, A. Mardinoglu, Å. Sivertsson, C. Kampf, E. Sjöstedt, A. Asplund, I. Olsson, "Proteomics. Tissue-based map of the human proteome," Science, vol. 347, article 1260419, 2015.

[18] M. G. Rees, B. Seashore-Ludlow, J. H. Cheah, D. J. Adams, E. V. Price, S. Gill, S. Javaid, M. E. Coletti, V. L. Jones, N. E. Bodycombe, C. K. Soule, "Correlating chemical sensitivity and basal gene expression reveals mechanism of action," Nature Chemical Biology, vol. 12 no. 2, pp. 109-116, DOI: 10.1038/nchembio.1986, 2016.

[19] Q. Zeng, S. Sun, Y. Li, X. Li, Z. Li, H. Liang, "Identification of therapeutic targets and prognostic biomarkers among CXC chemokines in the renal cell carcinoma microenvironment," Frontiers in Oncology, vol. 9, 2019.

[20] L. Zhou, Y. Li, Z. Li, Q. Huang, "Mining therapeutic and prognostic significance of STATs in renal cell carcinoma with bioinformatics analysis," Genomics, vol. 112, pp. 4100-4114, 2020.

[21] H. Zhang, Z. Dai, W. Wu, Z. Wang, N. Zhang, L. Zhang, W. J. Zeng, Z. Liu, Q. Cheng, "Regulatory mechanisms of immune checkpoints PD-L1 and CTLA-4 in cancer," Journal of Experimental & Clinical Cancer Research, vol. 40 no. 1,DOI: 10.1186/s13046-021-01987-7, 2021.

[22] A. E. Vilgelm, A. Richmond, "Chemokines modulate immune surveillance in tumorigenesis, metastasis, and response to immunotherapy," Frontiers in Immunology, vol. 10, 2019.

[23] L. Zhou, H. Tang, F. Wang, L. Chen, S. Ou, T. Wu, J. Xu, K. Guo, "Bioinformatics analyses of significant genes, related pathways and candidate prognostic biomarkers in glioblastoma," Molecular Medicine Reports, vol. 18, pp. 4185-4196, 2018.

[24] S. W. Huang, N. D. Ali, L. Zhong, J. Shi, "MicroRNAs as biomarkers for human glioblastoma: progress and potential," Acta Pharmacologica Sinica, vol. 39, pp. 1405-1413, 2018.

[25] Q. T. Ostrom, L. Bauchet, F. G. Davis, I. Deltour, J. L. Fisher, C. E. Langer, M. Pekmezci, J. A. Schwartzbaum, M. C. Turner, K. M. Walsh, M. R. Wrensch, "The epidemiology of glioma in adults: a "state of the science" review," Neuro-Oncology, vol. 16, pp. 896-913, 2014.

[26] S. Liu, C. Zhang, N. R. Maimela, L. Yang, Z. Zhang, Y. Ping, L. Huang, Y. Zhang, "Molecular and clinical characterization of CD163 expression via large-scale analysis in glioma," Oncoimmunology, vol. 8, 2019.

[27] A. M. Alves, L. F. Diel, M. L. Lamers, "Macrophages and prognosis of oral squamous cell carcinoma: a systematic review," Journal of Oral Pathology & Medicine, vol. 47, pp. 460-467, 2018.

[28] Z. Xu, L. Wang, J. Tian, H. Man, P. Li, B. Shan, "High expression of B7-H3 and CD163 in cancer tissues indicates malignant clinicopathological status and poor prognosis of patients with urothelial cell carcinoma of the bladder," Oncology Letters, vol. 15, pp. 6519-6526, 2018.

[29] S. Garvin, H. Oda, L. G. Arnesson, A. Lindström, I. Shabo, "Tumor cell expression of CD163 is associated to postoperative radiotherapy and poor prognosis in patients with breast cancer treated with breast-conserving surgery," Journal of Cancer Research and Clinical Oncology, vol. 144, pp. 1253-1263, 2018.

[30] C. Pusztai, M. V. Yusenko, D. Banyai, A. Szanto, G. Kovacs, "M2 macrophage marker chitinase 3-like 2 (CHI3L2) associates with progression of conventional renal cell carcinoma," Anticancer Research, vol. 39, pp. 6939-6943, 2019.

[31] L. Mu, C. Yang, Q. Gao, Y. Long, H. Ge, G. De Leon, L. Jin, Y. E. Chang, E. J. Sayour, J. Ji, J. Jiang, "CD4+ and perivascular Foxp3+ T cells in glioma correlate with angiogenesis and tumor progression," Frontiers in Immunology, vol. 8, 2017.

[32] N. Zhang, Z. Dai, W. Wu, Z. Wang, H. Cao, Y. Zhang, Z. Wang, H. Zhang, Q. Cheng, "The predictive value of monocytes in immune microenvironment and prognosis of glioma patients based on machine learning," Frontiers in Immunology, vol. 12,DOI: 10.3389/fimmu.2021.656541, 2021.

[33] A. L. Hung, R. Maxwell, D. Theodros, Z. Belcaid, D. Mathios, A. S. Luksik, E. Kim, A. Wu, Y. Xia, T. Garzon-Muvdi, C. Jackson, "TIGIT and PD-1 dual checkpoint blockade enhances antitumor immunity and survival in GBM," Oncoimmunology, vol. 7, article e1466769, 2018.

[34] C. Li, Y. Zhou, H. Deng, Y. Ye, S. Zhao, S. Liang, S. Cai, J. Lin, Y. Tang, Y. Wu, "Mining database for the therapeutic targets and prognostic biomarkers among STAT family in glioblastoma," Cancer Biomarkers, vol. 30, pp. 179-191, 2021.