Cellular landscape of primary CRC and liver metastasis

To understand the single-cell level landscape of primary colorectal cancer (CRC) and liver metastases, we reanalyzed the scRNA-seq dataset GSE23155914. This dataset initially comprised 26 samples from 18 patients with colorectal cancer. Following quality filtering of the single-cell sequencing data, the final analysis included 24 samples from 17 patients. These samples consisted of 4 colorectal cancer tissues (Col_T), 3 adjacent normal tissues (Col_N), 9 liver metastasis tissues (Liv_T), and 8 liver normal samples (Liv_N). For data processing, principal component analysis (PCA) was performed with n_comps = 50, followed by the construction of neighborhood graphs with the number of neighbors set to 10. Using the Louvain algorithm, we identified 23 unique cell clusters, which were relatively evenly distributed across all samples (Fig. 1A, Supplementary File Figure S1A), suggesting that there is no significant batch effect interference among the samples. Cell types were annotated using CellTypist based on known marker genes, categorizing 56,163 cells into 9 major cell types (Fig. 1B). T cells and myeloid cells had the highest cell numbers, with 36,566 and 6,549 cells, respectively (Supplementary File Table S1). T cells, NK cells, myeloid cells, mast cells, B cells, plasma cells, epithelial cells, fibroblast cells, and endothelial cells exhibited distinguishing canonical marker gene expression patterns (Fig. 1C, Supplementary File Figure S1B-D). Notably, CD3D, NKG7, S100A9, CPA3, MS4A1, MZB1, KRT18, DCN, and CLDN5 were highly specifically expressed in their respective cell types (Fig. 1D, Supplementary File Figure S1B). These results demonstrate that we have clearly delineated the distinct cell lineages in CRC liver metastasis patient samples.

Fig. 1

Analysis of the colorectal cancer single-cell atlas. (A) Uniform Manifold Approximation and Projection (UMAP) visualization of 56,163 cells categorized into 23 distinct clusters. (B) UMAP delineation of nine principal cell lineages (T cells, NK cells, Myeloid cells, Mast cells, B cells, Plasma cells, Epithelial cells, Fibroblasts, and Endothelial cells) derived from CRC patient samples. (C) Violin plots depicting the distribution of a selected canonical marker gene expression across each of the major cell lineages. (D) UMAP representations of the nine principal cell lineages, each highlighted by the expression of its respective canonical marker gene.

Characterization of different myeloid cell subsets

Immunosuppressive myeloid cells are an indispensable driver in tumor progression and immunotherapy resistance in CRC. Therefore, we further analyzed different myeloid cell subsets to identify key subgroups associated with metastasis. A total of 6,627 myeloid cells were included in the analysis and further classified into 9 subtypes (Fig. 2A-B). These myeloid cells were primarily composed of four cell types: macrophages, monocytes, dendritic cells, and a small number of mast cells. Distinct marker gene expression differences were observed among the four myeloid cell subsets (Fig. 2C-D, Supplementary File Figure S2A-C), and UMAP visualization of marker gene expression also revealed subset-specific expression patterns (Fig. 2E). We assessed the proportion of different myeloid cell subsets across various sample types. Notably, macrophages were significantly enriched in liver metastasis lesions, while they were less prevalent in normal liver tissue, with monocytes showing the opposite pattern (Fig. 2F). This prominent presence of macrophages in liver metastases highlights their immunosuppressive effect, aiding tumor metastasis and colonization. Collectively, characterizing tumor-associated macrophages in CRC may provide valuable insights for understanding tumor metastasis and predicting patient prognosis.

Fig. 2

Gene expression profiling of myeloid cell clusters in CRC patients. (A) A UMAP visualization of myeloid cells delineated into nine distinct clusters. (B) A detailed UMAP plot of 6,549 myeloid cells, categorized by four myeloid cell phenotypes. (C) Violin plots depicting the expression levels of three canonical marker genes for each myeloid cell type (Macrophages: APOE, C1QA, C1QB; Monocytes: S100A9, FCN1, S100A8; Dendritic Cells: IDO1, XCR1, CLEC9A; Mast Cells: CPA3, KIT, MS4A2). (D) A heatmap revealing the expression patterns of three representative genes across macrophages, monocytes, dendritic cells, and mast cells. (E) UMAP plot of the same representative genes highlight their expression within macrophages, monocytes, dendritic cells, and mast cells. (F) Proportional distribution of four myeloid cell types in each sample.

