Analysis of CSF2RA Expression and Its Functional and Clinical Significance in Gastric Cancer Based on Bioinformatics ()
1. Introduction
Gastric cancer is the fifth most common malignant tumor worldwide, and its high recurrence rate and heterogeneity lead to significant differences in the prognosis of patients [1]. Although the application of targeted therapy (such as anti-HER2 therapy) and immune checkpoint inhibitors (ICI) has improved the prognosis of some patients, molecular heterogeneity remains a key bottleneck in the improvement of efficacy [2] [3]. In recent years, the immunosuppressive properties of the tumor microenvironment (TME) have been confirmed as the core driver of gastric cancer progression [4].
Colony stimulating factor receptor α (CSF2RA) is a key subunit of GM-CSF receptor [5]. It is located in the Xp22.3/Yp11.3 pseudoautosomal region and regulates myeloid cell differentiation by activating the JAK-STAT pathway [6]. Previous studies have suggested that CSF2RA mRNA is significantly upregulated in cancer tissues, such as breast cancer [7], colorectal cancer [8], renal cancer, and liver cancer [9]. CSF2RA promotes immune escape in inflammatory diseases and hematological tumors [9], but its expression pattern, clinical significance and immune regulation mechanism in gastric cancer have not been systematically elucidated. This study was based on a multi-dimensional bioinformatics analysis. The core modules of gastric cancer were identified by TCGA differential expression analysis and WGCNA network, combined with machine learning algorithms to screen key genes, analyze the correlation between CSF2RA expression and clinicopathological features and prognosis, and explore its effect on TME immune cell infiltration, checkpoint molecules, and therapeutic response. The aim of this study was to provide new targets for the prognostic evaluation and immunotherapy strategies of gastric cancer.
2. Methods
2.1. Data Acquisition
We downloaded the GSE54129 and GSE65801 datasets from the GEO database [10] (https://www.ncbi.nlm.nih.gov/geo/), used GEOquery and limma to standardize the data, and merged them after correcting the batch effect with the sva package. RNA-seq (FPKM) and clinical data of GC were downloaded from The TCGA (https://portal.gdc.cancer.gov), Convert FPKM to TPM and standardize log2 (TPM + 1).
2.2. TCGA Data Analysis
Based on the RNA-seq data of TCGA-STAD, differential expression analysis was performed using the DESeq2 package. After standardization, the differential genes were screened with |log2FC| > 1 and P < 0.05. A volcanic map was drawn using the ggplot2 package to show the results of the differential gene analysis.
2.3. Weighted Gene Co-Expression Network Analysis (WGCNA) for Disease-Related Gene Screening
R language was used to analyze the combined data of GSE54129 and GSE65801. The linear model of GC and normal tissues was constructed using lmFit (), and eBayes () was adjusted using empirical bayesian. Differential genes were extracted with |log2FC| > 1 and P < 0.05. Based on the standardized expression data, a co-expression network was constructed using the WGCNA algorithm [11], and the Pearson correlation coefficient between genes was calculated [12]. The soft threshold β was transformed into the adjacency matrix, and the TOM matrix was then calculated. Genes were clustered by dynamic tree cutting algorithm [13] (the module contained 60 genes at least, and the cutting height was 0.25 to merge similar modules), and the characteristic genes were used to characterize the module and analyze its correlation with clinical features [14].
2.4. Intersection of WGCNA Module Genes and TCGA Differentially Expressed Genes (DEGs) with Expression Matrix Extraction
The module genes related to specific clinical features obtained by WGCNA were intersected with the differential genes of TCGA to obtain a set of key genes that were significant in both analyses, and their bulk transcriptome data in TCGA samples were extracted.
2.5. Feature Gene Selection Using Three Machine
Learning Algorithms
In this study, three machine learning methods, LASSO (Least Absolute Shrinkage and Selection Operator), Random Forest and GBM (Gradient Boosting Machine), were used to screen characteristic genes. LASSO [15] is a penalty-based regression method, which shrinks the unimportant feature coefficients to zero by L1 regularization to achieve feature selection. Random forest [16] is an ensemble learning method, composed of multiple decision trees. The tree model is constructed by random sampling and feature selection, and the importance of the features in the tree is screened. GBM [17] is an iterative decision tree algorithm that optimizes the loss function through gradient descent to gradually improve model performance and determine key features. The results were integrated, and the feature gene sets screened by the above three methods were integrated. Finally, the ggplot2 and the VennDiagram packages were used to visualize the results and draw the Venn diagram to obtain the final candidate gene.
2.6. Receiver Operating Characteristic (ROC) Curve Analysis for Evaluating the Diagnostic Value of CSF2RA in GC
The pROC package [18] was used to analyze the diagnostic value of CSF2 RA gene expression in 407 TCGA samples (including 32 normal samples and 375 tumor samples) for gastric cancer. The ROC curve was drawn by ggplot2, and the AUC value was marked. Bootstrap method (1000 times) was used to verify stability.
2.7. CSF2RA Gene Expression Analysis
In this study, the TIMER2.0 (http://timer.comp-genomics.org) online analysis tool) was directly used to analyze the expression of CSF2RA in pan-cancer. We also used GEPIA2 website (http://gepia2.cancer-pku.cn/#index) to draw CSF2RA gastric cancer-specific expression box plot.
2.8. Clinical Correlation Analysis of CSF2RA
The relationship between CSF2RA and patient survival and its clinical significance were further studied. Based on the information of CSF2RA patients in TCGA database, the K-M survival curve was drawn. The Complex Heatmap package of R language was used to draw heat maps to show the relationship between CSF2RA expression level (divided into Low and High) and clinical characteristics such as gender, age, TNM stage and overall stage. The box plot was constructed by ggplot2 package of R language to compare the difference in CSF2RA gene expression level between TNM staging groups.
2.9. Construction of a Nomogram
To establish a robust prognostic tool, we developed a nomogram incorporating independent predictors derived from multivariate Cox proportional hazards regression analysis. The model integrates clinically relevant variables, including age, sex, primary tumor (T) stage, and regional lymph node (N) stage,to provide individualized survival probability estimates for patients with gastric cancer.
2.10. CSF2RA-Based Differential Expression and Functional
Enrichment
In order to dissect the potential functions of CSF2RA-related genes, this study first divided the samples into CSF2RA high-expression and low-expression groups based on the median expression levels of CSF2RA genes in the TCGA-GC samples. Subsequently, the DESeq2 package was applied to directly compare the expression differences of all genes between the high- and low-expression groups to identify differentially expressed genes (|log2FC| > 1, p < 0.05) that were significantly associated with CSF2RA expression status. The screened significant DEGs were included in the subsequent analysis. Gene Ontology (GO) enrichment analysis and KEGG pathway enrichment analysis [19] were performed on this fraction of DEGs using the cluster profiler package with the aim of elucidating the biological processes, molecular functions, cellular components, and key metabolic and signaling pathways that may be involved in the relevant genes in the high CSF2RA expression state.
2.11. Protein-Protein Interaction (PPI) Network Analysis
We used the STRING database [20] to map the CSF2RA protein interaction network.
2.12. Immune Cell Differential Analysis and Immune Cell
Correlation Analysis
The analysis of immune cell difference was completed by R language (4.5.0): the limma package standardized the gene expression matrix; cIBERSORT package [21] estimated the proportion of 22 kinds of immune cells (permutations = 1000); differential cells (FDR < 0.05, |log2FC| > 1) were screened by edgeR package. Based on the data of immune cell ratio, Pearson correlation coefficient was calculated, and the correlation diagram of immune cells was drawn by ggplot2 package.
2.13. Tumor Microenvironment (TME) Analysis
Tumor microenvironment difference analysis was performed in R 4.5.0. The ESTIMATE package [22] calculates the matrix, immunity, and ESTIMATE scores; after grouping according to the median score, the limma package was analyzed for difference (|log2FC| > 1, P < 0.05). To investigate whether CSF2RA affects the prognosis of GC patients by acting on the tumor microenvironment, and whether there is a correlation between the expression of CSF2RA and tumor mutation load.
2.14. Immune Checkpoint and Immunotherapy Response Analysis
We used R language to draw heat maps presenting immune checkpoint-related genes [23] and CSF2RA gene association patterns. For immunotherapy response prediction, the immunophenotypic score (Immunophenoscore, IPS) of gastric cancer (STAD) samples was obtained through the Cancer Immunome Database (TCIA) (https://tcia.at/home), which was calculated based on the expression profiles of immune-related genes in the tumor microenvironment, and has been shown to be effective in predicting patients’ response to ICI therapy [24]. Based on the expression level of CSF2RA, the STAD samples were divided into high-expression and low-expression groups, and the limma package was used to perform difference analysis, combined with the ggpubr package to visualize the difference between the two groups in the therapeutic efficacy of PD-1 inhibitors and CTLA4 inhibitors [25], in order to clarify the effect of the expression level of CSF2RA on the response to immunotherapy.
3. Results
3.1. TCGA Differential Expression Analysis
11,999 differential genes were obtained after differential analysis, and the top 50 genes in terms of significance of the differential analysis results were shown in a volcano plot (Figure 1).
Figure 1. A difference analysis volcano chart.
3.2. WGCNA Co-Expression Network Analysis
In this study, 17350 genes were screened from the combined data set for WGCNA analysis, and the sample tree diagram and trait heat map were drawn (Figure 2(a)). By analyzing the scale independence graph and the mean connectivity graph, β = 6 is selected as the appropriate soft threshold for constructing scale-free networks (Figure 2(b), Figure 2(c)). Based on the TOM matrix, the modules were hierarchically clustered, and similar modules on the clustering tree (Figure 2(d)) were merged to obtain multiple co-expression modules related to Con and Treat traits (Figure 2(e)). The results showed that the MEmagenta module was closely related to the Treat trait (r = 0.59, P = 7e−20), and we selected this module for subsequent analysis. The MEmagenta module gene from WGCNA was crossed with the previously identified DEG and its cross-gene expression matrix was extracted.
3.3. Machine Learning-Based Feature Gene Screening
We used three machine learning algorithms to screen out disease-related characteristic genes. Among them, LASSO regression selected 16 genes (Figure 3(a)), GBM selected 50 genes (Figure 3(b)), and random forest selected 20 genes (Figure 3(c), Figure 3(d)). Moreover, the GBM gene screening criteria is importance > 0. The feature genes obtained by these three algorithms were subjected to intersection analysis, and finally overlapped with a biomarker (Figure 3(e)), which was CSF2RA.
Figure 2. WGCNA co-expression network. (a) Clustering dendrogram of the merged samples; (b) Scale-free fit indices (β) and average connectivity of various soft threshold powers; (c) Histogram of connectivity distributions and scale-free topology at β = 6; (d) Dendrogram of genes clustered by the metric of variance (1-TOM); and (e) Heatmap of the genes’ correlation with clinical traits.
Figure 3. Machine learning-based feature gene selection. (a) LASSO regression coefficient profiles of candidate genes. (b) Top 50 genes ranked by importance in GBM algorithm. (c), (d) Random forest gene importance ranking and error rate curve. (e) Venn diagram of overlapping genes identified by three algorithms, with CSF2RA as the sole intersection.
3.4. Receiver Operating Characteristic (ROC) Curve Analysis
The diagnostic value of CSF2RA was evaluated by ROC curve. The area under the ROC curve (AUC) of CSF2RA was 0.58 (Figure 4).
3.5. CSF2RA Gene Expression Profiling
Pan-cancer expression analysis of CSF2RA in the TIMER2 database showed that CSF2RA was differentially expressed in 13 cancers containing GC (Figure 5(a)). The expression of CSF2RA in 619 TCGA samples showed that the expression of CSF2RA in tumor tissues was higher than that in normal tissues (P < 0.5) (Figure 5(b)).
3.6. Association between CSF2RA and Clinicopathological
Features in GC
K-M survival analysis the overall survival rate of patients with high CSF2RA expression was lower than that of patients with low CSF2RA expression in GC patients (P < 0.01, Figure 6(a)). Heatmap analysis (Figure 6(b)) showed that the
Figure 4. ROC evaluating the diagnostic value of CSF2RA in GC (AUC = 0.58).
Figure 5. CSF2RA gene expression analysis. (a) CSF2RA expression in pan-cancer; (b) CSF2RA GC-specific expression box plot.
Figure 6. Correlation between CSF2RA and clinicopathological characteristics in gastric cancer. (a) Kaplan-Meier survival analysis; (b) Heatmap of clinical parameter correlations; (c)-(f) CSF2RA expression across TNM staging groups.
high expression of CSF2RA was significantly correlated with clinical stage and T stage, suggesting that CSF2RA may play an important role in GC progression. There were significant differences in CSF2RA expression between the different TMN stages and the overall clinical stages (Figures 6(c)-(f)).
3.7. Construction of the Prognostic Nomogram
Based on a multivariate regression analysis incorporating independent risk factors as predictive factors, a nomogram for overall survival (OS) in GC patients was developed (Figure 7(a)). Discrimination performance analysis of the model showed that the C-index of the nomogram was 0.665 (95% CI: 0.62 - 0.71), indicating moderate predictive discriminatory ability. After curve validation, the results indicated that the model’s predicted probabilities were highly consistent with the actual occurrence probabilities, suggesting that the nomogram has reliable predictive capability for OS across different time periods (Figure 7(b)). The two sets of results mutually corroborate each other, collectively supporting the model’s clinical applicability.
Figure 7. Prognostic nomogram for GC patients. (a) Nomogram integrating independent predictors (CSF2RA, TNM stage, etc.) for 1-/3-/5-year overall survival (OS). (b) Calibration curve of the predictive model.
3.8. Enrichment Analysis
The results of GO and KEGG enrichment analysis (Figure 8(a), Figure 8(b)) showed that these genes were mainly enriched in immune-related functions, covering immune response processes such as leukocyte-mediated immunity and lymphocyte-mediated immunity. It involves immune pathways such as hematopoietic cell lineage, cytokine-cytokine receptor interaction, etc. It is also associated with immune molecular functions such as antigen binding and immune receptor activity, which provides a functional basis for exploring the mechanism of CSF2RA in gastric cancer.
Figure 8. Enrichment analysis of CSF2RA-associated genes. (a) Gene Ontology (GO) terms; (b) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.
3.9. PPI Network Analysis
The protein-protein interaction network interacting with CSF2RA in Homo sapiens was obtained using the STRING database (Figure 9(a)), and then the data were imported into Cytoscape software for visual analysis of the protein-protein interaction map (Figure 9(b)). Genes that interact with CSF2RA include CSF2RB, JAK2, and IL3RA. Further studies can be carried out according to the specific mechanism of their interaction with CSF2RA.
3.10. Correlation Analysis of CSF2RA Expression with Immune Cell Infiltration
Immune cell infiltration and correlation analysis was carried out using the CIBERSORT algorithm combined with the R package. Box line plots showed (Figure 10(a)) that after grouping based on median CSF2RA expression, the distribution of multiple immune cell subpopulations differed between high and low expression groups (P < 0.05). In correlation analyses, CSF2RA expression was positively correlated with activated memory CD4⁺ T cells (R = 0.2, P = 0.018) (Figure 10(c)), M2-type macrophages (R = 0.22, P = 0.0096) (Figure 10(f)); with regulatory T cells (Tregs, R = −0.21, P = 0.015) (Figure 10(b)), plasma cells (R = −0.22, P = 0.0093) (Figure 10(d)) were negatively correlated, and neutrophils also showed a weak positive correlation (R = 0.18, P = 0.037) (Figure 10(e)), and the scatter plots visualized these correlation trends, suggesting that CSF2RA expression may be involved in the relevant biological processes through the modulation of immune cell infiltration .
3.11. Tumor Microenvironment (TME) Analysis
Analysis of CSF2RA, after grouping by median expression, the violin plots showed (Figure 11(a)) that the StromalScore, ImmuneScore, and ESTIMATEScore of the high-expression group were significantly higher than those of the low-expression group (all P < 0.001), suggesting that it was correlated with the elevated tumor microenvironment (TME) score; Correlation analysis showed that CSF2RA expression was significantly negatively correlated with tumor mutational load (TMB) (R = −0.16, P = 0.0026) (Figure 11(b)).
3.12. Immune Checkpoint and Immunotherapy Response Analysis
To investigate the predictive value of CSF2RA on immunotherapy response in gastric cancer patients, CSF2RA and immune checkpoint association analysis was carried out in this study. The correlation heatmap (Figure 12(a), Figure 12(b)) showed that CSF2RA was associated with many immune checkpoint genes, and the dominant red line indicated that positive correlation was the main pattern, such as the CD family on the surface of immune cells [26] (CD80, CD274, etc.), and members of tumor necrosis factor receptor superfamily [27] (TNFSF9, etc.), which suggests that CSF2RA may be involved in the immune gene regulatory network. Correlation heatmaps showed that CSF2RA was mostly positively
Figure 9. PPINnetwork of CSF2RA. (a) STRING database-derived interactions. (b) Cytoscape visualization highlighting core interactors (CSF2RB, JAK2, IL3RA).
Figure 10. Analysis of CSF2RA immune infiltration in GC. (a) Difference in immune cell expression in GC between CSF2RA high expression group and low expression group. (b)-(f) Correlation between CSF2RA expression level and immune cell infiltration level.
Figure 11. Tumor microenvironment (TME) analysis. (a) Correlation of CSF2RA expression with tumor microenvironment score; (b) Correlation of CSF2RA expression with tumor mutation load.
Figure 12. Immune checkpoint and therapy response. (a), (b) Heatmaps of CSF2RA correlations with checkpoint genes (c), (d) Analysis of CSF2RA expression and immunotherapy.
correlated with immune checkpoint genes. The results of violin box line plot showed that the low CSF2RA expression group had a good therapeutic effect when treated with CTLA4 inhibitor alone (P < 0.01, Figure 12(c)); and the low CSF2RA expression group had a good therapeutic effect when treated with neither PD-1 inhibitor nor CTLA4 inhibitor (P < 0.001, Figure 12(d)).
4. Discussion
CSF2RA is the α-subunit of the granulocyte-macrophage colony-stimulating factor (GM-CSF) receptor [5], which is located in the pseudoautosomal region of the X/Y chromosome [28], and plays a key role in hematopoiesis, immune regulation and inflammatory response. In recent years, its role in the tumor microenvironment [9] has garnered increasing attention. In this study, we revealed for the first time that CSF2RA is significantly highly expressed in GC by combined multi-omics analysis, and its expression level is tightly correlated with tumor progression and poor prognosis. This finding is consistent with the pro-carcinogenic [29] role of the GM-CSF signaling pathway in solid tumors: previous studies have shown that GM-CSF promotes tumor cell proliferation through activation of the JAK/STAT pathway [6], whereas CSF2RA, as its receptor subunit, may be involved in the malignant phenotypic regulation of GC by mediating downstream signaling. In addition, the MEmagenta module, which CSF2RA was located, was highly associated with the GC phenotype in the WGCNA, further supporting its central position in GC development.
Differential expression analysis (TCGA) combined with WGCNA co-expression network revealed that the MEmagenta module, where CSF2RA is located, is highly correlated with gastric cancer phenotype. Further cross-validated by three machine learning algorithms (LASSO, GBM, Random Forest), CSF2RA was identified as the only overlapping signature gene, highlighting its reliability as a gastric cancer biomarker. Despite its limited diagnostic efficacy, it has potential as a core component of multigene combination markers.
The expression of CSF2RA in gastric cancer tissues was significantly higher than that in normal tissues, and was further elevated in advanced clinical staging and T/N staging. Survival analysis confirmed that the overall survival of patients with high CSF2RA expression was significantly shorter, suggesting that it could be used as an independent prognostic factor. The CSF2RA high-expression group exhibited higher stromal scores, immunity scores, and ESTIMATE scores, suggesting that it was closely associated with immune-infiltrating TME [30].
Immune infiltration analysis revealed that CSF2RA was positively correlated with M2 macrophages and activated CD4⁺ T cells, whereas it was negatively correlated with regulatory T cells and plasma cells, and CSF2RA may create an immune escape microenvironment by synergistically promoting pro-tumorigenic M2 polarization [31] and suppressing plasma cell-mediated humoral immunity [32]. By GO/KEGG analysis, CSF2RA-related genes were significantly enriched in immune response regulation, such as leukocyte-mediated immunity [33]; and cytokine-receptor interactions, consistent with GM-CSF signaling function, promoted myeloid immune cell activation [34], corroborating its immune regulatory function. According to PPI network analysis, CSF2RA and CSF2RB may promote JAK2 activation by forming a complex [35], which activates the JAK-STAT pathway pathway to regulate immunity [6]; in the gastric cancer microenvironment, the complexes may abnormally activate or promote tumor-associated macrophage (TAM) polarization toward the pro-tumor M2 type [36], contributing to tumor immune escape and angiogenesis. Our results showed a negative correlation between risk scores and TMB in a prognostic model of gastric cancer, implying that patients in the high-risk group were not as effective as patients in the low-risk group in receiving immunotherapy [37]. CSF2RA was positively correlated with several immune checkpoint genes [23] (CD274 (PD-L1), CD80, TNFSF9, etc.), suggesting that it may be involved in immune escape. In the anti-PD-1/CTLA-4 treatment model, it was found that the efficacy was better in the CSF2RA low-expression group when treated with neither PD-1 nor CTLA4 inhibitors, and when treated with only CTLA4 inhibitors. This is informative for immunotherapy in patients with advanced GC suggesting that it may be predictive of immunotherapy response. The column-line graph model constructed in this study integrates CSF2RA expression and TNM stage, which can accurately predict the 1/3/5-year survival rate of patients (with good calibration curve fit) and provide a tool for individualized prognostic assessment.
In summary, this study confirmed that CSF2RA drives GC progression by remodeling the immunosuppressive microenvironment, and is a potential prognostic marker and immunotherapy target. In the future, it is necessary to further explore: the specific mechanisms by which CSF2RA regulates immune cell infiltration (e.g., through cytokine signaling); and its application value in targeted therapy or immune-combination therapy for gastric cancer.
Funding
Guangxi Natural Science Foundation (Grant No. 2025GXNSFHA069094);
Baise Science and Technology Program (“Technology Foundation” Initiative Special Project, Grant No. 2025ZJ0707);
National College Students innovation and entrepreneurship training program (202510599015);
National College Students innovation and entrepreneurship training program (S202410599049).