Introduction

Colorectal cancer (CRC) is a globally highly prevalent malignancy with a marked propensity for hepatic metastasis1,2. Liver metastasis is the most common form of distant metastasis from CRC due to the close anatomical and vascular connection between the colorectal and the liver3. As many as two-thirds of CRC patients eventually develop colorectal cancer liver metastases (CLM), and liver involvement is the leading cause of death in CRC patients4. Surgical resection is the standard treatment for liver metastases, but it is only feasible in a minority of cases due to various factors such as tumor ___location, size, number of metastases and presence of extrahepatic disease5. Furthermore, even after radical resection, approximately 75% of patients remain at risk of recurrence due to latent spread of tumor cells6. The underlying molecular mechanisms driving CLM occurrence and progression remain poorly understood, posing a significant challenge for effective clinical intervention. However, early detection and timely surgical intervention can significantly improve progression-free survival (PFS) and overall survival (OS) in patients with CLM. Therefore, there is a need to identify new prognostic markers to assess patient prognosis and effectively guide therapeutic decisions.

Immunotherapy has emerged as a promising avenue for the treatment of CLM and has attracted significant attention in the medical field7. By using the patient’s own immune system to target cancer cells, immunotherapy has shown significant efficacy in the treatment of various cancers8. Although PD-1/PD-L1 inhibitors, CTLA-4 inhibitors and CAR-T cell therapy are currently the main modalities of immunotherapy, their application in chronic lymphocytic leukaemia is still in the research phase9. A growing body of evidence highlights the critical role of the tumor immune microenvironment in tumor progression and treatment response10. A comprehensive understanding of the intricate cellular and phenotypic aspects of the tumor immune microenvironment has the potential to facilitate the development of effective cancer therapies, particularly for metastatic disease11. Notably, recent advances in single-cell RNA sequencing have greatly improved our understanding of the tumor immune microenvironment in a wide range of cancers, including hepatocellular carcinoma12, breast cancer13, lung cancer14 and renal cancer15. The integration of single-cell genomics and transcriptomics analyses has further enriched our understanding of genomic heterogeneity within and between tumors.

Despite significant progress in understanding CRC, the complexity and diversity of the tumor immune microenvironment in liver metastases and their primary sites remain incompletely elucidated. Prognostic models based on immune genes are expected to predict survival outcomes in a range of cancers, including colorectal cancer16, oral cancer17, bladder cancer18 and squamous cell carcinoma of the tongue19. However, the specific roles of immune-related genes in the development and progression of CLM need to be further explored and clarified.

Photodynamic Therapy (PDT) is a promising and innovative approach for the treatment of tumors, which uses specific photosensitizers and specific wavelengths of laser light to induce the death of cancer cells by generating reactive oxygen species (ROS)20. For chronic lymphocytic leukaemia, PDT can be used as an adjunctive treatment modality, especially when there is limited residual tumor tissue after surgical resection or other therapeutic interventions. Studies have shown that PDT can induce apoptosis and necrosis of tumor cells, leading to the release of large amounts of antigens and pro-inflammatory factors, which activate the body’s immune system21,22. This immune activation has the potential to enhance the patient’s immune response to the tumor and potentially improve the efficacy of immunotherapy. Conversely, PDT may also induce a degree of immunosuppression, especially at the end of treatment23. This immunosuppressive effect may compromise the efficacy of immunotherapy, and the timing and sequencing of treatments need to be carefully assessed to mitigate any adverse effects. The interaction between PDT and immunotherapy produces effects that are not completely invariant in oncological treatment. Indeed, the exact impact of PDT on the tumor immune microenvironment remains uncertain, emphasising the importance of identifying potential immune genes that may be modulated by PDT. This exploration is crucial for the full integration and application of PDT, elucidation of its mechanisms, development of optimal therapeutic strategies, and provision of a stronger scientific basis for clinical practice.

The integration of weighted gene co-expression network analysis (WGCNA) and differentially expressed gene analyses (DEGs) in The Cancer Genome Atlas (TCGA) database represents a reliable method to identify differentially expressed genes associated with CLM. In this study, we developed a prognostic model for immune-related genes by intersecting these findings with the immune genome and validating the predictive power of the model through least absolute shrinkage selection operator (LASSO) regression analysis. In addition, we used various algorithms to delve into the relationship between riskscores and immune infiltration levels, revealing the intricate interactions between immune gene expression and the tumor microenvironment. We also employed gene ontology (GO) enrichment analysis, DNA methylation analysis, drug sensitivity analysis and gene set variation analysis (GSVA) to elucidate the biological functions and signalling pathways involved in the identified risk prognostic models. Using the Cellchat package, we investigated intercellular communication in high- and low-risk scenarios, providing insights into the dynamics of cell-cell interactions in the tumor microenvironment. Furthermore, by single-cell RNA sequencing (scRNA-seq), we performed a comprehensive analysis of the cellular and molecular heterogeneity of normal colorectal tissue, primary tumors and corresponding liver metastases. Finally, reverse Transcription-Quantitative polymerase chain reaction (qRT-PCR) was used to assess the effect of PDT on the regulation of immune-related genes, while immunohistochemistry (IHC) was used to validate the expression of immune-related genes in normal colorectal tissues, primary tumors and liver metastases. This multifaceted approach holds significant promise for revealing novel biomarkers associated with the diagnosis and prognosis of CLM patients and for unravelling the underlying molecular mechanisms.

Methods and materials

Data acquisition

In this study, colon adenocarcinoma (COAD, a subtype of CRC) tissue chips were obtained from TCGA, including transcriptome data (TCGA-COAD-RNAseq), tissue mutation data (TCGA-COAD-Snp), and methylation data (TCGA-COAD-Meth), along with corresponding clinical information, including 465 tumor samples and 85 adjacent samples. Additionally, transcriptome data (GSE39582) including 562 tumor samples and 17 adjacent samples for colorectal cancer tissue were obtained from GEO. Transcriptome data (GSE41258) including 55 normal colon samples,49 polyp samples, 186 primary colon tumor samples and 47 colon tumor liver Metastasis samples for colorectal cancer tissue were also obtained from GEO. And, the single-cell transcriptome sequencing data of paired samples of adjacent tissue, colorectal cancer, and colorectal cancer liver metastasis from the same patient was obtained from GSE221575.

Selection of key immune genes related to metastasis

