Detection of immune-related cells and genes in the glioma tumor microenvironment from the CGGA database using bioinformatics strategies

Background: Glioma is one of the most common and aggressive types of primary malignant tumors in the neuroectoderm, with high mortality rates, and recurrence. Immune checkpoint inhibitors have provided new insights in the treatment of tumors. Yet, there are still a lack of effective immune markers and therapeutic targets for glioma immunotherapy. Material and Methods: First, The data of glioma and corresponding clinical characteristics obtain from the CGGA database. Secondly, we calculated the abundance of immune cells for each sample by the CIBERSORT algorithm and then identied immune cells associated with survival by Kaplan-Meier survival analysis. Critical immune genes were identied naturally by differential gene analysis, functional enrichment analysis, protein interaction analysis, and co-expression analysis. Eventually, the TIMER database was used to analyze further the level of immune inltration of key immune cells. Results: Here, we have discovered macrophages M0 and monocytes, are associated with the survival of glioma, and three immune genes associated with survival were also uncovered, of which BIRC5 is an immune gene-specic to monocytes. The BUB1B and NCAPH are immune genes specic to macrophages M0. Furthermore, three immune genes were positively correlated with the inltration levels of immune cells such as purity and dendritic cell and had a high response to anti-PD-1 treatment. Conclusions: In conclusion, we have identied macrophages M0 and monocytes, and three specic immune genes that may be of high research value in glioma immunotherapy. Our work may provide new research insight into the diagnosis and prognosis of glioma.


Introduction
Glioma is a primary intracranial cancer that is the most common and predominantly aggressive [1]. It appears to be a malignant tumor of several cellular neoplasm origins with neural stem cell-like characteristics and is characterized by high mortality and recurrence [2]. Radiotherapy after surgery combined with temozolomide is currently the main treatment for patients with glioma [3]. Nevertheless, the prognosis of glioma remains disappointing, with an average 5-year survival rate of 4.7% [4], and median overall survival (OS) of 14.6 months. In 2021, the World Health Organization (WHO) classi cation of central nervous system (CNS) tumors rede nes glioma classi cation by adding histologic classi cation, grade, and molecular genetics. It introduces new tumor types and subtypes that provide a basis for better understanding the prognosis and optimally treating patients with speci c CNS cancers [5]. To date, it remains that no effective treatment or diagnosis of glioma is available.
With the accelerated advancement of next-generation sequencing technology, the volume of bioinformatics data has exposed an explosive growth [6]. Great progress has been made in medical research, especially in tumor research [7]. Recently, tumor immunology has rapidly developed, while immune checkpoint inhibitors are also practiced in many tumor treatments [8], and the improvement of diagnosis and treatment of glioma is promising. Immunotherapy and clinical trials targeting glioma are also conducted on a large scale, but whose outcome is not remarkable. Since tumor immune escape is universal, tumor immunotherapy can also be mainly used to personalize diagnostic treatment for patients. To date, immunotherapy for glioma remains highly challenging, and prediction of immunotherapeutic biomarkers and the search for new effective therapeutic targets have been a critical concern.
The tumor microenvironment (TME) is the eld that tumor cells interact with immune cells and immunoregulatory substances in the immune system [9]. Various immune cell subsets are recruited into the tumor microenvironment through recognition of chemokines and chemokine receptors and play a remarkably essential role in tumorigenesis [10]. They also in uence the diagnosis and treatment of tumors [11]. Immune evasion ability is expected to provide novel insight into cancer diagnosis and treatment as a new feature to utilize cancer cells. Immune checkpoint inhibitors, represented with anti-CTLA-4 and anti-PD antibodies, and CAR-T, considered a form of secondary immune cell therapy, play an essential part in anti-tumor treatment [12]. TME plays an in uential role in the ght against tumors and immune cells and provides an opportunity for the immune escape of tumor cells [13]. Therefore, the effect of immune elements in TME on the expression of tumor cell genes and the effect of clinical drugs can also in uence the overall therapeutic outcome of tumors.
Cancer immunotherapy mainly involves the mutual recognition between immune substances in the tumor microenvironment and tumor cells to play a role in removing cancer cells by stimulating and enhancing the antibody immune response of the body [14]. Currently, targeting immunotherapy such as CTLA-4, PD-1, and PDL1 is the immunotherapy of tumors that has achieved great achievements [15]. In the current report, Two survival-associated immune cells and three immune cell-speci c prognosis-associated genes were identi ed and key immune cells and biomarkers may serve as new immunotherapeutic targets for glioma.