Construction of prognostic signature based on metastasis-specific macrophages

To identify key macrophage subgroups promoting tumor metastasis, we performed sub-cluster analysis on 4,319 macrophage cells. This analysis yielded four macrophage phenotypes, labeled as EEF1G+, C1QC+, SEPT2+, and S100A8 + macrophages (Fig. 3A-B, Supplementary File Figure S3A-B). Notably, EEF1G + macrophages were predominantly found in liver metastases (Fig. 3A), with their marker genes significantly more expressed in liver metastases than in primary lesions (Fig. 3C). C1QC + macrophage marker genes had similar expression levels in both liver metastases and primary lesions. SEPT2 + and S100A8 + macrophages were primarily expressed in primary lesions and reduced in liver metastases (Fig. 3C). Additionally, we compared our four macrophage subtypes with the previously reported seven pan-cancer TAM subtypes15. We found that EEF1G + macrophages predominantly expressed marker genes associated with Inflam-TAMs and IFN-TAMs, while also specifically exhibiting high expression of the Angio-TAMs marker gene VEGFA. Notably, TAM-derived VEGFA can facilitate metastasis by promoting various aspects of tumor progression. In contrast, C1QC + macrophages showed dominant expression of LA-TAMs marker genes, indicating a potential role in lipid metabolism and tissue remodeling. The SEPT2 + and S100A8 + macrophages exhibited moderate expression of marker genes from various oncogenic TAM subtypes (Fig. 3D). Taken together, we consider EEF1G + macrophages to be the primary metastasis-specific macrophage subgroup. We further explored currently available single-cell datasets for colorectal cancer and reanalyzed GSE178318, which comprises paired primary and liver metastasis samples from six patients (Supplementary File Figure S3C and S3D). In this dataset, we observed a significantly higher prevalence of the EEF1G + macrophage subpopulation in liver metastases compared to primary tumors, thereby validating the existence of the liver metastasis-specific macrophage subpopulation identified in this study. Differentially expressed genes (DEGs) analysis for EEF1G + macrophages, compared to other macrophage cell types, revealed 1,264 up-regulated and 107 down-regulated genes (Fig. 3E). We assessed the prognostic value of these DEGs using the expression and survival data of TCGA-COAD cohort. Then, we applied the LASSO-Cox proportional hazards regression model to these prognostic DEGs, ultimately identifying a prognostic signature composed of 8 genes. We defined a risk model based on the EEF1G + macrophage-related gene signature as the EMGS (EEF1G + Macrophage Gene Signature). Subsequently, for each patient, the EMGS was calculated using the formula: EMGS = (0.002 × HIF1A-AS3 expression) + (0.359 × AL645608.7 expression) + (0.022 ×INO80B expression) + (0.154 × HSPA1A expression) + (0.181 × SNAI1 expression) + (0.424 × IFITM10 expression) + (0.248 × ZNF736 expression) + (0.296 × SNHG7 expression) (Fig. 3F).

Fig. 3

Gene expression profiling of macrophage cell cluster in CRC patients. (A) UMAP visualization categorizes macrophage cells into normal colon (Col_N), colon tumor (Col_T), normal liver (Liv_N), and liver tumor (Liv_T) groups. (B) A refined UMAP plot of 4,319 macrophage cells, classified into four macrophage phenotypes (Mph-EEF1G, Mph-C1QC, Mph-SEPT2, and Mph-S100A8). (C) Violin plots illustrate the expression profiles of three key marker genes within the Col_T and Liv_T macrophage populations. (D) The stacked violin plot shows the expression levels of the marker genes, such as IL1B, CCL3, ISG15, CXCL8, ARG1, CD274, MKI67, CDK1, LYVE1, HES1, VEGFA, SPP1, APOC1, and APOE, associated with specific functional states like inflammatory TAMs (Inflam-TAMs), IFN-TAMs, regulatory TAMs (Reg-TAMs), proliferative TAMs (Prolif-TAMs), resident TAMs (RTM-TAMs), angiogenic TAMs (Angio-TAMs), and lipid-associated TAMs (LA-TAMs). Each row represents different macrophage subtypes identified in the study, including Mph-EEF1G, Mph-C1QC, Mph-SEPT2, and Mph-S100A8. (E) A volcano plot highlights differentially expressed genes between the Mph-EEF1G phenotype and other macrophage cell types. (F) Forest plots present hazard ratios (HRs, depicted by pink circles) and 95% confidence intervals (represented by horizontal lines) from Cox regression survival analyses, correlating the overall survival with eight macrophage signature genes.