This study first utilized the TCGA-COAD-RNAseq data to detect modules of highly correlated genes using the WGCNA package. In the WGCNA analysis process, a soft threshold was selected with a threshold of R² = 0.85. Based on this, hierarchical clustering was used to identify various co-expression modules. Subsequently, these modules were associated with clinical information (Tumor stage, venous invasion, pathological M, pathological N, pathological T, lymphatic invasion) to identify gene modules related to colorectal cancer metastasis through correlation analysis. Subsequently, the corresponding gene sets were intersected with the immune gene set from ImmPort database (https://www.immport.org/resources). Then, survival analysis was employed to filter gene sets with prognostic value with. There is a significant issue of multicollinearity among the genes, and LASSO regression can effectively address this problem, thereby enhancing the reliability of the model. Therefore, after obtaining the aforementioned intersecting genes, LASSO regression analysis was used to remove collinearity based on overall survival time, resulting in the final target gene set.

Prognostic evaluation of key immune genes related to metastasis

After LASSO regression, we developed an immune-related gene prognostic model and calculated riskscores for colorectal cancer patients in the TCGA dataset based on the final target gene set. Then, Samples were then divided into high-risk and low-risk groups based on the mean riskscore of all samples. Survival curves were used to evaluate differences in overall survival between the groups, and a corresponding clinical predictive nomogram was constructed using the RMS package. The predictive performance of the model was evaluated by the AUC area of the ROC curve. Additionally, colorectal cancer tissue RNAseq from GEO (GSE39582) were used as an external validation set to demonstrate the universality and reproducibility of the gene set in predicting overall survival prognosis. In addition to overall survival, differences in disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) between the high-risk and low-risk groups were analyzed. This study also assessed the value of the model in predicting DSS, DFI, and PFI.

Prognostic model and drug sensitivity analysis

Half maximal inhibitory concentration (IC50) values for 481 small molecules in 1001 cell lines, along with their corresponding mRNA gene expression data from the prognostic model’s gene set, were obtained from the Genomics of Therapeutics Response Portal (CTRP). The mRNA expression data was merged with drug sensitivity information and subjected to Pearson correlation analysis to evaluate the association between gene mRNA expression and drug IC50 values. Subsequently, the p-values were adjusted using the False Discovery Rate (FDR). Subsequent to this, pertinent data from colon cancer cell lines were collected to compute riskscore for each cell line using the prognostic model. The correlation between the riskscore of colon cancer cell lines and drug IC50 values was examined. Furthermore, based on the mean riskscore, the colon cancer cell lines were categorized into high-risk and low-risk groups, and the differences in drug IC50 between the two groups were assessed using violin plots.

Analysis of genetic background differences in different risk groups

Based on the patient ID, mutation data samples (TCGA-COAD-Snp) and methylation data samples (TCGA-COAD-Meth) were obtain from TCGA, matching the patients of transcriptome data. Then, mutation data samples (TCGA-COAD-Snp) and methylation data samples (TCGA-COAD-Meth) obtained the corresponding grouping information for their transcriptome based on sample IDs. The mutation data and methylation data from TCGA were divided into high-risk and low-risk groups for sample analysis. The differences in mutation burden and methylation between the different groups were then analyzed. Maftool package was used for visual analysis of the mutation data to compare and analyze gene mutation differences and the affected relevant pathways in different groups. Using the ELMER package to analyze methylation data, we obtained information on the differences in methylation between different groups.

Analysis of the immune microenvironment based on key immune genes

Utilizing Cibersoft package to analyze the immune cell infiltration in colorectal cancer from TCGA based on transcriptome data, and using SingR package for cell annotation on single-cell datasets, thereby combining colorectal cancer tissue and single-cell data to analyze the immune infiltrating microenvironment. Scoring the single-cell data based on the key immune gene set related to metastasis, dividing the single cells into high-risk and low-risk groups using the median as a threshold, and subsequently using Cellchat package for cell interaction analysis to identify differential cellular pathways between the two groups. Employing the Monocle package for pseudo-time series analysis to compare the differences in cell developmental trajectories between the two groups. Evaluating the distribution characteristics of active oxygen in the immune microenvironment of colorectal cancer using ROS-generated gene sets (Table S1), and analyzing the relationship between active oxygen and the key immune gene set related to metastasis through co-localization analysis.

Photodynamic therapy in HCT15 and HCT116 cell line

HCT15 and HCT116 cells were seeded in medium in 6-well plates and incubated overnight, and after washing twice with PBS, they were incubated with 0.5 mM sodium Hua porphyrin (DVDMS) in 2 mL medium containing 10% fetal bovine serum (FBS). After the addition of photosensitizer and incubation for 6 h, the 6-well plate was irradiated with a 630 nm wavelength laser until the total energy density reached 10 J/cm2.

Study ethics declarations

Animal care and interventions were performed in accordance with the National Institute of Health Guide for the Care and Use of Experimental Animals, and the animal protocol was approved by the Animal Ethics Committee of the Second Xiangya Hospital of Central South University. The study was conducted in compliance with ARRIVE guidelines.

Mouse subcutaneous graft tumor model construction and treatment

Eight 6 ~ 8-week-old male C57BL/6J mice (purchased from Hunan SJA Laboratory Animal Co., Ltd, Changsha, China) were taken, and 5 × 105 MC38 cells were resuspended with 200 µL of PBS and inoculated into the subcutaneous tissues on the right side of the mice. When the tumor volume reached 80 ~ 100 mm3, the mice were divided into 2 groups: the control group and the PDT group (15 mg/kg mouse body weight, light dose of 100 J/cm2, laser wavelength of 660 nm). The photosensitiser was injected into the tail vein on days 0, 3, 6 and 9, respectively, and the tumors in the PDT group were irradiated 8 ~ 10 h after the photosensitiser injection, while the NC group was not subjected to additional treatment. After anesthesia with 1.5% isoflurane, mice were executed by cervical dislocation and subcutaneous tumors were obtained. tumor volume was calculated as:0.52*L*W2.

Reverse transcription-quantitative PCR

Total RNA was isolated from the HCT15 cell line and HCT116 cell line using TRIzol reagent (Invitrogen) and then reverse transcribed using the PrimeScript RT-PCR kit (Takara, Japan) according to the manufacturer’s instructions. qPCR was used to detect the relative expression of mRNA. GAPDH was used as an internal control. The relative expression of relevant regulatory genes was normalized using the expression of GAPDH as a standard and calculated using the 2-∆∆Ct method.

Immunohistochemistry

Tumors were fixed in 10% neutral buffered paraformaldehyde. Then, the samples were embedded in paraffin. The numbered slices were placed on a roasting machine and baked at 65 °C for 45 min. Sections were sequentially immersed in xylene I for 10 min, xylene II for 8 min, xylene III for 8 min, absolute ethanol for 5 min, 90% ethanol for 3 min, 80% ethanol for 3 min, and 70% ethanol for 3 min. They were washed 3 times with PBS shaking for 3 min each time, followed by 3 times with distilled water shaking. The antigen repair solution was then performed. The slides were then submerged in PBS for five minutes, repeated twice, and submerged in distilled water for three minutes, repeated twice. Endogenous catalase was removed by immersion in 3% H2O2 for 10 min. The tissue portion of the sections was blocked by adding 10% goat serum blocking solution for 1 h at room temperature. The blocking solution was shaken off, and the tissue was covered by dropping primary antibody, then the wet box was incubated overnight in a cold room at 4 °C. The next day, the wet box was removed and rewarmed for 30 min at room temperature. The cells were washed 3 times with PBS shaking for 3 min each time, followed by 3 times with distilled water shaking for 3 min each time. The water was dried by shaking, the beads were blotted, secondary antibodies were added dropfold and incubated for 10 min at room temperature. The cells were washed 3 times with PBS shaking for 3 min each time, followed by 3 times with distilled water shaking for 3 min each time. Each slide was stained by dropping an appropriate amount of DAB solution. Hematoxylin staining was followed for 1 min. Rinse with tap water for 3 min. The sections were further immersed in 70% alcohol for 1 min, 80% alcohol for 1 min, 95% alcohol for 2 min, absolute ethanol for 4 min, xylene I for 3 min, and xylene II for 3 min. Cool and dry the slices, then add appropriate amount of neutral gum to seal the slices, cover the cover glass, and be careful of bubbles. After drying completely, the specimens were observed under a microscope and photographed.

