Identification of macrophage-related molecular subgroups and risk signature in colorectal cancer based on a bioinformatics analysis

Abstract

Macrophages play a crucial role in tumor initiation and progression, while macrophage-associated gene signature in colorectal cancer (CRC) patients has not been investigated. Our study aimed to identify macrophage-related molecular subgroups and develop a macrophage-related risk model to predict CRC prognosis. The mRNA expression profile and clinical information of CRC patients were obtained from TCGA and GEO databases. CRC patients from TCGA were divided into high and low macrophage subgroups based on the median macrophage score. The ESTIMATE and CIBERSORT algorithms were used to assess immune cell infiltration between subgroups. GSVA and GSEA analyses were performed to investigate differences in enriched pathways between subgroups. Univariate and LASSO Cox regression were used to build a prognostic risk model, which was further validated in the GSE39582 dataset. A high macrophage score subgroup was associated with poor prognosis, highly activated immune-related pathways and an immune-active microenvironment. A total of 547 differentially expressed macrophage-related genes (DEMRGs) were identified, among which seven genes (including RIMKLB, UST, PCOLCE2, ZNF829, TMEM59L, CILP2, DTNA) were identified by COX regression analyses and used to build a risk score model. The risk model shows good predictive and diagnostic values for CRC patients in both TCGA and GSE39852 datasets. Furthermore, multivariate Cox regression analysis showed that the risk score was an independent risk factor for overall survival in CRC patients. Our findings provided a novel insight into macrophage heterogeneity and its immunological role in CRC. This risk score model may serve as an effective prognostic tool and contribute to personalised clinical management of CRC patients.

Keywords:

• Colorectal cancer
• macrophage-related genes
• risk score model
• tumour immunity
• macrophage-related subtypes

Introduction

Colorectal cancer (CRC) is one of the most malignant tumours and has a major impact on the normal lives of millions of people [1]. Approximately 25% of CRC patients have stage IV tumours. There are still some patients who can be diagnosed at an early stage, but the cancer is still metastasises [2]. According to statistics, there are more than 2.2 million new CRC patients worldwide, and the incidence is increasing [2]. Despite advances in the treatment of CRC in recent years, the 5-year survival rate of patients remains low [3–4]. Currently, surgery is the main treatment for CRC, but due to the high risk of recurrence and metastasis after surgery, patient survival is generally low [5]. Therefore, to improve the prognosis of CRC patients, it is necessary to identify potential prognostic markers for the therapeutic targets and prognostic assessment of CRC.

Tumour-infiltrating immune cells are involved in the initiation and development of CRC. Alteration of T cell phenotypes may not only affect the regulation of chronic mucosal inflammation but also colitis-associated carcinogenesis [6]. T cell receptor stimulation and co-stimulatory signalling enhance T cell activation, offering a promising therapeutic option for immunotherapy in patients with CRC [7]. In addition, immune cell infiltration has been shown to be associated with overall survival in CRC. High levels of plasmacytoid dendritic cells have been associated with prolonged overall survival in CRC patients [8]. In addition, innate immune responses, including tumour-associated macrophages, have been shown to be associated with overall survival in patients with CRC patients [9]. A recent study showed that a high macrophage M1: macrophage M2 density ratio was associated with better overall survival in CRC patients [10]. In addition, the tumour-associated macrophage-related biomarkers, including CD68 + CD163+ and CD86+, were potential prognostic markers for CRC patients [11]. Therefore, the characterisation of tumour-infiltrating immune cell levels and the identification of immune cell infiltration-related genes are important for the diagnosis and treatment of CRC, as they will help us to further explore the mechanism of tumour immune infiltration and provide theoretical guidance for the immunotherapeutic response for CRC patients.

In the present study, a comprehensive systematic research was performed to screen immune cell infiltration-associated prognostic genes. An immune cell infiltration-related prognostic model was then constructed to evaluate its performance in predicting of the prognosis of CRC patients.

Methods

Data collection

RNA sequencing data and corresponding clinical information of CRC samples (including 41 normal samples and 454 cancer samples) were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/repository). The TCGA cohort was used as the discovery set. RNA sequencing data and corresponding overall survival information from 585 cancer samples were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/). This cohort (GSE39582) was used as the validation set.

Identification of subtypes associated with immune cell infiltration

First, the TIMER algorithm was used to assess the level of immune cell infiltration in CRC. Then, the association between overall survival of CRC patients and immune cells was evaluated using Kaplan-Meier analysis. The p-value < 0.05 indicated a significant difference. Based on the results of the Kaplan-Meier analysis, the immune cells associated with survival were identified. CRC patients were then divided into two subtypes based on the median level of survival-associated immune cells.