Performance of the EMGS Signature in training and validation cohorts

After calculating the EMGS for each patient, we used the median EMGS as a cutoff to divide patients into EMGS-high and EMGS-low groups. In the TCGA-COAD training set, the EMGS-high group demonstrated significantly higher mortality risk and lower overall survival rates (Fig. 4A, HR = 2.5, 95% CI: 1.77–3.54, P < 0.0001). Time-dependent ROC curves were employed to estimate the prognostic accuracy of the EMGS model, revealing 1-, 3-, and 5-year AUC values of 0.762, 0.682, and 0.693, respectively (Fig. 4B). Additionally, an overview illustrating risk score, grouping, survival time and status, along with the expression levels of model genes, showed that an increase in risk score was associated with a higher number of patient deaths and an increase in the expression levels of model genes (Fig. 4C). Furthermore, to validate the generalizability of the prognostic model, we calculated EMGS for patients in three independent GEO cohorts. The EMGS-high group had a significantly lower overall survival rate compared to the EMGS-low group in all three datasets: GSE29621 (Fig. 4D, HR: 2.49, 95% CI: 1.29–4.82, P = 0.027), GSE17536 (Fig. 4E, HR: 2.04, 95% CI: 1.38–3, P = 0.0024), and GSE39582 (Fig. 4F, HR: 1.4, 95% CI: 1.1–1.78, P = 0.022). Finally, we performed multivariate Cox proportional hazards analysis using EMGS along with clinical factors such as age, gender, and tumor stage. The results indicated that EMGS is an independent predictor of overall survival (OS) in the TCGA-COAD cohort (P = 0.01, Fig. 4G). In summary, these results demonstrate the superior efficacy of the EMGS risk score as a prognostic signature.

Fig. 4

Survival analysis in high-risk and low-risk groups based on EMGS. (A) Kaplan–Meier curves compared the overall survival of TCGA-COAD patients between high-risk and low-risk groups. (B) ROC curves of the EMGS score for predicting the risk of death at 1, 3, and 5 years. (C) The distribution of risk score (top), survival status (middle), and expression (bottom) of the identified eight EMGS genes. Validation of EMGS in (D) GSE29621 (n = 65) cohort, (E) GSE17536 (n = 177) cohort and (F) GSE39582 (n = 556) cohort. (G) Forest plot shows HRs and 95% confidence intervals (horizontal ranges) derived from multivariate cox regression analyses including EMGS group, age, gender, and stage.

Functional pathways and metabolic alterations associated with the EMGS signature

To elucidate the biological mechanisms for the poorer prognosis associated with a high EMGS score, we initiated a differential expression analysis (DEG) between the high and low-risk groups. We found significant upregulation of 2779 genes and downregulation of 16 genes in the high-risk group (Supplementary File Table S2). Gene Ontology annotation of the DEGs revealed enrichment in biological processes associated with external encapsulating structure organization and extracellular matrix organization, with cellular components and molecular functions closely linked to the cell membrane and extracellular matrix (Fig. 5A). Gene Set Enrichment Analysis (GSEA) indicated that the high-risk group predominantly enriches in pathways related to epithelial-mesenchymal transition (EMT) and angiogenesis, key processes involved in metastasis, while downregulating pathways like the interferon alpha signaling pathway (Supplementary File Fig. 5B-C). Additionally, using gene set variation analysis (GSVA) to calculate sample-wise gene set enrichment scores, we observed significantly enhanced activity of the EMT and TGFβ pathways in the high-risk group compared to the low-risk group (Fig. 5D). These results suggest a greater metastatic potential in the high-risk group. Furthermore, we calculated the metabolic fluxes for 13,082 reactions and 142 metabolic pathway activities using METAFlux (METAbolic Flux balance analysis). We identified 166 fluxes with significant differences between the high and low-risk groups (Supplementary File Table S3), involving processes such as glycolysis, purine metabolism, nucleotide metabolism, and beta oxidation of fatty acids. By integrating fluxes of reactions associated with specific pathways, we found that at the metabolic pathway level, the high-risk group significantly downregulated beta oxidation of fatty acids, folate metabolism, and sulfur metabolism, while upregulating inositol phosphate metabolism and ether lipid metabolism (Fig. 5E, Supplementary File Table S4). These findings highlight the association between the EMGS signature and lipid metabolism reprogramming, reflecting the potential role of targeting metabolic pathways in cancer treatment.