Statistical analysis

Continuous variables between two groups were compared using the Mann-Whitney U test; continuous variables among three groups were compared using the Kruskal-Wallis test, followed by pairwise comparisons between groups. Survival analysis with Kaplan-Meier (K-M) curves and log-rank tests were conducted to assess the impact of various factors on patient survival. Assess the impact of a single factor on survival time using the univariate Cox proportional hazards model, and calculate the hazard ratio for that factor. Time-dependent receiver operating characteristic (ROC) curves were utilized to evaluate the prognostic value of the prognostic model.

Results

Identification of metastasis-associated hub modules in CRC and identification of immune prognostic genes

The overall design and flowchart of this study are shown in Fig. 1. WGCNA analysis was performed using TCGA-COAD RNA-seq data to establish gene co-expression network (Fig. 2A). Correlation analysis was performed between gene modules and clinical information, and modules related to metastasis were selected for further analysis. Based on the P value (P < = 0.05), grey modules were identified as metastasis-related modules (Fig. 2B), containing a total of 1765 genes (Table S2). We took the intersection of these genes with the immune gene set (immune dataset from immport: https://www.immport.org/shared/home), resulting in 199 genes (Fig. 2C). We then performed survival analyses of these genes using survival curves and COX regression to screen for genes of prognostic value and ultimately identified 11 immune genes associated with metastasis (Table S3). Subsequently, GSVA scores for each sample were calculated based on these 11 genes, and the results showed that the GSVA scores of colorectal cancer tumor samples were significantly higher than those of normal samples (Fig. 2F). Additionally, diagnostic analysis based on the GSVA scores revealed that the area under the ROC curve (AUC) was 0.77 (Fig. 2G). Therefore, these 11 genes effectively distinguished the tumor samples from adjacent non-cancer samples, suggesting that they are characteristic genes for colorectal cancer. The results of univariate COX regression indicates that these 11 genes can be used as candidate genes for prognostic analysis of colorectal cancer (Fig. 2H).

Fig. 1
figure 1

Flow chart of the study design.

Full size image
Fig. 2
figure 2

Construction of gene co-expression modules related to clinical features of colorectal cancer. (A) Cluster analysis of gene modules. (B) Correlation analysis between gene modules and tumor clinical information. Grey modules were identified as transfer-related. (C) Intersection of metastasis-related grey module genes and immune-related genes. (D) Principal component analysis of colorectal cancer samples based on prognostic genes. (E) Heatmap of 11 genes (TCGA). (F) GSVA scores of 11 genes in colorectal cancer tumors and normal colon groups (TCGA). (G) ROC curve for colorectal cancer diagnosis based on GSVA scores of 11 genes (TCGA). (H) Forest plot of univariate COX regression for prognostic genes.

Full size image

Construction and validation of immune-related gene prognostic (IRGP) model

We performed LASSO regression analysis on the 11 genes to generate an immune-related gene prognostic model, including 7 genes (HSPA1A, ULBP2, RBP7, OXT, SLC11A1, INHBB, and ICOS) (Fig. 3A, B). The riskscore calculation formula for (IRGP) model is riskscore = 0.0268305149779617 * SLC11A1 + 0.0989521440194191 * OXT + 0.0975401267404972 * ULBP2 + 0.0836439907393528 * RBP7 + 0.0189408929664606 * INHBB − 0.136824461009379 * ICOS + 0.146741890849542 * HSPA1A. Based on the median riskscore, colorectal cancer patients were categorised into low-risk and high-risk groups. We validated the accuracy of the IRGP model in the TCGA (Fig. 3C, D) and GEO (Fig. 3E, F) databases using K-M survival analysis and receiver operating characteristic curves (ROC). The results showed that the survival of CRC patients in the low-risk group was better than that of the high-risk group in both the TCGA cohort (p = 2.6e-5, Fig. 3C) and the GEO cohort (p = 4.7e-5, Fig. 3E). In addition, ROC showed that the IRGP model had high accuracy in predicting OS in the TCGA cohort (Fig. 3D) and the GEO cohort (Fig. 3F): in the TCGA cohort, the AUC values for OS at 1, 3, and 5 years were 0.70, 0.68, and 0.58, respectively (Fig. 3D). The 1-, 3-, 5-, 7-, and 10-year AUC values for OS were 0.62, 0.62, 0.63, 0.59 and 0.66, respectively (Fig. 3F). Based on the expression of HSPA1A, ULBP2, RBP7, OXT, SLC11A1, INHBB and ICOS, we constructed a nomogram using the TCGA-COAD dataset. By scoring the above genes, we found that higher total scores were associated with lower 1-, 3- and 5-year survival probabilities (Fig. 3G).

Fig. 3
figure 3

Construction and verification of the immune prognosis model. (A) Illustration for LASSO coefficient spectrum of overlapping gene. (B) Adjusted parameters of LASSO regression model. (C) K-M curve analysis of OS in high-risk and low-risk groups base on TCGA databases. (D) ROC curve of riskscore in the TCGA databases. (E) K-M curve analysis of OS in high-risk and low-risk groups base on GEO databases (GSE39582). (F) ROC curve of riskscore in the GEO databases (GSE39582). (G) Nomogram of the prognostic model.

Full size image

To validate the accuracy of the model, we next assessed the disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) of CRC patients in the TCGA database using K-M analysis. We found that disease-free survival (PFI), disease-free progression interval (DFI), and disease-free survival (DSS) were significantly longer in the low-risk group compared with the high-risk group (Supplementary Figs. 1 A, 1B, and 1 C). Nomograms based on the IRGP model successfully predicted DSS, DFI, and DFI in CRC patients (Supplementary Figs. 1D, 1E, and 1 F). The AUC values of 1-year, 3-year, and 5-year PFI were 0.52, 0.62, and 0.69, respectively (Supplementary Fig. 1G). The AUC values of 1-year, 3-year, and 5-year DFI were 0.58, 0.70, and 0.83, respectively (Supplementary Fig. 1H)0.1-year, 3-year and 5-year AUC values for DSS were 0.55, 0.66 and 0.66, respectively (Supplementary Fig. 1I). These results suggest that the prognostic risk model has a strong ability to predict the survival of CRC patients.