Materials And Methods
Data capture and advanced processing RNA-Seq gene expression pro les of glioma patients were rst downloaded from the Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/) database. The clinical data, mainly on gender, age, tumor grade, clinical stage, and survival time, were also downloaded.

The identi cation of immune cells with prognostic signi cance
The CIBERSORT is a tool for deconvolution of expression matrices of human immune cell subtypes using principles from linear support vector regression [16]. The method is based on a known reference set of gene expression features for 22 immune cell subtypes. The calculation of immune cell abundance was performed as follows. In the present study, we ran CIBERSORT in local R software to calculate the abundance ratio matrix of immune cells in gliomas. P ≤ 0.05 was used as a criterion for ltering out samples. The correlation analysis was then performed on the abundance of immune cells in the ltered samples. At the same time, Kaplan Meier analysis of the overall survival of immune cells in 22 was performed based on their abundance, with the median as the cut-off level, and tested with Log-Rank. Then, immune cells associated with survival were identi cation.

The recognition of candidate immune cells in association with clinical information
The abundance of candidate immune cells was analyzed in the 159 samples with the WHO classi cation, age (>45 and <45), and IDH mutation, in combination with immune cell abundance and clinical characteristics. To explore the relationship between candidate immune cells and clinical, independent samples t-test was used for two variables, while an one-way ANOVA test was used for variables with item variables >2.
The identi cation of immune-related genes Candidate immune cells were selected and classi ed into high and low groups according to their median. To obtain the differentially expressed genes of the candidate immune cells, we used the Limma R language package for performing differential expression gene analysis. Here, criteria for differential analysis were set to the condition of |logFC|>0.8 and P<0.05.

