CSF1R is a Prognostic Biomarker and Correlated with Immune Cell Infiltration in the Gastric Cancer Microenvironment

Introduction

Gastric cancer (GC) is the fifth most commonly diagnosed cancer and the third leading cause of cancer mortality worldwide.¹ As patients at the early stage of GC are usually asymptomatic, most of them are diagnosed with advanced stages.^2,3 Although progress has been made in treatments, including surgery, chemotherapy, molecular targeted therapy, and radiotherapy, the prognosis of advanced GC patients remains pessimistic.³ Therefore, it is urgently needed to improve GC diagnosis, treatment, and prognosis.

The tumor microenvironment (TME) has recently attracted more and more attention from the cancer research community. A growing body of evidence has elucidated that the TME could influence gene expression of tumor tissues and is gradually recognized as a critical determinant in tumor development and therapeutic responses.^4–6 The TME refers to the environment where tumor cells originate and develop. In addition to tumor cells, the TME also consists of stromal cells, the tumor vascular system, immune cells, the extracellular matrix (ECM), as well as the acidic and hypoxic environment of the tumor.^7–9 Stromal cells and immune cells, as two major nontumor cell types in the TME, are considered to be of great value in tumor diagnosis and prognosis.^10–13 Thus, an overall understanding of the TME may provide crucial vision into GC initiation and progression, and contribute to the development of promising diagnostic and prognostic biomarkers.

With the advent of the biological big data era, bioinformatics technique provides a more convenient platform for researchers to guide basic experiments.¹⁴ ESTIMATE was developed to evaluate stromal and immune scores for reflecting the level of stromal and immune cells infiltrated in the TME according to the transcriptional profiles of tumor tissues.¹⁵ This algorithm has been applied to colon cancer and glioblastoma, demonstrating the effectiveness of the big-data based algorithm.^16–18 Besides, CIBERSORT algorithm was performed to evaluate the fraction of 22 infiltrating immune cells in the TME using the LM22 signature.¹⁹ CIBERSORT algorithm also has been validated and applied to analyze the association between immune cell landscapes and treatment outcome in lung and breast cancer.^20,21 In this study, the ESTIMATE and CIBERSORT algorithms were utilized to calculate the stromal and immune scores of GC cases and the proportion of 22 immune cells infiltrated in TME based on The Cancer Genome Atlas (TCGA) database. Furthermore, we explored the TME-related differentially expressed genes (DEGs) and found CSF1R correlated with the prognosis of GC and immune cells infiltrated in TME. We hypothesized that CSF1R might play an essential role in the GC microenvironment and could be a promising therapeutic target and prognostic biomarker.

Materials and Methods

Data Preparation

A flowchart of the study is shown in Figure 1. Gene expression data of GC samples were downloaded from the TCGA database. The corresponding patients’ clinical information, including survival time, survival status, grade, and tumor stage, were also obtained from the TCGA database. Then, the ESTIMATE algorithm was applied to calculate the stromal, immune, and ESTIMATE scores for GC samples. The acquirement and application of all data involved in this study were performed following TCGA publication guidelines and data access policies.

Figure 1 Flowchart of this study. Abbreviations: TCGA, the Cancer Genome Atlas; DEGs, differentially expressed genes; PPI, protein-protein interaction; GSEA, gene set enrichment analysis.

Clinicopathological Characteristics Analysis and Survival Analysis

Based on the stromal, immune, and ESTIMATE scores, GC patients were classified into two groups: high-score group and low-score group. We evaluated the association between clinicopathological characteristics and the stromal/immune/ESTIMATE scores.²² Furthermore, the Kaplan-Meier method was utilized to assess the association between the survival rate and the above scores. A Log rank test tested the association. The analyses were carried out through the R packages “survival” and “survminer”, and p < 0.05 was regarded as statistical significance.

DEGs Identification

Gene expression profiles were screened using the R package “limma”, and DEGs between the low-score and high-score groups were identified. The criteria for determining DEGs were an absolute value of log₂fold change (|log₂FC|) >1 and false discovery rate (FDR) <0.05. Heatmaps of DEGs were drawn by the package “pheatmap” in R, and the overlapping DEGs were determined by the package “Venn diagrams” in R.

Enrichment Analysis