Characterisation of immune status and immune infiltration landscape

The ESTIMATE algorithm was used to evaluate the ESTIMATE score, immune score, stromal score and tumour purity between the two subgroups. Meanwhile, the CIBRSORT algorithm was used to compare the level of immune cell infiltration between the two subgroups.

Enrichment analysis of potential pathways

Gene set variation analysis (GSVA) was used to calculate the score of signalling pathways between the two subgroups [12]. We applied the Gene set enrichment analysis (GSEA) software (v4.1.0) was used to evaluate the potential biological processes that were mainly enriched between the two subgroups [13]. The p-value < 0.05 indicated a significant difference. The results of GSVA and GSEA were visualised using the ClusterProfiler, enrichplot, and ggplot2 packages of R.

Identification of differentially expressed macrophage-related genes (DEMRGs)

Macrophage-related genes (MRGs) between high and low macrophage infiltration groups were identified using the limma package of R. The limma package was also used to identify differentially expressed genes (DEGs) between tumour and normal groups. The | log fold change (FC)| ≥1 and p < 0.05 were set as thresholds to perform the differentially expressed analysis. The differentially expressed macrophage-related genes (DEMRGs) derived from the DEGs and MRGs were then screened using the Venn tool.

Development of a macrophage-related prognostic model

We screened the prognosis-related genes from the DEMRGs using the Univariate Cox analysis. Then, the prognosis-related DEMRGs were used to perform the least absolute shrinkage and selection operator (LASSO) analysis and establish the prognostic model [14,15]. The Kaplan-Meier curve was used to analyse the prognostic differences between the two subgroups. The accuracy of the prognostic model was evaluated through the time-dependent ROC curve, and the “timeROC” package of R was used to perform ROC curve analysis.

Establishment of the nomogram

A nomogram was constructed using the “rms” and “survival” packages of R. In addition, the nomogram predicted survival probability and the observed survival probability were plotted to generate the evaluation calibration curve.

Quantitative real-time polymerase chain reaction (qRT-PCR) analysis

In 2021-2022, a total of eight pairs of cancerous tissues and peritumoral tissues were collected from CRC patients who underwent surgery at the Department of General Surgery, Heyuan People’s Hospital. This research was approved by the Ethics Committee of Heyuan People’s Hospital. Total RNA was extracted from the tissues using TRIzol reagent (Invitrogen, CA, USA). The cDNA was synthesised using reverse transcriptase (Invitrogen, CA, USA). qRT-PCR was performed using a Light Cycler 480 system (Roche). Primers are listed in Table S1.

Statistical analysis

Statistical analysis was performed using R software (version 4.2.1). The Wilcoxon rank-sum test was used to compare the two groups, and the results were visualised using the ‘ggplot2’ package. The ggplot2 package of R was used to assess the correlation between the risk score and the level of immune cell infiltration. A p < 0.05 was reported as statistical difference.

Results

Screening for prognostic immune cells in CRC

We used the TIMER algorithm to quantify the level of immune cell infiltration in CRC samples in the TCGA database. As shown in Figure 1, high macrophage infiltration was associated with poorer prognosis in CRC patients. The CRC patients were then divided into two subtypes according to the median value of the macrophage infiltration level.

Figure 1. Association of B cell (A), dendritic cells (B), macrophage (C), neutrophil (D), T cell CD4 (E), T cell CD8 (F) with overall survival of CRC. CRC patients were divided into two subtypes based on the median level of immune cells.

Comparison of immune status and immune cell infiltration in macrophage-associated subgroups

As shown in Figure 2A–D, ESTIMATE score, immune score and stromal score were significantly decreased in the low macrophage infiltration group, while tumour purity was significantly increased in the low macrophage infiltration group (p < 0.05). In addition, infiltration of T cells CD4 memory resting, T cells CD4 memory activated, NK cells activated, dendritic cells activated, and mast cells activated were significantly increased in the low-macrophage infiltration group (p < 0.05). On the contrary, infiltration of T cells gamma delta, macrophages M0, macrophages M2, and mast cells resting were significantly decreased in the low macrophage infiltration group (p < 0.05) (Figure 2E).

Figure 2. Characterisation of immune landscape between the low- and high-macrophage subgroups. Boxplot shows the level of ESTIMATE score (A), immune score (B), stromal score (C), and tumour purity (D) between the low- and high-macrophage subgroups. Boxplot shows significantly different immunocyte infiltration between the low- and high-macrophage subgroups (E).