Then, we utilized 14 lists of tumor progression-related genes obtained from CancerSEA (http://biocc.hrbmu.edu.cn/CancerSEA) which encompass genes associated with angiogenesis, apoptosis, cell cycle, differentiation, DNA damage, DNA repair, EMT, hypoxia, inflammation, invasion, metastasis, proliferation, quiescence, and stemness for Gene Set Variation Analysis (GSVA) scoring. The results demonstrate that the riskscore effectively captures these key biological processes. Combined with GSEA analysis, it was found that biological processes such as colorectal cancer metastasis and angiogenesis are significantly enhanced in the high-risk group (Supplementary Fig. 2). These results suggest that patients with high-riskscores have a higher degree of tumor malignancy. Consistently, through the analysis of TCGA clinical data, we found that patients with higher riskscores had higher stages of colorectal cancer tumors, with more evident tumor invasion and metastasis (Supplementary Fig. 3). These results to some extent explain why colorectal cancer patients with high-riskscores have a poorer prognosis.

Tumor genetic landscapes of different risk subgroups in the IRGP model

To enhance the clinical application of the IRGP model, we further investigated the expression and corresponding functions of differential genes in different risk groups in the model. The transcriptome data matched 375 corresponding mutation data samples, with 172 cases in the high-risk group and 203 cases in the low-risk group. We analysed the functions of the genes in the low-risk group (Fig. 4A) and the high-risk group (Fig. 4B), respectively, and counted the functions in which these genes were involved or participated. In order to clarify the differences between the high-risk and low-risk groups, we also performed GO pathway analysis on genes that were differentially expressed between the two groups. The results showed that these immune-related prognostic genes were mainly involved in the PTK-RAS pathway, NOTCH signalling pathway, Wnt signalling pathway and TP53 (Fig. 4C and D). Figures 4E and F show the co-occurring mutations and mutually exclusive mutations studied and observed in the low-risk and high-risk groups. In the low-risk group, OBSCN mutations were almost co-mutated with TTN, PIKA3A, SYNE1, FAT4 and MUC16 (P < 0.05), and MUC16 was co-mutated with TTN, SYNE1 and FAT4 (P < 0.05). In addition, in the low-risk group, APC and OBSCN mutations were almost mutually exclusive (P < 0.05). In the high-risk group, TP53 mutations were almost mutually exclusive with OBSCN, MUC16, and SYNE1 mutations (P < 0.05). In the high-risk group, FAT4 mutations were almost co-mutated with APC, TTN, PIKA3A, SYNE1, MUC16 and MUC16 (P < 0.05). These interrelationships and expression differences may reveal the reasons why CRC patients move towards different regression outcomes.

Fig. 4
figure 4

Genetic landscape of tumors with different risk groupings in the IRGP model (high-risk group: 180 cases, low-risk group:107 cases). (A) Drugs affected by mutated genes in the low-risk group. (B) Drugs affected by mutated genes in the high-risk group. (C) GO enrichment analyses of low-risk groups. (D) GO enrichment analyses of high-risk groups. (E) The heatmap illustrates the co-occurrence and mutually exclusive mutations in low-risk groups. (F) The heatmap illustrates the co-occurrence and mutually exclusive mutations in high-risk groups. (G,H) Comparison of TMB between patients in the high- and low-riskscore groups. (I) The mutational landscape of low-risk groups and high-risk groups.

Full size image

To clarify whether patients in different risk subgroups were associated with oncogene mutations, we compared the tumor mutational burden (TMB, here we identified TMB as the total number of somatic non-synonymous mutations detected per megabase pair (Mb) in patients with high and low-riskscores. Unsurprisingly, the median TMB was higher in patients in the high-risk scoring group relative to the low-risk group (1.59/Mb vs. 1.44/Mb, Fig. 4G, H). Meanwhile, the mutation rates of some classical oncogenes differed between the two groups, and we listed the mutation rates and mutation types of the ten genes with the highest mutation frequencies (Fig. 4I), such as the frequency of TP53 mutation was 55% in the high-risk group and 40% in the low-risk group.

To further reveal differences in the tumor genetic landscape among CRC patients in different risk subgroups, we also analysed gene expression and methylation in both groups. The transcriptome data matched 287 corresponding methylation data samples, with 180 cases in the high-risk group and 107 cases in the low-risk group. We found that the genomes of patients in the high-risk group were more methylated and distal genes were less expressed relative to the low-risk group (Supplementary Fig. 4A). There are 95,536 differential methylated sites, which cover 920 immune-related genes, including six target genes: RBP7, INHBB, ICOS, SLC11A1, HSPA1A, and OXT (Table S4). Also, promoter methylation was negitively correlated with gene expression scores (Supplementary Fig. 4B). Next, we functionally enriched the intersection of differentially methylated and differentially expressed genes and found that the enriched pathways contained lots of about immune-related pathways (Supplementary Fig. 4C). The intersection genes include 197 immune-related genes, covering four target genes: SLC11A1, RBP7, INHBB, and HSPA1A ( Table S5). Together, these data reveal differences in the intrinsic genetic landscapes of CRC patients in different risk subgroups.

The IRGP model predicts the drug treatment response in colorectal cancer patients

In order to explain the reason for the model’s ability to predict patient prognosis, an attempt was made to analyse the relationship between the model and chemotherapeutic drugs. Based on the drug sensitivity analysis, seven genes - HSPA1A, ULBP2, RBP7, OXT, SLC11A1, INHBB, and ICOS were found to broadly influence drug sensitivity in tumor treatment (Fig. 5A). The top 30 drugs show a broad correlation with these six genes (Table S6). This result suggests that the IRGP model can predict the drug treatment response in colorectal cancer. Inclusion of 46 colorectal cancer cell lines, the IRGP model riskscores were calculated using the transcriptome data of each cell line. Subsequently, based on the mean riskscores, the cell lines were grouped, with 19 colorectal cancer cell lines in the high-risk group and 27 in the low-risk group. The IC50 values of Oxaliplatin, Cisplatin, Irinotecan, and 5-Fluorouracil - four chemotherapeutic drugs, were obtained for these 46 cell lines ( Table S7). The IC50 of Oxaliplatin in colorectal cancer cell lines in the high-risk group was higher than that in the low-risk group (Fig. 5B). The IC50 of Irinotecan in colorectal cancer cell lines in the high-risk group was higher than that in the low-risk group (Fig. 5C). The IC50 of 5-Fluorouracil in colorectal cancer cell lines in the high-risk group was higher than that in the low-risk group (Fig. 5D). The IC50 of Cisplatin in colorectal cancer cell lines in the high-risk group was higher than that in the low-risk group (Fig. 5E). Consistently, clinical response analysis based on TCGA clinical information indicated that patients with lower riskscores had better clinical response to drug treatments (Fig. 5F). These results suggest that the riskscores of the IRGP model can effectively predict the clinical response to chemotherapy in patients.

Fig. 5
figure 5

IRGP model and drug sensitivity. (A) Drug sensitivity analysis of IRGP model genes. (B) Relationship between riskscores of colorectal cancer cell lines based on the IRGP model and Oxaliplatin IC50. (C) Relationship between riskscores of colorectal cancer cell lines based on the IRGP model and Irinotecan IC50. (D) Relationship between riskscores of colorectal cancer cell lines based on the IRGP model and 5-Fluorouracil IC50. (E) Relationship between riskscores of colorectal cancer cell lines based on the IRGP model and Cisplatin IC50. (F) Relationship between riskscores of colorectal cancer patients and clinical response to drug treatment.

Full size image

Relationship between the IRGP model riskscore and immune status

To further gauge the influence of other factors in the tumor microenvironment between the different subgroups, we further analysed the relationship between the model and the infiltration of 21 immune cells. Analysis of immune infiltration using cibersort package showed differences in the abundance of immune cell infiltration in the CRC high-risk group, the low-risk group, and in the surrounding normal colorectalic tissue. The infiltration of CD4+ T cells, plasma cells, neutrophils, CD8+ T cells, and memory-activated CD4+ T cells was significantly higher in the low-risk group than in the high-risk group (Fig. 6A, B). In the high-risk group, only M0 macrophages had a higher level of infiltration (Fig. 6A, B). In the normal colorectal tissue group, the infiltration of B cells, CD4+ T cells, plasma cells, memory-resting activated CD4+ T cells, M2 macrophages and activated mast cells was significantly higher than that in the high-risk group. Infiltration of M0 macrophages and resting mast cells was significantly lower in the normal paraneoplastic group than in the high-risk group. However, memory B cells and resting status myeloid dendritic cells were not significantly correlated with the riskscore.

Fig. 6
figure 6

Immune status in different subgroups ( normal group:41 cases, high-risk group: 217 cases, low-risk group:214 cases). (A,B) Immune cell infiltration in different subgroups. (C) TIDE score between the three groups.

Full size image

The TIDE (Tumor Immune Dysfunction and Exclusion) score is a biomarker tool based on transcriptome sequencing data that assesses the tumor’s immune evasion capability and predicts the response to immune checkpoint inhibitors (ICIs)24,25. By utilizing contact indicators such as MDSC, CAF, and TAM M2, one can obtain the Immune Dysfunction Score, Immune Exclusion Score, and the combined score, known as the TIDE score, which includes both the Immune Dysfunction Score and the Immune Exclusion Score. The TIDE score can be used to predict responses to immunotherapy, inform personalized treatment decisions, and investigate mechanisms of immune evasion. We next evaluated the potential clinical efficacy of immunotherapy in different groups by TIDE method. Normal colorectalic tissues had higher TIDE scores than the high- and low-risk groups (Fig. 6C). The high-risk group had a higher exclusion score than the low-risk group, but did not differ from the low-risk group in terms of T-cell dysfunction (Fig. 6C). These differences highlight microenvironmental differences between colorectal cancer and normal tissue adjacent to the cancer, and between different risk subgroups within the cancer, which may contribute to the different outcomes.

Inferring intercellular communication in the tumor immune microenvironment

The tumor immune microenvironment contains different types of cells and its complex intercellular communication is closely related to tumor metastasis26. Therefore, we inferred the interactions between different cell populations in the high-risk and low-risk groups based on CellChat (Fig. 7). Collagen signaling plays an important role in tumor metastasis. Previous studies have indicated that the collaboration between macrophages, fibroblasts, and collagen in promoting tumor function is primarily attributed to their facilitation of immune evasion and metastasis in pancreatic cancer27,28,29. The results above show that compared to the low-risk group, the high-risk group exhibited significant macrophage infiltration; therefore, we first considered the possibility of fibroblast-macrophage interactions in collagen signaling. Unsurprisingly, in the high-risk group, most of the fibroblast interactions occurred between epithelial cells and macrophages, whereas in the low-risk group, collagen signalling between fibroblasts and epithelial cells was markedly reduced (Fig. 7A, B). This phenomenon suggests that collagen released by fibroblasts may interact with macrophage-associated receptors to promote extracellular matrix (ECM) formation and ultimately tumor growth. This prompted us to further focus on IL-1, a regulator of the collagen signalling pathway, and we found that in the IL1 signalling pathway, the main interacting cells in the high-risk group were epithelial cells and macrophages, whereas in the low-risk group, it changed to fibroblasts and macrophages (Fig. 7C, D). This shift may have affected the secretion of IL-1 by macrophages, which led to tumor regression.

Fig. 7
figure 7

Inference and analysis of cell-cell communication using CellChat. (A) The Cell–cell communications between different cell types regarding the collagen signaling pathway network. (B) The heatmap shows the communication probability of the collagen signaling pathway network. (C) The Cell–cell communications between different cell types regarding the IL1 signaling pathway network. (D) The heatmap shows the communication probability of the IL1 signaling pathway network. (E) The Cell–cell communications between different cell types regarding the MHC-I signaling pathway network. (F) The heatmap shows the communication probability of the MHC-I signaling pathway network. (G) The Cell–cell communications between different cell types regarding the MHC-II signaling pathway network. (H) The heatmap shows the communication probability of the MHC-II signaling pathway network.

Full size image

Tumor cells can evade major histocompatibility complex (MHC)-mediated immune surveillance through a variety of mechanisms, such as down-regulating the expression of MHC molecules, or altering their mechanisms of antigen processing and presentation30. Therefore, we also analysed cellular interactions between MHC I/II signalling pathways. We found that signalling interactions between epithelial cells and macrophages were attenuated in the low-risk group relative to the high-risk group, but almost all other intercellular signalling interactions were enhanced, especially between epithelial cells and T and B cells, as well as tissue stem cells (Fig. 7E-H).

To further reveal the impact of signalling pathways of intercellular communication on tumors, we quantified the strength of output interactions of each signalling pathway in the high- and low-risk groups. We found that the output interaction strengths of IL1 and IL10 pathways were higher in the high-risk group than in the low-risk group, whereas the output interaction strengths of CD99 were significantly higher in the low-risk group than in the high-risk group (Supplementary Fig. 5A, B). In addition, the expression levels of immune checkpoint genes (CD45, IL1, IL4, IL6, IL2, CD40, CD22, and CD23) were also higher in macrophages in the high-risk group than in the low-risk group. Among them, macrophages were the main cell type expressing ligands and receptors that are actively involved in cell-cell interactions. These results suggest that macrophages play a central role in the tumor immune microenvironment. The riskscore of B cells, endothelial cells, epithelial cells, fibroblasts, macrophages and T cells were statistically significant (P < 0.01) in normal colorectal tissues, primary tumors and liver metastases (Supplementary Fig. 5C).

Single-cell landscape analysis of heterogeneous tumor microenvironments in colorectal normal tissue, colorectal cancer and CLM

We used UMAP to downscale the dataset of primary colorectal cancer, colorectal normal tissue and CLM cell populations in CLM (GSE221575) from the GEO database to illustrate the distribution of single-cell sequencing profiles (Fig. 8). Based on CellMarker expression patterns, we obtained seven cell types (macrophages, fibroblasts, B cells, T cells, epithelial cells, tissue stem cells and endothelial cells). Following our IRGP model, we first divided colorectal cancer tissues into high-risk and low-risk groups. We found that the high-risk group had more macrophages and fibroblasts. Epithelial cells, T cells and B cells were more abundant in the low-risk group (Fig. 8A). Compared to the low-risk group, the high-risk group has a higher proportion of epithelial cells, while the proportion of other cells in the corresponding microenvironment, including fibroblasts, epithelial cells, and tissue stem cells, increases. In contrast, the proportions of immune cells such as B cells and T cells decrease (Fig. 8B). These factors all contribute to the metastasis of tumor cells and immune evasion.

Fig. 8
figure 8

Single-cell sequencing analysis. (A) UMAP showing seven cell types (macrophages, fibroblasts, B cells, T cells, epithelial cells, tissue stem cells, endothelial cells) based on CellMarker expression. High-risk groups have more macrophages and fibroblasts, while low-risk groups show higher levels of T cells, B cells, and epithelial cells. (B) Cell-type proportion comparison between high-risk and low-risk groups. High-risk groups exhibit increased epithelial cells and fibroblasts, with decreased T and B cells. (C) Cell composition across normal tissue, primary tumors, and liver metastases. Macrophages decrease progressively, fibroblasts are highest in primary tumors, and B cells are lowest in tumors and metastases. (D) Cell-type proportion across normal tissue, primary tumors, and liver metastases. (E,F) Riskscore increases progressively from normal tissue to primary tumors and liver metastases, with macrophage riskscores following a similar trend.

Full size image

We further analysed the monocyte landscape in different tissue types. We found a progressive decrease in the number of macrophages in liver metastases, primary tumors and normal colorectal tissues (Fig. 8C). The proportion of B cells and T cells were significantly lower in primary tumors compared with normal colorectal tissues and liver metastases and (Fig. 8D). In addition, the proportion of fibroblasts and macrophages was higher in the primary tumor and liver metastases than in normal colorectal tissues, and it was higher in the liver metastases than in the primary tumor (Fig. 8D). The riskscore shows an increasing trend in the colon, primary colorectal cancer lesions, and colorectal cancer liver metastases, with statistically significant differences (Fig. 8E). The results of riskscore showed a trend of increasing riskscore of macrophages in liver metastases, primary lesions and normal colorectal cancer tissues (Fig. 8F).

To further investigate the process of cell lineage differentiation in liver metastases, we used RNA velocity analysis combined with PAGA trajectory analysis and quasi-temporal trajectory analysis to describe the differentiation pathways of tumor cells in detail. Notably, during the determination of developmental trajectories, liver metastasis cells were found to be at a later stage of differentiation than primary tumor cells (Supplementary Fig. 6A). The differentiation status of cells in the high-risk group was also later compared to cells in the low-risk group (Supplementary Fig. 6B).

The relationship between risk score and colorectal cancer liver metastasis

The analysis of the colonic epithelium through single-cell sequencing found that the riskscores in the epithelium of primary colorectal cancer and liver metastatic lesions of colorectal cancer are higher than those in normal colonic epithelium, and the riskscore in the epithelial cells of liver metastatic lesions is higher than that in the primary colorectal cancer epithelium (Fig. 9A). Subsequently, the analysis of the GSE41258 dataset revealed no significant differences in riskscores between normal colonic tissue and colonic polyps, while the riskscores in colorectal cancer tissues and liver metastatic tissues of colorectal cancer significantly increased, with the riskscore in liver metastatic tissues being higher than that in colorectal cancer tissues (Fig. 9B). These results indicate that the riskscore is closely related to liver metastasis of colorectal cancer.

Fig. 9
figure 9

PDT inhibits the metastasis of colorectal cancer by regulating the riskscore. (A) Riskscores in the epithelium of normal colon, colorectal cancer primary lesions, and liver metastatic lesions from single-cell sequencing (GSE221575). (B) Riskscores in tissues from normal colon, colorectal cancer primary lesions, and liver metastatic lesions (GSE41258). (C) KEGG pathway information related to colorectal epithelial riskscore. (D) GSVA of ROS gene set in colorectal cancer tumor tissues compared to normal tissues (TCGA). (E) Correlation between colorectal cancer riskscore and GSVA of ROS gene set (TCGA). (F) Relationship between GSVA of ROS gene set and prognosis of colorectal cancer (TCGA). (G) Results of ROS flow cytometry before and after PDT treatment in HCT15 cells (blue: ROS in HCT15 without PDT treatment, red: ROS in HCT15 after PDT treatment). (H) Changes in seven immune genes associated with riskscore in the HCT15 cell line before and after PDT intervention measured by -∆T. (I) Changes in riskscore based on ∆T in the HCT15 cell line before and after PDT intervention (∆T = T of the target gene - T of GADPH). (J) Changes in seven immune genes associated with riskscore in the HCT116 cell line before and after PDT intervention. (K) Changes in riskscore -∆T in the HCT116 cell line before and after PDT intervention (∆T = T of the target gene - T of GADPH).

Full size image

The functional enrichment results of differentially expressed genes between high-risk and low-risk groups of epithelial cells showed that the KEGG analysis of differential genes based on high-risk and low-risk scores in colorectal cancer epithelial tissues demonstrated that, apart from affecting lymphocyte migration, the results are consistent with single-cell immune infiltration findings. Furthermore, the riskscore is closely associated with necroptosis, focal adhesion, tight junctions, apoptosis, and the HIF-1 signaling pathway (Fig. 9C). Focal adhesion and tight junctions are closely related to the metastasis of tumor cells, while necroptosis and apoptosis are closely related to the survival of tumor cells. These results suggest that the riskscore is associated with liver metastasis in colorectal cancer.

The relationship between risk score and colorectal cancer liver metastasis

The HIF-1 signaling pathway is closely related to ROS; therefore, PDT, which generates ROS, is likely to influence the riskscore. GSVA analysis using a gene set that promotes ROS showed that the GSVA of ROS in TCGA colorectal cancer samples is lower than that in normal tissues (Fig. 9D). Correlation analysis revealed a negative correlation between the riskscore and GSVA of ROS (Fig. 9F), and colorectal cancer patients with higher GSVA of ROS scores have better DFI survival (Fig. 9E).

PDT generates ROS through photosensitizers. After incubating colorectal cancer cells HCT15 with the photosensitizer, ROS was detected, and the experimental results confirmed that PDT can promote ROS production in the colorectal cancer cell line HCT15 (Fig. 9G). Subsequently, QPCR was used to detect the expression of seven immune-related genes in HCT15 and HCT116 after PDT treatment. The results showed that PDT upregulates the expression of the ICOS gene in both cell lines, while downregulating the expression of HSPA1A, ULBP2, RBP7, OXT, SLC11A1, and INHBB (Fig. 9H and J). Correspondingly, the riskscore calculated based on QPCR-derived -∆T results indicated that PDT can downregulate the riskscore in colorectal cancer cell lines HCT15 and HCT116 (Fig. 9I and K). Combining these results suggests that PDT may lower the riskscore of colorectal cancer through ROS.

PDT affects tumor immune function by modulating riskscore

Previous studies have shown that immune cell populations promote and limit tumorigenesis through the production of ROS31,32. In addition, ROS are involved in various steps of tumor cell metastasis, including migration and invasion, detachment from normal cell contact, survival in the circulation and establishment of growth in distant tissues33,34. To further investigate the relationship between riskcore and metastasis, we used the set of genes associated with ROS generation (GOBP_POSITIVE_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS.v2024.1.Hs.gmt) to assess the distribution characteristics of ROS in the immune microenvironment of normal colorectal tissues, primary lesions and liver metastases. Through co-localisation analysis, we investigated the relationship between ROS and key immune gene sets associated with metastasis. We found that the content of ROS scores was significantly increased in liver metastatic lesions compared with normal colorectal tissues and primary lesions, and the co-localisation of ROS with metastasis-associated key immune genes riskscore in liver metastatic lesions showed that ROS scores were mutually exclusive of riskscore, with higher ROS scores corresponding to lower riskscores, a phenomenon that was particularly evident in macrophages (Supplementary Fig. 7). Based on the above results, ROS affects tumor immunity in colorectal cancer, particularly in macrophages. Therefore, PDT also has the potential to regulate tumor immunity, and this effect is closely related to the changes in riskscores of various cell types in colorectal cancer.

The above results indicate that the proportion of macrophages in the immune microenvironment gradually increases in the order of normal tissue, primary colorectal tumors, and colorectal cancer liver metastases. Additionally, the results of functional enrichment analysis based on single-cell sequencing revealed a close correlation between macrophages and the activation of T cells (Fig. 10A). Previous studies have shown that macrophages can promote tumor development by suppressing cellular immunity to tumors via the PD-1 pathway35,36,37,38. The results of the present study showed that the trend of PD-1 and PD-L1 in the three types of tissues was consistent with the trend of macrophage proportion, with the highest expression and proportion of expression in liver metastases (Fig. 10B). This result suggests that macrophage-mediated PD-1 immunosuppressive pathway plays an important facilitating role in colorectal cancer ductal liver metastasis. Analysis of transcriptomic data by TCGA revealed that riskscores were positively correlated with the expression of PD-1 and PD-L1 in colorectal cancer (Fig. 10C). Whereas PD-1 was found to be predominantly expressed in T cells by single-cell microenvironmental analysis, PD-L1 was predominantly expressed in fibroblasts and macrophages. In addition, the range and expression of PD-1 and PD-L1 were higher in the high-risk scoring group than in the low-risk group, regardless of the type of cells (Fig. 10D, E).

Fig. 10
figure 10

PDT affects the immune function of CRC. (A) Functional enrichment results of marker genes of various cells in the colorectal cancer immune microenvironment. (B) The trends of PD-1 and PD-L1 in the three types of tissues were consistent with the trends of macrophage proportion. (C) Riskscores were positively correlated with the expression of PDCD1 and CD274 in TCGA colorectal cancer mRNA sequencing data. (D,E) Relative to the low-risk group, the range and amount of PDCD1/CD274 expression were higher in the high-risk scoring group and were consistent with the expression of some immune cells. (F) IHC results of mouse MC38 cell line subcutaneous transplantation tumors. The expression and range of CD68 were significantly lower in the PDT group compared with the control group.

Full size image

In addition, previous immune infiltration results showed higher macrophages in the high-risk group than in the low-risk group, especially more pronounced in the M0 and M2 classes of immunosuppressive or non-functional macrophages. To further validate the effect of PDT on the whole tumor microenvironment, we constructed a subcutaneous tumor model of the mouse colorectal cancer cell line MC38 and treated it with PDT. The results showed that the expression and range of CD68 in the tumor tissues were significantly reduced in the photodynamic treatment group compared with the control group, suggesting that photodynamic therapy could reduce the proportion of macrophage infiltration in colorectal cancer (Fig. 10F). Together, these results suggest that PDT can improve the tumor immune function in colorectal cancer by reducing the expression and infiltration of PD-L1 in macrophages.

Analysis of immune-related genes regulated by PDT

Interaction and synergy between different genes are essential to understand the complexity of the disease and the molecular mechanisms of disease development39. We subjected seven genes in the IRGP model to interaction analysis, of which SLC11A1 is the core gene (Supplementary Fig. 8A). We show the relationship between gene expression and cells. In macrophages, SLC11A1 had the highest expression compared to other genes. In T cells, ICOS had the highest expression compared to other genes. HSPA1A was highly expressed in macrophages, fibroblasts, B cells, T cells, epithelial cells, tissue stem cells and endothelial cells (Supplementary Fig. 8B). Next, we performed gene expression analysis by different tissue subgroups. We found that SLC11A1 and ICOS were most highly expressed in CLM compared to normal colorectal tissues and primary tumors (Supplementary Fig. 8C). The highest expression of HSPA1A was found in primary tumors (Supplementary Fig. 8C). Therefore, we wanted to verify the expression levels of these three genes in different tissues. We verified the expression levels of SLC11A1, ICOS and HSPA1A in primary colorectal cancer and paracancerous tissues, and liver metastases by immunohistochemistry. We found that HSPA1A expression was highest in primary tumors, whereas SLC11A1 and ICOS expression was highest in liver metastases (Supplementary Fig. 8D), which was consistent with the database analysis.

To further elucidate the intrinsic effect of PDT on tumor growth, we performed experimental validation. EdU assay was performed to detect the effect of PDT on colorectal cancer cell proliferation. The results showed that PDT reduced the proliferative ability of HCT15 and HCT116 cells compared with the control group (Fig. 11A, B). In addition, the PDT-treated cells exhibited a significant wound healing rate (Fig. 11C-E). We also measured the mouse subcutaneous tumor model as previously described. The results showed that subcutaneous tumors were significantly reduced in mice treated with PDT compared to controls (Fig. 11F-H). Immunohistochemical assays showed significantly higher ICOS expression and significantly lower SLC11A1 and HSPA1A expression in the PDT group (Fig. 10I). Together, these results suggest that PDT may treat CLM by regulating immune-related prognostic genes.

Fig. 11
figure 11

PDT inhibited the proliferation and migration of colorectal cancer cells, increase the expression of ICOS and decrease HSPA1A and SLC11A1 expression. (A,B) The proliferation capacity was measured by EdU assay in HCT15 and HCT116 cells with PDT. (C–E) Wound healing was used to detect the ability of migration of HCT15 and HCT116 cells after PDT. (F) The mice were randomly divided into control group and PDT group, with 4 mice in each group. Intervention methods are described in Materials and methods. After 14 days, mice were sacrificed and tumors were shown. (G,H) The weight and tumor volume of the mice were recorded every 2 days. (I) The expression of HSPA1A, ICOS and SLC11A1 was analyzed by IHC staining (Bar, 50 μm). *p < 0.05, ***p < 0.001, ****P < 0.0001.

Full size image

Discussion

With the rapid development of technologies such as whole transcriptome sequencing, single-cell sequencing and artificial intelligence, in-depth study of the pathogenesis of cancer has become a reality40. More and more studies have focused on the tumor immune microenvironment, and immune-related genes have been found to serve as prognostic markers for survival outcomes or response to immunotherapy in tumor patients41. Zhang et al.42 successfully established a prognostic model consisting of nine mast cell-related genes by analysing scRNA-seq and bulk RNA-seq data. This model provides an important guide for selecting immunotherapy strategies and predicting treatment efficacy in patients with lung adenocarcinoma. Ren et al.43 found that a composite signature based on cancer-associated fibroblast features could more accurately stratify patients and predict their responsiveness to immunotherapy for lung adenocarcinoma. In addition, Zhang et al.44 constructed an immunogenic cell death-related signature that could not only effectively predict the prognosis of lung adenocarcinoma patients and their response to immunotherapy, but also guide clinical treatment to a certain extent. Another study reported a prognostic model based on nine cancer-associated fibroblast-associated genes, which independently predicted the overall survival of colorectal cancer patients, and provided new genomic evidence for immunotherapeutic strategies targeting cancer-associated fibroblasts45. In this study, we constructed a prognostic model of immune-related genes, which included seven differentially expressed metastasis-associated immune genes: HSPA1A, ULBP2, RBP7, OXT, SLC11A1, INHBB, and ICOS. This model accurately predicted DSS, DFI, and PFI in CLM patients, showing excellent predictive accuracy. In addition, based on the seven-gene immunoprognostic model, a nomogram was established to predict the prognosis of CLM patients, which is expected to provide a theoretical basis for the accurate assessment of the prognosis of CLM. These results not only enrich our understanding of the tumor immune microenvironment, but also lay a solid foundation for future personalized treatment strategies.

The tumor immune microenvironment contains different types of cells, and its complex intercellular communication is closely related to tumor metastasis46,47. We inferred interactions between different cell populations in high and low groupings based on CellChat. Most interactions occurred between fibroblasts, epithelial cells and macrophages. In the low-risk group, collagen signalling between fibroblasts and epithelial cells was reduced. Notably, the strength of interactions between macrophages and fibroblasts and between fibroblasts and epithelial cells was significantly reduced in the low-risk group. Collagen signalling also plays a key role in tumor metastasis. Studies have shown that type IV collagen-activated signalling is a key and widely relevant cue for tumor cell survival and growth in the liver48,49. We also found that macrophages in the low-risk group interacted less with epithelial IL1 signalling than those in the high-risk group. Consistent with our study, Koncina et al. found that IL-1R1+ cancer-associated fibroblasts (CAFs) influence colorectal cancer development and immune responses through the IL-1 receptor signalling pathway. The use of the IL-1 receptor antagonist Anakinra reduced colorectal tumor growth in vivo50. This suggests a possible inhibitory effect on tumor development by interfering with the IL-1 signalling pathway. This may be due to fibroblast-macrophage interactions that interfere with macrophage IL1 secretion, leading to tumor regression. In addition, MHC-I signalling and MHC-II signalling were reduced in the high-risk group compared to the low-risk group. Indeed, tumors can evade immune surveillance by reducing the expression of MHC molecules30. Indeed, tumors can evade immune surveillance by reducing the expression of MHC molecules. Specifically in CLM, it was found that MHC-I expression was strongly associated with tumor outcome after complete resection51. Thus, interactions between fibroblasts, epithelial cells and macrophages play a crucial role in liver metastasis and may contribute to the discovery of promising therapeutic targets.

Single-cell sequencing analysis showed significant differences in the tumor immune microenvironment between the groups. Macrophages were more abundant in the high-risk group, whereas B cells, T cells and epithelial cells were more abundant in the low-risk group. Notably, macrophages were also more abundant in liver metastases compared to primary and normal tissues. Normal colorectalic tissues had higher levels of B cells. B cells are an important component of the humoral immune response and are responsible for the production and release of antibodies against specific antigens and play an important role in the anti-tumor response52. Macrophages are key in the production of inflammatory cytokines and activation of the immune system53. These data demonstrate the accuracy of our model, and therefore, our IRGP model will enable researchers to develop more effective therapeutic regimens to enhance the immune response. It also provides new information for individualised immunotherapy in CLM and improves clinical prognosis.

Gene causality inference revealed that SLC11A1, ICOS and HSPA1A are key immune prognostic genes, with SLC11A1 being the most central. The highest expression of SLC11A1 was found in macrophages. ICOS had the highest expression in T cells compared to other genes. HSPA1A was highly expressed in macrophages, fibroblasts, B cells, T cells, epithelial cells, tissue stem cells and endothelial cells. The expression of SLC11A1 and ICOS was higher in liver metastases than in normal colorectal tissue and primary tumors. HSPA1A expression was highest in primary tumors. HSPA1A expression was highest in primary tumors, whereas SLC11A1 and ICOS expression was highest in liver metastases, which was consistent with the results of database analysis. Based on these analyses, we wanted to verify whether our data could be combined and applied with PDT. PDT is a potential treatment for CLM based on the local activation of photosensitisers within the tumor, inducing chemical damage that leads to tumor cell death. The immune prognostic score we constructed is closely related to ROS, and PDT is a novel tumor therapy that kills cancer cells by generating ROS, which may modulate key immune genomes associated with metastasis34. In addition, our results confirmed that SLC11A1, ICOS and HSPA1A were significantly altered in CRC cells after PDT. SLC11A1, also known as natural resistance-associated macrophage protein-1, was initially reported to be used to fight against a variety of pathogens. As a member of the solute carrier family, SLC11A1 plays a role in innate immunity, autoimmune diseases and infections54,55. Susceptibility to infections and autoimmune diseases has been associated with SLC11A1. Shahzad et al. reported that in patients with tuberculosis, low expression and variants of SLC11A1 may impair the immune response to tuberculosis56. A recent study showed that SLC11A1 has multiple effects on macrophage activation and plays an important role in the immune response57. In the present study, we also found high expression of SLC11A1 in macrophages, and the content of macrophages in liver metastases was significantly higher than that in normal colorectal tissues and primary tumors, suggesting that SLC11A1 in macrophages plays an important role in the occurrence and development of liver metastases, and may be a new therapeutic target for CLM. ICOS, a co-stimulatory receptor that promotes the proliferation of T-cells57, is closely related to inhibitory activity of regulatory T cells (Tregs) and can influence the functional activity of a wide range of immune cells that regulate the adaptive immune system58. Notably, the ICOS/ICOSL axis plays a dual role, which may be related to its involvement in anti-tumor T-cell responses and promotion of tumor responses. Thus, both antagonist and agonist antibodies may be therapeutic targets for tumor therapy against the ICOS/ICOSL pathway. Indeed, ICOS activation enhances the effect of inhibitory checkpoint blockade, and monoclonal antibodies against ICOS not only reduce the function of immunosuppressive Tregs but also inhibit lymphoid tumor cells expressing Tfh markers59. Two agonist and one antagonist monoclonal antibodies are currently being evaluated in phase I/II trials. However, the efficacy and safety of these therapies, as well as the combination strategy with anti-ICOS agonists or antagonists, remain unclear. HSPA1A is a protective protein in the heat shock family of proteins that is upregulated after stress60. As a molecular chaperone that regulates stress and tumor cell survival, HSPA1A also localises and embeds itself in the plasma membrane of stressed and tumor cells61. Interaction of HSPA1A with specific lipid membranes allows HspA1A to modulate many important functions such as immunity, membrane stability, autophagy and apoptosis62. Therefore, understanding and manipulating the surface expression of HSPA1A is a promising therapeutic approach.

Our study still has some limitations. Firstly, our study is only a bioinformatics analysis and the prognostic model needs to be further validated in clinical practice. Second, further experiments are needed to verify the specific roles of these 7 immune-related genes in CLM. In conclusion, through comprehensive analysis, we established a prognostic model of CLM consisting of 7 immune-related genes. This model not only assesses the prognosis of CLM patients, but also determines the level of immune cell infiltration and intercellular communication mechanisms. Importantly, some immune-related genes showed significant expression differences after photodynamic therapy treatment, suggesting that these genes may play a promotional or antagonistic role in the treatment of CLM with photodynamic therapy, but further experimental validation is needed. This not only provides new ideas for potential immunotherapy strategies, but also lays a solid foundation for the study of CLM and tumor microenvironment.