Functional annotation of immune-related genes
The Database for Annotation, Visualization, and Integrated Discovery [17] (DAVID, https://david.ncifcrf.gov/) and Metascape (http://metascape.org/gp/) are the most frequently used tools for GO term analysis and classi cation and gene enrichment analysis, respectively. GO terms and KEGG pathways were identi ed by biological process (BP), molecular function (MF), and cellular component (CC). P<0.05 was used as a criterion for screening. M is the number of all genes annotated for a speci c GO term; M is the GO term annotation for a certain number of differentially expressed genes.
The calculated p-values were Bonferroni-corrected to a p-value of 0.05 threshold, and eligible GO terms were de ned as GO terms signi cantly enriched in differentially expressed genes, and the main biological functions of exercising differentially expressed genes could be determined by GO function signi cant enrichment analysis.

The construction of PPI network and identi cation of hub genes
To better understand protein-protein interactions in immune-related genes, a PPI interaction network was established using the STRING web server [18] (https://cn.string-db.org/). PPI node pairs with a combined score ≥ 0.4 were selected for further analysis. The cytoHubba plugin [19] in Cytoscape software was used to identify the genes with the highest connectivity as hub genes for identi cation.
is the product of all positive integers containing v and (i -1)! is the product of all positive integers that are less than CI. If there are no edges between neighbors of node v, then MCC(v) is equal to its degree.

Detection and validation of immune prognosis-related genes
To investigate the prognostic effect of hub genes, we combined the dataset of hub genes with clinical data and used one-way Cox proportional risk regression analysis (hazard ratio (HR) ≠ 1, p<0.05) to screen for genes associated with prognosis. The training set was divided into high and low expression groups according to the expression of each gene. The survival R language package was used for performing Kaplan-Meier (K-M) analysis and a one-sample log-rank test to show the relationship between candidate gene expression and OS in patients with glioma. The nal prognosis-related gene selection genes were obtained.

( )( ) ( )
Immunological correlation analysis of hub immune-related genes Tumor Immunity Estimation Resource [20] (TIMER, http://timer.cistrome.org/) allows for a systematic and comprehensive analysis of immune in ltration levels of immune cells in different cancers. In our work, TIMER was applied to detect the correlation of key immune gene expression in six types of in ltrating immune cells, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells in glioma. Simultaneously, the Submap algorithm was performed to predict the clinical response to PD-1 and PD-L1 immune checkpoint blockade therapy in high and low expression subgroups of key immune genes [21]. P-value < 0.05 was statistically signi cant.
Identi cation of small-molecule drugs with hub immunerelated gene effect The CMAP database (https://portals.broadinstitute.org/CMap/) consists of 7056 gene expression pro les induced by 1309 small molecules widely used to explore the potential uncharted effects of existing drugs on the disease [22]. In this investigation, the drug development tool was exploited to mine small molecule drugs that worked with critical genes in glioma. A cutting criterion of p <0.05 was adopted. Enrichment scores (from -1 to +1) were calculated to evaluate the similarity between genes and drugs. When enrichment scores were >0 were considered drugs with potential synergistic effects on glioma, indicating that they could mimic the biological state of cancer cells. On the contrary, the enrichment score < 0 indicates a potentially antagonistic effect of small molecule drugs, which are considered to reverse the biological process of cancer cells and can be used as therapeutic agents.

Statistics Analysis
In this research, R (version: 4.1.0) software has been used for statistical tests and graph drawing.
Wilcoxon test has been used for statistical differences between the two groups, and Kruskal-Wallis H test has been used for comparing multiple groups. Statistical signi cance was de ned as p < 0.05.

Results
Data sources and preprocessing Acquiring the immune cells associated with the prognosis of glioma After CIBERSORT ltering by the deconvolution algorithm, we got 159 from 693 samples. The proportion of 22 immune cell types in each sample was also calculated ( Figure 1A, B). Monocytes, type 2 macrophages (M2), and mast cells activated in the TME of gliomas represent a higher proportion of all cells, while the percentage of cells such as B cells, T cells CD4 naïve, and Eosinophils is lower. Subsequently, the Kaplan-Meier analysis was used to correlate the abundance ratios of 22 immune cell types with overall patient survival. It was discovered that the abundance ratios of monocytes (p = 0.15 × 10 −3 ) and macrophages M0 (p = 0.16× 10 −1 ) were correlated with patient survival (Figure 1C, D).

Identi cation of critical immune-related genes in macrophages M0
When we identi ed macrophages M0 as a critical immune cell in glioma, which we analyzed for correlation with clinical characteristics, we found is that macrophages M0 correlates strongly with WHO mutation, methylation, and classi cation While M0 was not statistically correlated with codeletion, gender and primary recurrence status. (Figure 2A). Subsequently, 548 up-regulated immune-related genes and 289 down-regulated immune-related genes were identi ed in differential analysis of macrophages M0 sub-cell population ( Figure 2B). Gene function annotation demonstrated that immune-related genes were majorly involved in positive regulation of cell proliferation, leukocyte migration, cell division, and other functions ( Figure 2C, Supplement Table 2); KEGG pathway enrichment analysis revealed that immunerelated genes were mainly involved in Ras signaling pathway, PI3K-Akt signaling pathway, p53 signaling pathway, and other classical biological pathways of glioma ( Figure 2D, Supplement Table 3). We performed protein-protein interaction network analysis to explore the interactive relationship of immunerelated genes and to get hub genes. The most connected network nodes were identi ed by the Cytoscape plugin cytoHubba, we obtained 22 hub genes, their connectivity was consistent (MCC score = 1.93E+70).

Identi cation of critical immune-related genes in monocytes
The monocytes are another critical immune cell for glioma that we identi ed. We found a strong correlation between monocytes and WHO mutations ( Figure 3A). Next, 246 up-regulated immune-related genes and 646 down-regulated immune-related genes were identi ed during differential analysis of monocytes ( Figure 3B). Gene function annotation showed that immune-related genes mainly function with cell proliferation, mitochondrial nuclear division, in ammatory response, etc. ( Figure 3C, Supplement Table 4). KEGG pathway enrichment analysis demonstrated that differential immune-related genes were mainly involved in biological pathways such as Rap1 signaling pathway, PI3K-Akt signaling pathway, etc.
( Figure 3D, Supplement Table5). As similar analysis as macrophages M0 for monocytes was used to identify the 21 hub genes (MCC score = 5.51E+67) with the greatest connectivity ( Figure 3E, F).

Recognition immune-related genes speci c to macrophages M0 and monocytes
To further identify immune-related genes speci c to macrophages M0 and monocytes, we overlapped the identi ed core immune-related genes in two immune cells ( Figure 4A), and there was a large overlap of critical immune-related genes in both classes of immune cells A total of 23 genes, 20 of which are intersected, and only 3 of which are speci c to both immune cell classes. We discovered that BIRC5 is an immune-related gene speci c to monocytes. BUB1B and NCAPH are crucial immune-related genes speci c to macrophages M0. Gene function enrichment analysis showed that all these hub genes are mainly involved in cell division, mitochondrial cell cycle process, meiotic nuclear division, DNA replication, and other functions ( Figure 4B, Supplement Table6). We performed a Kaplan-Meier survival analysis of immune-related genes speci c to two classes of immune cells. We noted a signi cant correlation between these three genes (BIRC5, BUB1B, and NCAPH) and clinical signatures ( Figure 4C-D).
Critical immune-related genes are involved in tumor immune in ltration in glioma Three hub immune-related genes are all potential immunotherapeutic targets and are important for immune-related studies of glioma. Meanwhile, we examined the correlation between three immune cellspeci c survival-related genes and the in ltration level of immune cells by TIMER, and the results are shown in Figure 5A. The expression levels of BIRC5 were positively correlated with the in ltration levels of purity, dendritic Cell. The expression levels of NCAPH were positively correlated with the in ltration levels of purity, neutrophil. The expression levels of NCAPH were positively correlated with the in ltration levels of purity, neutrophil, and dendritic cell. The expression levels of BUB1B were positively correlated with the in ltration levels of purity, CD8+ T dell, and dendritic Since PD1 and PD-L1 play an essential role in tumor immunosuppression and immunotherapy, we investigated the relationship between the expression levels of three hub immune-related genes and PD1 and PD-L1, respectively. We found that the expression levels of NCAPH and BIRC5 were signi cantly in positive correlation with the expression of PD-L1, and the expression levels of BUB1B were signi cantly correlated with the expression of PD-1 and PD-L1. These data indicate that three key immune genes can be used as accurate predictors of immune features of glioma ( Figure 5B).

CMAP analysis identi es small molecule drugs for glioma treatment
To acquire small molecule drugs that interact with a key immune gene, we uploaded three key immune genes to the CMAP website. Overall, eight small molecule drugs could be potential drugs for the treatment of glioma patients by shear criteria.

Discussion
Glioma is a malignant tumor of the central nervous system with a high mortality rate and a high risk of recurrence [23]. The treatment for glioma includes surgical resection, adjuvant radiotherapy, chemotherapy combined with temozolomide, etc. [24]. Nevertheless, the lack of precise understanding of the pathogenesis of glioma has resulted in a low survival rate. Recently, new therapeutic approach immune checkpoints have been discovered, which provide new ideas for immunotherapy of tumors by speci cally blocking immunosuppressive effects [25]. In the present work, we screened and identi ed cells and genes signi cantly associated with immune in ltration and clinical features in the glioma microenvironment by bioinformatics approach. This research provides a new idea for immunotherapy of glioma.
In the present research, we discovered that two types of immune cells, macrophages M0 and monocytes, are closely associated with the survival of patients with glioma. During immune editing, immune cells have the function of withdrawing and evoking tumorigenic responses [26]. Macrophages are the most abundant cells in the tumor microenvironment and are highly plastic, making them perform multiple functions in the tumor microenvironment [27]. Macrophages are commonly de ned as two subtypes activated and inactivated, where activated macrophages include classically activated M1 cells and alternatively activated M2 cells [28]. The other subtype is an inactivated macrophage called M0, and this type of cell typically is neither stimulated nor receives signals that promote activation and functional polarization [29]. Generally, the macrophages M0 can be reprogrammed to polarized macrophages once the macrophages M0 receive M1 or M2 speci c signals and in turn shift to an anti-in ammatory M1 phenotype, including radiotherapy and immune checkpoint blockade, such as anti-PD-L1/PD-1 strategies [27]. Monocytes are innate immune cells of the mononuclear phagocyte system and have an in uential role in regulating cancer development and progression [30]. It effectively suppresses tumor-speci c T-cell immunity by secreting tumor necrosis factor and stimulating monocytes to express PD-L1, PD-L1+ monocytes promote tumor growth in vivo [31]. In summary, monocytes and macrophages M0 are signi cantly associated with survival in patients with glioma and probably play an invaluable role in immune in ltration and immunotherapy.
GO functional annotation and KEGG pathway enrichment analysis of macrophages M0 and monocytes showed shared involvement in classical immune pathways such as PI3K-Akt signaling pathway, p53 signaling pathway, and ECM-receptor interaction. In the tumor microenvironment, immune cells interact through cell-cell communication 31063754. the PI3K-Akt signaling pathway, a classical signaling pathway in tumorigenesis, has an important role in activating numerous immune cells [32]. The PI3K-AKT-mTOR signaling pathway has been reported to contribute to the development and activity of lymphocytes [33]. The mutation of the tumor suppressor p53 is critical for tumor evasion of apoptosis and senescence, and p53 dysfunction also promotes in ammation and supports tumor immune evasion, thereby serving as an immunological driver of tumorigenesis [34]. Targeting the p53 pathway in TME can be used to reverse immunosuppression, enhance the antitumor killed effect, and harness tumor-speci c, persistent, and systemic antitumor immunity with minimal toxicity [34]. ECM-receptor interaction is instrumental in the process of tumor cell shedding, adhesion, degradation, motility, and proliferation [35]. It has been reported that ECM can promote epithelial-mesenchymal transition (EMT) of cancer cells in colorectal cancer [36]. The ECM can inhibit angiogenesis in glioma through endothelial-macrophage interactions and altered extracellular matrix in glioma [37]. Meanwhile, Platelet activation, Rap1 signaling pathway, in ammatory response, leukocyte migration, and other functions and signaling pathways are very tightly linked to tumor immunity [38,39].
Ultimately, we characterized three cell-speci c immune-related genes, with BIRC5 being a monocytesspeci c immunity gene and BUB1B and NCAPH being macrophages M0-speci c immune genes. Immunological associations with the 3 genes have not been reported in previous studies of gliomas. The BIRC5 (Survivin), which belongs to the family of apoptosis suppressor genes, is abundantly expressed in tumors [40]. It has been investigated that Brexpipazole can increase the sensitivity of glioma stem cells to osimertinib by reducing the expression of BIRC5 [41]. BUB1B is a mitotic spindle checkpoint molecule. It is associated with tumor cell aggressiveness and drug resistance in gliomas [42]. NCAPH is a structural component of chromosomes during mitosis [43]. High expression of NCAPH promotes the proliferation of colorectal cancer cells and increases radiation resistance, and NCAPH gene silencing stops the cell cycle at the G2/M phase, subsequently inducing apoptosis. In prostate cancer, knockdown of NCAPH decreased cell proliferation and inhibited the aggressiveness of prostate cancer cells.
In conclusion, we identi ed two immune cells and three speci c immune genes associated with survival in glioma patients and demonstrated that the expression of three immune genes correlated with the overall survival of glioma patients in the present study. Macrophages M0 and monocytes, two classes of immune cells, and three immune genes, BIRC5, BUB1B, and NCAPH, can be used as prognostic glioma markers and guide immunotherapy glioma. The results of this work are mainly from bioinformatics experiments and remain to be further investigated by cell biology experiments. This study provides new and more valuable ideas to study the tumor microenvironment and immunotherapy and prognosis of glioma.

Declarations Author Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
Publicly available datasets were analysed in this study. The gene expression RNA-Seq of glioma from the Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/).

Competing interests
The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest.
Funding    Supplementary Files