Enrichment analysis of Gene Ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway for DEGs were performed by R language loaded with packages “ggplot2”, “enrichplot”, and “clusterProfiler”. The GO analysis showed the feature of DEGs in the biological process, cellular component, and molecular function and the KEGG analysis revealed the enrichment of DEGs in the signaling pathways. P and q <0.05 were considered statistically significant.

PPI Network and Cox Regression Analysis

The PPI network of DEGs was constructed based on the STRING database. Nodes with the confidence of interactive relationship > 0.7 were regarded as the threshold. Subsequently, CytoHubba, which is a plug-in of Cytoscape software, was utilized to identify the top 20 hub genes. Cox regression analysis was utilized to determine prognostic DEGs. DEGs with p < 0.05 for survival comparison were considered as prognostic DEGs. Venn diagrams were used to visualize the common gene CSF1R between prognostic DEGs in univariate Cox regression and the top 20 hub genes in the PPI network.

Gene Set Enrichment Analysis (GSEA) and CIBERSORT Algorithm

GC patients were classified into low-expression and high-expression groups based on CSF1R expression. To explore the role of CSF1R in regulating the GC microenvironment, Gene set enrichment analysis for the KEGG pathway was performed using GSEA version 3.0 downloaded from Broad Institute. KEGG gene sets for enrichment analysis were obtained from the Molecular Signatures Database. NOM p < 0.05 and FDR < 0.25 were regarded as the threshold.

Then, CIBERSORT algorithm was applied to assess the proportion of immune cells infiltrated in the TME using the LM22 signature. LM22 signature consisting of 547 genes was used to identify 22 types of infiltrating immune cells. Furthermore, we performed the difference and correlation analyses to explore the correlation between CSF1R expression and the immune cells infiltrated in the TME (P < 0.05).

Statistical Analysis

Statistical analyses were performed via R software (version 4.0.2) and Perl (version 5.30.1) in this study. All statistical methods and appropriate R packages were described throughout the study. P < 0.05 was set as the cutoff criterion.

Results

The Correlation of Stromal, Immune, and ESTIMATE Scores with Survival and Clinicopathological Characteristics

We downloaded gene expression information and clinical data of 330 GC cases from TCGA. All the patients diagnosed as adenomas or adenocarcinomas were included in our study. The clinicopathological characteristics of GC patients are listed in Table 1. Subsequently, we calculated the stromal, immune, and ESTIMATE scores of GC samples by the ESTIMATE algorithm. The stromal and immune scores ranged from −726.16 to 941.62 and −287.49 to 1418.40, respectively. The ESTIMATE scores, which are the combination of immune and stromal scores, were distributed between −788.22 to 2118.38.

Table 1 Clinicopathological Characteristics of GC Patients from TCGA Database

We compared the stromal/immune/ESTIMATE scores among the GC patients under grade, TNM classification, and stage. A significant correlation was observed between the above scores and grade, stage, and T classification. However, all scores were not associated with N classification and M classification (Figure 2). Furthermore, to evaluate the correlation between the stromal/immune/ESTIMATE scores and survival of GC patients, we divided them into the low-score and high-score groups and compared the overall survival rate in the two groups. Intriguingly, the survival curves revealed that the stromal high-score group had a lower survival rate than the low-score group. Patients with high immune/ESTIMATE scores had a lower survival rate than those with low scores, although no statistical differences were found (Figure 3).

Figure 2 The correlation of stromal, immune, and ESTIMATE scores with clinicopathological characteristics (A–E).

Figure 3 The correlation of stromal, immune, and ESTIMATE scores with the survival of GC patients (A–C). Abbreviation: GC, gastric cancer.

DEGs Identification and Enrichment Analysis

We performed the differential analysis to determine TME-related DEGs. Heatmaps showed the differential gene expression profiles according to the stromal and immune scores (Figure 4A and B). A total of 1055 DEGs were identified between the stromal low-score and high-score groups. Among them, 866 genes were significantly upregulated, while the other 189 genes were downregulated (|log₂FC| >1 and FDR <0.05). Simultaneously, 840 DEGs were obtained from the immune score groups, including 617 upregulated genes and 223 downregulated genes (|log₂FC| >1 and FDR <0.05). Finally, 560 DEGs in both the stromal and immune score groups were identified as TME-related DEGs (Figure 4C and D).