Functional enrichment analysis in macrophage-associated subgroups

GSVA was performed to investigate the signalling pathways that differed significantly between the two subgroups. As shown in Figure 3A, a total of 3512 differentially expressed signalling pathways were identified between the two macrophage-associated subgroups. There were 3472 up-regulated pathways and 40 down-regulated pathways. In addition, the autophagy, apoptosis, inflammation, and immune-related pathways were activated in the high macrophage-associated subgroup, such as regulation of autophagy, regulation of autophagosome assembly, autophagosome maturation, apoptotic process, positive regulation of neuronal apoptotic process, macrophage cytokine production, regulation of gamma delta T cell activation, positive regulation of B cell receptor signalling pathway, positive regulation of CD8 positive alpha beta T cell differentiation, positive regulation of macrophage activation, etc (Figure 3B). In addition, GSEA results showed that the regulation of autophagy, chemokine signalling pathway, T cell receptor pathway, B cell receptor pathway, cytokine-cytokine receptor interaction, and natural killer cell mediated cytotoxicity were significantly upregulated in the high macrophage infiltration group (p < 0.05) (Figure 4).

Figure 3. GSVA in macrophage-associated subgroups. (A) The volcano plot depicts the differentially expressed pathways between low- and high-macrophage subgroups. (B) The heatmap plot exhibits the representative pathways between low- and high-macrophage subgroups.

Figure 4. GSEA in the low- and high-macrophage subgroups.

Identification of DEMRGs in CRC

As shown in Figure 5A, a total of 1029 MRGs were screened between the high and low macrophage infiltration groups. Among them, there were 34 upregulated genes and 995 downregulated genes. A total of 7272 DEGs were screened between the tumour and normal groups. Among them, there were 5109 upregulated genes and 2163 downregulated genes (Figure 5B). Finally, a total of 547 DEMRGs were identified using a Venn tool (Figure 5C). The subsequent analysis focused on these DEMRGs.

Figure 5. Identification of DEMRGs in CRC. (A) The volcano plot depicts the MRGs between low- and high-macrophage subgroups. (B) The volcano plot depicts the DEGs between tumour and normal groups. (C) The intersection of genes between the DEGs and MRGs.

Development of a prognostic model based on the DEMRGs

Univariate Cox analysis was performed to identify prognosis-related DEMRGs. As shown in Figure 6A, a total of 25 prognosis-related DEMRGs were screened and used to perform the LASSO COX analysis. Then, seven DEMRGs (RIMKLB, UST, PCOLCE2, ZNF829, TMEM59L, CILP2 and DTNA) were screened and applied to construct the prognostic model (Figure 6B,C). CRC patients were then divided into two risk subgroups based on the median value of the risk scores. Kaplan-Meier analysis showed that the high risk score was associated with a worse prognosis in CRC patients (Figure 6D). In addition, a lower proportion of deaths were observed in the low-risk score group. The expression of RIMKLB, UST, PCOLCE2, TMEM59L, CILP2 and DTNA was downregulated in the low-risk group, whereas the expression of ZNF829 was upregulated in the low-risk group (Figure 6E). The results of time-dependent ROC analysis showed that the AUC of this risk model was 0.63, 0.63 and 0.6 for 1, 3 and 5 years, respectively (Figure 6F). Furthermore, the above results were validated in the GSE39582 dataset (Figure 7A–C).

Figure 6. Establishment of the prognostic model based on macrophage-associated subgroups. (A) Identification of prognosis-associated DEMRGs using the Univariate Cox analysis. (B-C) LASSO analysis identifies seven genes used to construct the prognostic model. (D) Survival analysis of the risk score in the TCGA dataset. (E) Distributions of the risk score, survival status, and the expression level of seven genes in the low- and high-risk groups. (F) ROC curves of the prognostic model.

Figure 7. Validation of the prognostic model. (A) Survival analysis of the risk score in the GSE39582 dataset. (B) Distributions of the risk score, survival status, and the expression level of seven genes in the low- and high-risk groups. (C) ROC curves of the prognostic model.

Risk score was an independent prognostic factor in CRC patients

The independent prognostic performance of the risk model was evaluated by univariate and multivariate COX analyses. Our results showed that the risk score was an independent prognostic factor for CRC patients (Table 1). As shown in Figure 8A, a nomogram based on age, gender and risk score was established to predict the prognosis of CRC patients. Furthermore, the calibration plot indicated that the predictions of the nomogram were consistent with the actual observations (Figure 8B–D).

Figure 8. Development of the nomogram for CRC patients. (A) The overall survival nomogram plot. (B) The calibration plot of the nomogram.

