The gene expression data and related clinical information for The Cancer Genome Atlas (TCGA)-colon adenocarcinoma (COAD) were sourced from TCGA database (https://portal.gdc.cancer.gov/). Immune-related genes (IRGs) were downloaded from the IMMPORT (https://www.immport.org/home) database, while PCD-related genes (PRGs) were referenced from the recent study by Zou et al (7). Genetic mutation data were obtained from the Genomic Data Commons (GDC; https://portal.gdc.cancer.gov/) database and analyzed using the ‘maftools’ package. Transcriptome datasets GSE17538, GSE29621, GSE39582, GSE44076 and GSE74602 were retrieved from the Gene Expression Omnibus (GEO) database. Additionally, GSE161277 was downloaded for single-cell RNA sequencing (scRNA-seq) data. Due to the small sample sizes of GSE39582, GSE44076 and GSE74602, they were merged and normalized into the GEO-Meta dataset for model validation.

Zhou et al (8) developed the apoptosis index and epithelial-mesenchymal transition (EMT) index to assess the prognosis of COAD. Based on their study, the Immune Potential Index (IPI) and PCD potential index (PPI) were introduced in the present study. Using the TCGA-COAD dataset, immune and PCD genes were utilized as independent variables, with overall survival (OS) as the target variable, to identify significant genes associated with prognosis through univariate Cox regression analysis. In the TCGA-COAD dataset, differentially expressed genes (DEGs) between tumors and normal tissues were first identified. These DEGs were then combined with significant genes from the univariate Cox regression analysis to pinpoint IRGs and PRGs. At the same time, these genes were intersected with those from the GSE17538, GSE29621, GSE39582, GSE44076 and GSE74602 datasets to identify IRGs and PRGs that present across all datasets.

Principal component analysis (PCA) was employed to quantify the levels of immunity and PCD in each tumor tissue, with main components 1 and 2 serving as the primary dimensions. Referring to previous research (9), the IPI and PPI were defined respectively, where the formula for calculating IPI or PPI is ∑ (|PC1i| + |PC2i|), with i representing the gene expression matrix of IRGs or PRGs. Patients in the TCGA-COAD and GSE17538 datasets were divided into two molecular subtypes based on the median values of IPI: IPI_high and IPI_low, while they are also divided into another two molecular subtypes based on the median values of PPI: PPI_high and PPI_low. Combing IPI and PPI, patients were divided into four molecular subtypes: IPI_high + PPI_high, IPI_high + PPI_low, IPI_low + PPI_high, and IPI_low + PPI_low. Kaplan-Meier (K-M) curve analysis was used to compare OS differences among these four molecular subtypes.

Based on the merged expression matrix of IRGs and PRGs, clustering analysis on the GSE17538 and TCGA-COAD datasets was conducted using the R package ‘ConsensusClusterPlus’. The association between the identified clustering categories and OS was then investigated.

To develop a high-precision and stable risk assessment model, the study of Liu et al (10) was referenced, and a combination of 101 algorithms was constructed by integrating 10 machine learning algorithms. The machine learning algorithms include random survival forest, elastic network, least absolute shrinkage and selection operator, Ridge, stepwise Cox, CoxBoost, partial least squares regression for Cox, supervised principal components, generalized boosted regression modeling and survival support vector machine. Details of the hyperparameters used in the machine learning models are provided in Table SI. A leave-one-out cross-validation framework was employed to fit these 101 prediction models, using a merged gene set of IRGs and PRGs as independent variables, with the TCGA-COAD dataset serving as the training set. At the same time, the consistency index (C-index) for each model was calculated across all validation datasets (GSE17538 and GSE29621). The model with the highest weighted average C-index in the validation set was selected as the best model. The code for model training and validation has been uploaded to https://github.com/wuqiumei/101-ML.git.

Using this optimal model, a risk score was calculated for each sample, referred to as ‘Score’. Samples from different datasets were divided into high and low groups based on the median value of Score. The value of Score was validated in the GSE17538, GSE29621 and TCGA-COAD datasets. Based on the TCGA-COAD dataset, a nomogram integrating the Score with other clinical factors was constructed to evaluate its predictive ability for 1-year, 3-year and 5-year survival rates. Additionally, the Score was explored in relation to various aspects such as TME, immune markers and genetic mutations based on the TCGA dataset.

The ‘estimate’ R package is used to calculate the Immune Score and Stromal Score of patients. The ‘xCell’ R package is used to assess the infiltration of various immune cell types in both groups with high and low scores. Additionally, Bagaev et al (11) proposed 29 functional gene expression signatures (Fges) to represent the cellular and functional attributes of the TME.

Based on the TCGA-COAD dataset, the Tumor Immune Dysfunction and Exclusion (TIDE) tool (http://tide.dfci.harvard.edu/login/) was employed to assess the potential individual response to ICI treatment (12). Individual scores such as Exclusion score, myeloid-derived suppressor cells' (MDSCs) score, tumor-associated M2 macrophages' (TAM) score and cancer-associated fibroblasts' (CAFs) score were also obtained.

Based on the TCGA dataset, the differential expression of genes included in the model between tumors and normal tissues was investigated. Underlying reasons for these differential expression levels from the perspectives of genetic mutations, copy number variations (CNV) and methylation were explored. Additionally, the expression association of the included genes with tumor immune genes was explored from multiple aspects such as chemokines, immune inhibitors, Major Histocompatibility Complex (MHC), immune stimulants and receptor-related genes.

The single-cell sequencing dataset GSE161277 was retrieved from the GEO database. Low-quality cells were excluded based on the criteria of nFeature_RNA >500, nFeature_RNA <10,000, and percent.mt <25. The ‘harmony’ package was used to preprocess and integrate four single-cell samples. The CellCycleScoring function was used to calculate the cell cycle scores, and the influence of the cell cycle on clustering was subsequently removed. The FindVariableFeatures function in the ‘Seurat’ R package was used to identify 2,000 highly variable genes for subsequent dimension reduction based on PCA. According to principal components 1 to 20, the t-distributed stochastic neighbor embedding method was employed to visualize visualizing single-cell clustering. Subsequently, the ‘singleR’ R package was used to annotate cells based on different markers on the cell surface.

Patients who received treatment for CRC in the Department of Gastrointestinal Surgery at Zhanjiang Central People's Hospital from November 2021 to November 2022 were selected as research subjects. The patients ranged in age from 34 to 78 years, with male and female patient percentages of 72.5 and 27.5%, respectively. Inclusion criteria were as follows: i) Age ≥18 years; ii) Confirmed diagnosis of CRC through pathological examination, without other tissue-originated malignant tumors; iii) no peripheral organ metastasis before surgery and underwent radical tumor surgery with 1-year follow-up. Patients with clear peripheral organ metastasis or recurrence within 1 year of follow-up, as evidenced by imaging or pathology, were classified into the transfer group. Otherwise, they were classified into the non-transfer group. From the transfer group, 10 patients were randomly selected and matched with 10 patients from the non-metastasis group based on propensity scores considering factors such as sex, age, body mass index and tumor staging. The operation process of propensity score matching was completed using SPSS 21.0 software (IBM Corp.). CRC specimens from these two groups were obtained retrospectively from the specimen library for immunohistochemistry (IHC) and western blotting to assess LTB4R expression. The present study adhered to the Helsinki Declaration and was approved by the Institutional Review Committee of Zhanjiang Central People's Hospital (approval no. 2023-002-013; Zhanjiang, China). Written informed consent was provided by all participants.

CRC tissue samples were collected, including 10 cases from the transfer group and 10 matched cases from the non-transfer group. The tissue was typically embedded in paraffin and sectioned (4-µm thick). Deparaffinization was performed using xylene, followed by rehydration with a gradient of alcohols. The sections were then placed in TRIS-EDTA antigen retrieval solution, covered, and microwaved for 10 min on high power and 5 min on medium power. After stopping the heating, the sections were cooled to room temperature and washed several times with distilled water. Endogenous peroxidase activity was blocked using 3% H₂O₂ at room temperature for 10 min. The sections were incubated overnight at 4˚C with a primary antibody against LTB4R (1:25; cat. no. R30118; Zen-Bio, Inc.). The next day, the sections were rewarmed at room temperature for 45 min, washed three times with PBS, and incubated with a biotinylated secondary antibody. The staining was visualized using DAB, followed by counterstaining with hematoxylin and differentiation with hydrochloric acid alcohol. Finally, the sections were dehydrated in a gradient of ethanol, cleared in xylene, and mounted for microscopic examination.

Total protein was extracted from CRC tissue. Protein concentration was determined using the BCA method, and the samples were stored at -80˚C. A total of, 20 µg proteins were separated by 10% SDS-PAGE gel electrophoresis and transferred onto a PVDF membrane using a constant current of 300 mA. The membrane was blocked with a blocking solution at room temperature for 1 h. The blocking solution was prepared by dissolving non-fat milk in TBST (containing 0.05% Tween-20) to a final concentration of 5%. The membrane was then incubated overnight at 4˚C with diluted primary antibodies: LTB4R (1:1,000) and GAPDH (1:50,000; cat. no. 60004-1-Ig; Proteintech Group, Inc.). After washing with TBST, the membrane was incubated with an HRP-conjugated secondary antibody (1:50,000; cat. no. BL003A; Biosharp Life Sciences) at room temperature for 1-2 h. Following additional washes, chemiluminescent detection was performed using a chemiluminescence substrate (cat. no. BL520A; Biosharp Life Sciences). Protein band intensities were analyzed using ImageJ software (version: 1.8.0; National institutes of Health), with GAPDH as the internal reference for comparing protein expression differences between different groups.

The Wilcoxon test was used to compare differences between continuous variables. while the chi-square test was employed for inter-group comparisons of categorical variables. The log rank test was utilized to analyze inter-group differences in survival rates. Depending on the data characteristics, Pearson correlation coefficient or Spearman correlation coefficient was used for correlation analysis. P<0.05 was considered to indicate a statistically significant difference (^*P<0.05; ^**P<0.01; ^***P<0.001; ^****P<0.0001).

A total of 2,483 IRGs were downloaded from the IMMPORT database (Table SII), and referred to the research by Zou et al (7) to obtain 1277 genes related to 12 types of PCD (Table SIII). Using the TCGA-COAD dataset, univariate Cox analysis was performed with OS time as the target variable and 135 immune genes (Table SIV) and 102 PCD genes (Table SV) that are related to prognosis were identified. The final IRGs and PRGs were obtained by selecting genes that exist in all datasets including GSE17538, GSE29621, GSE39582, GSE44076 and GSE74602.

Based on the TCGA-COAD dataset, the heatmaps of DEGs were generated for IRGs and PRGs between tumor tissues and adjacent normal tissues (Fig. 1A and B). Similarly, heatmaps of DEGs between tumors and adjacent tissues were obtained from the GEO-Meta dataset (Fig. 1C and D). After combining the gene expression matrix for IRGs and PRGs, consensus clustering analysis on the GSE17538 dataset was performed (Fig. 1G). The results revealed significant differences in OS rates among different clusters (Fig. 1H). The same analysis conducted on the TCGA-COAD dataset (Fig. 1E and F) also demonstrated significant differences in OS among different clusters (P<0.05). These findings suggested that the gene set comprising IRGs and PRGs can effectively distinguish the overall prognosis of colon cancer to a certain extent.

To further investigate the crosstalk between IRGs and PRGs, PCA was performed to quantify the immune level and PCD level in each tumor tissue, which are defined as IPI and PPI, respectively. The distributions of IPI and PPI are shown in Fig. 2A and B; the correlation analysis between IPI and PPI is presented in Fig. 2C. Survival analysis of the TCGA-COAD dataset revealed that the IPI_high group had a significantly lower survival rate compared with the IPI_low group (P<0.01; Fig. 2D). Similarly, the PPI_high group had a significantly reduced survival rate compared with the PPI_low group (P<0.01; Fig. 2E), indicating that both high PPI and high IPI are associated with poorer survival outcomes.

Based on IPI and PPI, patients in the TCGA-COAD dataset were categorized into four molecular subtypes: IPI_high + PPI_high, IPI_high + PPI_low, IPI_low + PPI_high, and IPI_low + PPI_low. K-M analysis demonstrated that the survival rates of the IPI_high + PPI_high, IPI_high + PPI_low, and IPI_low + PPI_high subtypes were significantly lower than those of the IPI_low + PPI_low subtype (P<0.001; Fig. 2F). After merging the first three subtypes into a single category called the ‘Others group’, a comparison with the IPI_low + PPI_low subtype revealed consistent results (Fig. 2G). The results of the multivariate Cox regression analysis (Fig. S1) based on TCGA-COAD data indicated that this molecular classification is an independent prognostic factor, unaffected by variables such as age, tumor stage, and microsatellite status. Survival analysis results for IPI and PPI based on the GSE17538 dataset are depicted in Fig. 2H and I, respectively. As shown in Fig. 2J, the IPI_low + PPI_low group has a higher survival rate. However, due to limitations in the dataset samples, this difference has not yet reached statistical significance. The analysis across various datasets indicates a crosstalk between immune and PCD. The combination of IPI and PPI effectively distinguishes OS rates and other prognostic factors in patients with CRC.

After merging IRGs and PRGs, an expression matrix containing 104 genes was constructed. Using this combined gene set as the independent variable and the TCGA-COAD dataset as the training set, 101 prediction models were fitted, and the C-index was calculated for each model in the validation datasets (GSE17538 and GSE29621). The performance of these prediction models is presented in Fig. 3A, which reveals that the CoxBoost + Ridge model achieved the highest weighted average C-index in the validation set. This model includes 13 gene features: VEGFA, CCRL2, LTB4R, RORC, INHBB, FABP4, NOS2, C8G, CD24, HSPB1, SFRP2, STK25 and NOL3. The risk score for each sample was calculated. Using the median value of Score, samples from different datasets were classified into high and low groups. The K-M analysis results for the TCGA-COAD, GSE17538 and GSE29621 datasets are depicted in Fig. 3B-D respectively, demonstrating that the prognosis of the high group was significantly worse. To explore the performance of Score in predicting OS in the GSE17538 dataset, it was found that AUC for 1-year, 3-year and 5-year OS was 0.792, 0.716 and 0.758, respectively (Fig. 3F). In the GSE29621 dataset, the AUCs for 1-year, 3-year and 5-year survival rates were 0.645, 0.604 and 0.623, respectively (Fig. 3E). Simultaneously, a multivariate Cox regression analysis based on the TCGA-COAD dataset revealed that the Score is a significant prognostic factor, independent of tumor grading, age and MSI (Fig. 3G). A nomogram incorporating the Score and other clinical factors was also constructed (Fig. 3H). The receiver operating characteristic curve areas for predicting 1-year, 3-year, and 5-year OS were 0.759, 0.797 and 0.791, respectively (Fig. 3I), indicating strong predictive performance.

TME is a complex environment comprising tumor cells, immune cells and stromal cells, all of which significantly affect tumor cell invasion and metabolism (13). Analysis of the TCGA-COAD dataset using the ESTIMATE algorithm revealed the Immune Score of the low Score group was significantly higher than that of the high Score group (P<0.05, Fig. 4A), suggesting greater immune cell infiltration in the low Score group. As shown in Fig. 4B, the Stromal Score of the low Score group was significantly lower than the high group (P<0.01). Using the XCELL algorithm to assess and compare immune cell infiltration between the two groups, it was observed that the low Score group had higher abundances of activated dendritic cells, B cells, CD4⁺ memory T cells, CD8⁺ naive T cells, CD8⁺ T cells, CD8⁺ Tcm cells, CD4⁺ T cells, CD4⁺ Tem cells, CD8⁺ memory T cells, natural killer (NK) cells, and others (Fig. 4C). Further correlation analysis indicated that the Score was negatively associated with the infiltration of 12 types of immune cells, including macrophages, activated dendritic cells, memory B cells, T cells, Th2 cells, CD8⁺ Tcm cells and NK cells (Fig. 4D). This analysis supports the conclusion that a higher Score is linked to reduced immune cell infiltration.

The scores of 29 Fges were compared between the two groups, and it was found that the low Score group exhibited significantly higher scores in antitumor immune-related Fges, including M1 signature, myeloid immune suppression, checkpoint molecules, Th2 signature, neutrophil signature, antitumor cytokines, NK cells, Th1 signature and T cells (Fig. 4E). The TCGA-COAD dataset defines three immune phenotypes: Desert phenotype (phenotype 1 and 2), excluded phenotype (phenotype 3 and 4) and inflamed phenotype (phenotype 5 and 6). As depicted in Fig. 4F, the proportions of desert phenotype (40%) and excluded phenotype (33%) in the high Score group were higher than those in the low Score group (48 and 22%, respectively), while the proportion of inflamed phenotype (31%) was lower compared with the low Score group (26%). A chi-square test performed on the two groups of immune phenotypes yielded a P-value of 0.036, indicating a significant difference in the proportion of immune phenotypes between the two groups (P<0.05).

The present study further investigated the association between Score and immunotherapy response. Given the correlation between immunotherapy response and MSI, the association between Score and MSI was explored. As demonstrated in Fig. 4G, there was no significant difference in Score between the microsatellite stable (MSS)/MSI low subtype and the MSI high (MSI-H) subtype in the TCGA-COAD dataset. The high Score group exhibited significantly higher scores for MDSC, CAF, TAM.M2, Exclusion and Dysfunction compared with the low Score group (P<0.05; Fig. 4H-L). This indicates that the high Score group is more prone to tumor immune evasion. Furthermore, the TIDE score of the high Score group was significantly higher than that of the low Score group (P<0.001; Fig. 4M), suggesting a lower response to ICI treatment in the high Score group.

Somatic mutations in the high and low Score groups in the TCGA-COAD dataset were further investigated. The waterfall plots for the two groups of mutated genes are shown in Fig. 5A and B, respectively. The results revealed that somatic mutations were present in 93.97% of samples in the low Score group and 98.02% of samples in the high Score group. The top 8 genes with the highest mutation frequencies were identical in both groups, including APC, TP53, PIK3CA, TTN, CRAS, MUC16, SYNE1 and FAT4. As depicted in Fig. 5C, there were notable differences in somatic gene mutations between the two groups. In the high Score group, mutations in TP53 (P<0.001), SNX1 (P<0.01), ZNF214 (P<0.01), ELAPOR1 (P<0.01), AP2M1 (P<0.01), LMX1A (P<0.01) and NACC2 (P<0.01) were notably higher compared with the low Score group. On the other hand, mutations in PCDH10 (P<0.001), WDR19 (P<0.001), SCN5A (P<0.01), ATP2 (P<0.01) and EFCAB13 (P<0.01) were greatly increased in the low Score group. The tumor-related signaling pathways were enriched by the mutated genes in the high and low Score groups, respectively 9 Fig. 5D and E). Both groups exhibited enrichment in largely the same signaling pathways, including RTK-RAS, NOTCH, WNT, Hippo, PI3K, Cell Cycle, MYC, TGF-β, TP53 and NRF2. Regarding tumor mutation burden, although the high Score group had a higher burden than the low Score group, this difference did not reach statistical significance (Fig. 5F).

In total, 13 feature genes included in the model were analyzed. Using the expression matrix of these genes, the consensus clustering method was applied to the TCGA-COAD dataset. The optimal number of clusters k was determined to be 2, resulting in the division of the dataset samples into 2 clusters (Fig. 6A). The comparison of the OS time between different clusters is presented in Fig. 6B. Similarly, the GSE17538 dataset was clustered using the same method, and the optimal number of clusters k was also found to be 2 (Fig. 6C). Although differences in OS time between the clusters were observed, they did not reach statistical significance (Fig. 6D). In the TCGA-COAD dataset, K-M analysis of the 13 feature genes indicated that all were related to patient survival rate (Fig. 6E). Univariate Cox analysis revealed that VEGFA, LTB4R, INHBB, FABP4, C8G, HSPB1, STK25, and NOL3 were significant risk factors, whereas the remaining 5 feature genes were identified as significant protective factors (Fig. 6F). In the GSE17538 dataset, univariate Cox analysis showed that VEGFA, FABP4, HSPB1 and SFRP2 were significant risk factors, while RORC and NOS2 were significant protective factors (Fig. 6G). The differential expression of feature genes between tumor tissues and normal tissues in the TCGA-COAD dataset were further investigated. In cancer tissues, genes such as VEGFA, LTB4R, INHBB, STK25 and NOL3 were significantly upregulated, while protective factors including CCRL2, RORC, CD24 and SFRP2 were significantly downregulated (Fig. 6H). In the GSE17538 dataset, risk factors VEGFA and HSPB1 were significantly upregulated in cancer tissues, while protective factors CCRL2, RORC, NOS2, C8G, CD24 and STK25 were significantly downregulated (Fig. 6I). To explore the driving factors of abnormal expression of these feature genes in CRC tissues, the correlation was analyzed between gene expressions and mutations, methylation, and CNVs using the cBioPortal database. The results of this correlation analysis are shown in Fig. S2A-C. Additionally, the correlation analysis of feature gene expression is presented in Fig. S2D. The results indicated a broad positive correlation among VEGFA, LTB4R, INHBB, STK25 and NOL3, while CCRL2, RORC and CD24 also showed a positive correlation.

In the present study, the correlation between the expression of the 13 selected feature genes and various immune-related gene categories was analyzed, including chemokines, immune inhibitors, MHC-related genes, immune stimulation-related genes and receptor-related genes. The correlation analysis results, depicted in Fig. 7A-E, reveal a significant association between the expression of these 13 feature genes and these five types of immune genes. Specifically, genes such as SFRP2, NSO2 and FABP4 exhibit a positive correlation with the expression of most chemokines, immune inhibitors, MHC-related genes, immune stimulation-related genes and receptor-related genes. Conversely, the expression of genes like LTB4R shows a negative correlation with the expression of the majority of these IRGs.

Preprocessing, integration, dimensionality reduction and clustering analysis was performed on GSE161277. Using the R package ‘singleR’, cells were annotated based on various cell surface markers, ultimately identifying 5 distinct cell types: epithelial cells, T cells, B cells, monocytes and endothelial cells (Fig. 8A). Comparison of the Score for each cell type (Fig. 8B) reveals obvious differences. Endothelial cells scored the highest, while B cells scored the lowest. This suggests that endothelial cells and B cells may play crucial roles in the crosstalk between cellular immunity and PCD. The expression levels of 13 feature genes in each cell subgroup are illustrated in Fig. 8C. Notably, VEGFA, STK25 and NOL3 are predominantly expressed in epithelial cells, while CD24 and other genes are mainly expressed in B cells and epithelial cells.

IHC analysis revealed that the staining of LTB4R was markedly more intense in colon cancer tissues in the transfer group (Fig. 9A), while the staining was lighter in the non-transfer group (Fig. 9B). The results of the IHC analysis showed that the expression of LTB4R in the transfer group was significantly higher than that in the non-transfer group. Western blotting analysis showed that the LTB4R protein expression level was significantly higher in the transfer group compared with the non-transfer group (Fig. 9C and D), with a statistically significant difference (P<0.01) (Fig. 9E). These findings suggested that LTB4R expression levels may be useful in distinguishing between CRC recurrence and metastasis.