Figure 4 Gene expression profiles of GC patients with different stromal and immune scores. (A) Heatmap of DEGs between the stromal high-score and low-score groups. (B) Heatmap of DEGs between the immune high-score and low-score groups. (C) Venn diagram of commonly downregulated DEGs in the stromal and immune score groups. (D) Venn diagram of commonly upregulated DEGs in the stromal and immune score groups. Abbreviations: GC, gastric cancer; DEGs, differentially expressed genes.

Furthermore, we performed functional enrichment analysis on 560 TME-related DEGs to identify their primary functions. GO enrichment analysis showed that DEGs were significantly enriched in the immune-related terms, including regulation of immune effector process, immune response−activating signal transduction, and immune response-activating cell surface receptor signaling pathway (Figure 5A). The KEGG enrichment analysis also indicated that most of the DEG-related pathways were correlated with immune response (Figure 5B).

Figure 5 Functional enrichment analysis of 560 DEGs. (A) GO enrichment analysis of DEGs. (B) KEGG enrichment analysis of DEGs. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI Network and Cox Regression Analysis

The PPI network was constructed through the STRING database to analyze the interaction of DEGs (Figure 6A). Using the CytoHubba plug-in, we identified the top 20 hub genes by calculating each node’s connection degree in the network (Figure 6B). Univariate Cox regression analysis for GC patients’ survival was conducted to identify 22 prognostic DEGs (Figure 6C). Furthermore, the intersection analysis between the top 20 hub genes in the PPI network and 22 prognostic DEGs was carried out, and CSF1R was found to be the most critical gene (Figure 6D).

Figure 6 PPI network and univariate Cox regression analysis. (A) PPI network of 560 DEGs. (B) The top 20 genes ordered by the number of nodes in PPI network. (C) Forest plot of univariate Cox regression analysis(D) Venn diagram showing the common DEGs shared by the top 20 genes in PPI network and prognostic genes in Cox. Abbreviations: DEGs, differentially expressed genes; PPI, protein-protein interaction.

The Correlation of the Survival and Clinicopathological Characteristics with CSF1R Expression

To evaluate the role of CSF1R in the GC patients’ survival, we classified all samples into low-expression and high-expression groups based on CSF1R expression. The Kaplan-Meier survival curve revealed that GC patients with CSF1R low expression had a better overall survival rate than those with CSF1R high expression (Figure 7). Furthermore, we evaluated the relationship between the expression of CSF1R and clinicopathological parameters. The result revealed that CSF1R expression was associated with grade and stage (Figure 8A and B). However, there were no significant differences in CSF1R expression among the GC patients with T classification, N classification and M classification (Figure 8D and E).

Figure 7 The correlation of CSF1R expression with the survival of GC patients. Abbreviation: GC, gastric cancer.

Figure 8 The correlation of CSF1R expression with clinicopathological characteristics of GC patients (A-E). Abbreviation: GC, gastric cancer.

GSEA

To explore the role of CSF1R in regulating the GC microenvironment, GSEA was performed based on the expression of CSF1R. We found that CSF1R was involved in immune-related signaling pathways, including T cell receptor signaling pathway and natural killer cell mediated cytotoxicity (Figure 9). The results implied that CSF1R might take part in the modulation of immune-active status in the GC microenvironment.

Figure 9 Gene set enrichment analysis of CSF1R.

The Correlation of the Fraction of Infiltrating Immune Cells with CSF1R Expression

The fraction of immune cells infiltrated in the TME was analyzed through the CIBERSORT algorithm (Figure 10). Then, the difference and correlation analyses were performed to explore the correlation between the immune cells infiltrated in the TME and CSF1R expression. The results indicated that a total of 14 kinds of infiltrating immune cells were related to the expression of CSF1R (Figure 11).

Figure 10 Infiltrating immune cells profile in gastric cancer samples and correlation analysis. (A) The proportion of 22 kinds of infiltrating immune cells. (B) The correlation between 22 kinds of infiltrating immune cells.

Figure 11 The association between the proportion of infiltrating immune cells and the expression of CSF1R. (A) Violin plot showed the ratio differentiation of infiltrating immune cells with the expression of CSF1R. (B) Scatter plot showed the correlation of infiltrating immune cells proportion with CSF1R expression. (C) Venn plot displayed 14 kinds of infiltrating immune cells associated with CSF1R expression by difference and correlation analysis.

Discussion