Table 1. Univariate and multivariate COX analyses for risk score and clinicopathological features.

Characteristics	Total(N)	Univariate analysis		Multivariate analysis
		Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
Age	454	1.024 (1.007-1.043)	0.007	1.037 (1.018-1.055)	<0.001
Sex	454
FEMALE	214	Reference
MALE	240	1.138 (0.763-1.697)	0.528
T	454
T1	11	Reference
T2	77	0.455 (0.088-2.355)	0.348	0.318 (0.060-1.675)	0.176
T3	309	1.260 (0.309-5.148)	0.747	0.538 (0.068-4.235)	0.556
T4	56	3.665 (0.856-15.690)	0.080	1.086 (0.132-8.902)	0.939
N	454
N0	267	Reference
N1	105	1.690 (1.017-2.808)	0.043	0.601 (0.201-1.803)	0.364
N2	82	4.017 (2.534-6.369)	<0.001	0.973 (0.340-2.789)	0.960
M	454
M0	333	Reference
MX	50	2.073 (1.121-3.836)	0.020	1.710 (0.879-3.326)	0.114
M1	64	4.457 (2.831-7.017)	<0.001	6.576 (2.173-19.897)	<0.001
Stage	454
Stage II	176	Reference
Stage III	128	1.815 (1.077-3.061)	0.025	2.621 (0.790-8.699)	0.115
Stage IV	64	4.788 (2.825-8.116)	<0.001
Stage I	75	0.569 (0.217-1.488)	0.250	0.933 (0.147-5.919)	0.941
Risk score	454	4.363 (2.334-8.156)	<0.001	3.444 (1.679-7.063)	<0.001

Bolded text indicates that p-values are statistically different.

Association between the risk score and clinicopathological features

No significant differences were observed between the risk score and clinicopathological features, including sex and age (Figure 9A,B). As shown in Figure 9C, CRC patients in the N1 group had a lower risk score than those in the N2 group (p < 0.05), and patients in the N0 group had a lower risk score than those in the N2 group (p < 0.001). CRC patients in the M0 group had a lower risk score than those in the M1 group (p < 0.05) (Figure 9D). CRC patients in the T2 group had a lower risk score than those in the T4 group (p < 0.001) (Figure 9E). CRC patients in the early stage had a lower risk score than those in the advanced stage (p < 0.05) (Figure 9F).

Figure 9. The association between the risk score and clinicopathological features, including sex (A), age (B), N stage (C), M stage (D), T stage (E), and TNM stage (F). An asterisk (*) indicates p < 0.05, two asterisks (**) indicates p < 0.01, three asterisks (***) indicates p < 0.001, ns indicates no significant difference.

Comparison of immune status and immunocyte infiltration between the low and high risk groups

The ESTIMATE score (Figure 10A), immune score (Figure 10B), and stromal score (Figure 10C) were significantly declined in the low risk group, while the tumour purity (Figure 10D) was significantly increased in the low risk group (p < 0.05). In addition, infiltration of T cells CD4 memory resting, T cells CD4 memory activated, and dendritic cells activated was significantly increased in the low risk group (p < 0.05). In contrary, infiltration of T cells CD8, T cells follicular helper, macrophages M0, macrophages M2, and neutrophils was significantly decreased in the low risk group (p < 0.05) (Figure 10E). Furthermore, the risk score was positively correlated with T cells CD8, T cells follicular helper, macrophages M0, macrophages M2, and neutrophils (p < 0.05). On the contrary, the risk score was negatively correlated with T cells CD4 memory resting, T cells CD4 memory activated, NK cells resting, and dendritic cells activated (p < 0.05) (Figure 11).

Figure 10. Characterisation of immune landscape between the low- and high-risk groups. Boxplot shows the level of ESTIMATE score (A), immune score (B), stromal score (C), and tumour purity (D) between the low- and high-risk groups. Boxplot shows significantly different immunocyte infiltration between the low- and high-risk group (E).

Figure 11. The association of risk score with immunocyte infiltration. An asterisk (*) indicates p < 0.05, two asterisks (**) indicates p < 0.01, three asterisks (***) indicates p < 0.001.

Gene expression analysis by qRT-PCR

We collected clinical samples to validate the expression of signature genes in CRC patients. As shown in Figure 12, CILP2 expression was significantly upregulated, but TMEM59L, DTNA, RIMKLB, PCOLCE2, UST and ZNF829 were significantly downregulated in CRC patients compared to normal samples (p < 0.05 or p < 0.01 or p < 0.001). This result is consistent with the findings in the TCGA-COAD dataset (Figure S1).