Fig. 5

Functional differences between EMGS high and low groups. (A) Gene Ontology analysis of differentially expressed genes between EMGS high and low groups. (B) GSEA analysis of hallmark pathways. (C) GSEA enrichplot displaying significantly upregulated and downregulated top pathways in the EMGS high group. (D) GSVA activity analysis of hallmark EMT and TGF-β signaling pathway. (E) Differences in metabolic pathways between EMGS high and low groups.

Immune microenvironment shaped by EMGS

To determine the differences in the immune microenvironment between the high and low EMGS score groups, we first deconvoluted bulk expression matrices to infer the proportions of various immune and stromal cell types. CIBERSORT indicated a notably higher proportion of M0 and M2 macrophages compared to the low-risk group, reflecting a positive association between the EMGS score and M2-type macrophages. Conversely, the low-risk group exhibited an elevated percentage of CD8 + T cell (Fig. 6A). MCPcounter also found that the high-risk group was significantly enriched in monocytic lineage, neutrophils, endothelial cells, and fibroblasts (Fig. 6B). Consistently, ESTIMATE indicated that the high-risk group had significantly higher stromal and ESTIMATE scores and lower tumor purity, although there were no significant differences in the immune score between the two groups (Supplementary File Figure S4A). Additionally, we characterized several gene signatures associated with immune evasion through T-cell exclusion. Our analysis revealed that patients in the high-EMGS group had notably elevated Cancer Associated Fibroblasts (CAF) scores (Fig. 6C), T-cell exclusion characteristics (Fig. 6D), and a higher TIDE score (Fig. 6E). These observations were consistent across three other independent validation datasets from GEO (Supplementary File Figure S4B-D). Collectively, these findings suggest that the higher EMGS score correlates with an immunosuppressive tumor microenvironment, which may contribute to reduced responsiveness to immune checkpoint blockade therapies and poorer clinical outcomes.

Fig. 6

Differences in the immune microenvironment between EMGS high and low groups. (A) Proportional distribution of 22 immune cell types obtained through CIBERSORT deconvolution. (B) Proportional distribution of eight immune and two stromal cell populations calculated by MCPcounter. For EMGS high and low groups, scores calculated using the ssGSEA algorithm for (C) CAF score, (D) various immune exclusion signature scores, and (E) TIDE score.

Prediction of drug response in EMGS-high patients

Patients with a high EMGS score, due to their greater metastatic potential and poorer prognosis, urgently require effective drug treatments. Here, we predicted drugs with differential responses between high and low EMGS patients based on drug sensitivity data from colorectal cancer cell lines. The analysis revealed that the EMGS-high group is more sensitive to three anticancer drugs compared to the EMGS-low group: docetaxel (Fig. 7A), a chemotherapy drug used to treat multiple metastatic and non-resectable tumor types, and two targeted cancer drugs, pazopanib and sorafenib (Fig. 7B-C), which are widely used in the treatment of hepatocellular carcinoma and kidney cancer patients, respectively. Additionally, we utilized the DrugMap database (http://drugmap.idrblab.net/) to search for Drug Therapeutic Targets among genes related to EMGS. We discovered that for the gene HSPA1A, which encodes Heat shock protein 70, several drugs are currently in development or clinical trials (Fig. 7D). These drugs may bring clinical benefits for the long-term survival of EMGS-high patients.

Fig. 7

Predicted effective anticancer drugs for the EMGS high group. The IC50 values of three anticancer drugs were predicted by the pRRophetic algorithm and showed significant differences in sensitivity between EMGS high and low groups: (A) Docetaxel, (B) Pazopanib, (C) Sorafenib. (D) Targeted drugs for the EMGS gene HSPA1A identified using the DrugMap database.