During the past decades, advances in high-throughput sequencing techniques and bioinformatics have improved our understanding of transcriptional alterations in GC. More and more biomarkers have been identified to be correlated with the diagnosis and prognosis of GC. However, the crucial regulators in the GC microenvironment and the change of the TME in GC progression have not been elucidated. In this study, the stromal, immune, and ESTIMATE scores were significantly related to the clinicopathological parameters of GC, such as grade, stage, and T classification. The results revealed that stromal and immune cells infiltrated in the TME might play a crucial role in GC progression. To explore the underlying mechanisms of TME change, we identified TME-related DEGs and performed enrichment analysis, construction of the PPI network, and COX regression analysis. GO and KEGG enrichment analysis indicated that most DEGs were remarkably involved in immune responses. Importantly, CSF1R was also found to be related to the clinicopathological characteristics (stage and grade) and the survival of GC patients. The results implied that CSF1R might play an essential role in modulating the TME and could be a promising therapeutic target and prognostic marker of GC.

CSF1R is a member of the type III protein tyrosine kinase receptor family and plays a central role in many diseases, such as chronic inflammatory diseases, bone disorders, and tumors.²³ CSF1R can be activated either in a paracrine manner or in an autocrine manner, then enhancing the progression and metastasis of tumors.^24,25 For example, the expression of CSF1R and CSF1 in neoplastic epithelial cells was predictive for ipsilateral recurrence and significantly correlated with the survival of breast cancer patients.^26,27 In addition, CSF1R was found to be upregulated in clear cell renal cell carcinoma, and its expression was positively associated with TNM stage, metastasis, and poor survival.²⁸ The elevated expression of CSF1R also was reported in Hodgkin/Reed–Sternberg (HRS) cell lines, and CSF1R signaling was demonstrated to be necessary for the growth of HRS cells.²⁹ Recently, the upregulated expression of CSF1R or CSF1 was greatly correlated with the progression and prognosis of patients with GC. Moreover, high coexpression of CSF1R and CSF1 was also an independent prognostic factor for overall survival.³⁰ Our results also showed that CSF1R was significantly related to the development and survival of GC patients. These findings further implied that CSF1R might be a potential biomarker for GC diagnosis and prognosis.

Furthermore, we further explored the potential role of CSF1R in regulating the TME. GSEA showed that CSF1R participated in immune-related pathways. The results implied that CSF1R might take part in the regulation of immune-active status in the GC microenvironment. Then, using the CIBERSORT algorithm, we estimated the proportion of immune cells infiltrated in the TME. The difference and correlation analyses were performed to explore the correlation between the expression of CSF1R and the immune cells infiltrated in the TME. The results revealed that the expression of CSF1R correlated with several kinds of infiltrating immune cells, including tumor-associated macrophages (TAMs), NK cells, B cells, Neutrophils, Eosinophils, T cells, Dendritic cells, and so on. Recently, several studies have indicated that inhibition of CSF1R can modulate both the number and the function of TAMs in different types of cancer.^31–35 For example, Mitchem et al revealed that targeting TAMs by inhibiting CSF1R decreased the numbers of infiltrating immune cells, relieved immunosuppression, and improved chemotherapeutic efficacy of pancreatic ductal adenocarcinoma.³¹ Pyonteck et al showed that an inhibitor of CSF1R was utilized to target TAMs in a mouse model of proneural glioblastoma multiforme, which significantly regressed established tumors and increased survival.³³ Furthermore, a recent study revealed that CSF1R inhibition decreased the abundance of TAMs within cervical tumors and the associated stroma and delayed the growth of established neoplasms.³⁵ In addition to TAMs, CSF1R can also be detected on other cells within the TME, including neutrophils, dendritic cells, and myeloid-derived suppressor cells.^36–38 All data suggested that CSF1R might be a promising crucial regulator for the microenvironment of GC and participate in regulating immune-active status in TME. However, further studies are needed to identify the exact immune cell atlas in the GC patients’ microenvironment with different CSF1R expression levels.

Conclusion

In summary, we applied the ESTIMATE algorithm to estimate the stromal, immune, and ESTIMATE scores. Then, we concluded that stromal and immune cells infiltrated in the TME might play a crucial role in GC progression. We further explored the TME-related DEGs and found CSF1R associated with GC development and survival. More interestingly, CSF1R might participate in regulating immune-active status in TME and could be a promising biomarker for GC therapy and prognosis. However, more experiments are needed to clarify the potential mechanisms of CSF1R in the microenvironment of GC patients.