Figure 12. Validation of signature genes by qRT-PCR analysis. An asterisk (*) indicates p < 0.05, two asterisks (**) indicates p < 0.01, three asterisks (***) indicates p < 0.001.

Discussion

CRC is a malignancy with a poor prognosis and is notably characterised by high mortality and high morbidity [16]. Therefore, there is a need to develop sensitive and reliable prognostic markers to identify patients who may benefit from immunotherapy. Tumour initiation and progression are closely linked to changes in the tumour microenvironment. Tumour cells can shape their microenvironment by producing various chemokines and cytokines [17]. In addition, the tumour microenvironment has been shown to be highly relevant to the prognosis of various cancers [18–20]. Therefore, the tumour microenvironment as a therapeutic target for tumours has attracted the research enthusiasm of scientists. Based on the special role of the tumour microenvironment in tumorigenesis, we identified important genes closely related to the prognosis of CRC patients by screening of prognosis-related immune cells.

The present study was the first to calculate the infiltration of immune cells in CRC samples. And our findings showed that high macrophage infiltration was associated with poor prognosis in CRC patients. Macrophages are innate immune cells, involved in tumour development, metastasis, invasion, and drug resistance [21]. Macrophage-associated therapeutic therapies may complement and synergize with currently available therapies in oncology [22–23]. A recent study showed that high levels of macrophages were associated with a poorer prognosis in lung cancer patients [24]. In addition, tumour-associated macrophages, the active class of inflammatory cells, play an important role in the initiation and progression of CRC [9]. For postoperative adjuvant chemotherapy with stage II colon cancer, tumour-associated macrophages were potential prognostic markers [25]. Consistent with these studies, our findings showed that macrophage infiltration may accelerate the progression of CRC.

Two macrophage-related subtypes were then identified based on the median macrophage infiltration level. Our results showed that the prognosis of the high-macrophage subgroup was worse than that of the low-macrophage subgroup. The high-macrophage subgroup had higher immune, stromal and ESTIMATE scores, but tumour purity was lower. In addition, we found that immune, inflammatory and autophagy-related pathways were activated in the high-macrophage subgroup, indicating that the high-macrophage subset is an immune-friendly tumour. We identified differential genes between the subgroups and established a prognostic risk model by LASSO regression analysis to better apply macrophage-related subgroups in the clinical management of CRC patients. Among these prognostic genes, rimK-like family member B (RIMKLB) is an enzyme involved in mammalian reproductive function [26]. Procollagen C proteinase enhancer 2 (PCOLCE2) is a therapeutic candidate for improving host defence by enhancing of neutrophil oxidative bursts [27]. Transmembrane protein 59-like (TMEM59L) is a brain-specific membrane-anchored protein that plays an important role in the central nervous system function and is involved in autophagy [28,29]. A recent study shows that the TMEM59L gene is associated with prognosis in several tumour types, including bladder urothelial carcinoma, colon adenocarcinoma, and kidney renal clear cell carcinoma [30]. Cartilage intermediate layer protein 2 (CILP2) is involved in neurogenesis through Wnt signalling pathways and ameliorates defects induced by RNF20 depletion [31]. In addition, CILP2 is an independent prognostic biomarker and is associated with advanced stages in CRC patients [32]. Dystrobrevin-α (DTNA) is involved in the maintenance and formation of synapses and thus involving in the regulation of the blood-brain barrier [33]. DTNA regulates TGFβ1 and P53 signalling and activates STATS, thereby contributing to the progession of HBV-evoked hepatocellular carcinoma [34]. In the present study, the risk model including RIMKLB, UST, PCOLCE2, ZNF829, TMEM59L, CILP2, and DTNA genes showed good prognostic value in CRC patients. Thus, this risk model is useful for the diagnosis and treatment of CRC.

Conclusion

Overall, our study identified macrophage-associated subtypes in patients with CRC. In addition, a prognostic model was developed and validated to predict the prognosis of CRC patients. Our findings will help us to better understand the characteristics of macrophage infiltration and provide novel therapeutic options for individualised treatment of CRC patients.

Ethics approval and consent to participate

This study was approved by the Ethics Committee of Heyuan People’s Hospital.

Authors’ contributions

Qi Liu drafted the manuscript and was responsible for the acquisition of data; Li Liao participated in the data analysis and modified the manuscript.

Disclosure statement

All authors declare that they have no competing interests.

Data availability statement

All dataset in the present study is available in TCGA (https://portal.gdc.cancer.gov/) database.

Additional information

Funding

This work did not receive any funding support.