Abstract Background Tumor microenvironment (TME), particularly immune cell infiltration, programmed cell death (PCD) and stress, has increasingly become a focal point in colorectal cancer (CRC) treatment. Uncovering the intricate crosstalk between these factors can enhance our understanding of CRC, guide therapeutic strategies, and improve patient prognosis. Methods We constructed an immune-related cell death and stress (ICDS) prognostic model utilizing machine learning methodologies. Furthermore, we performed enrichment analyses and deconvolution algorithms to elucidate the complex interactions between immune cell infiltration and the processes of PCD and stress within a substantial array of transcriptomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus data base (GEO) related to CRC. Single-cell sequencing and biochemical experiments were used to validate the interaction between the model genes and programmed cell death in tumor cells. Results The ICDS prognostic model exhibited robust predictive performance in seven independent cohorts, revealing an inverse correlation between model scores and patient prognosis. Meanwhile, the ICDS index was positively correlated with clinical stage. Model analysis indicated that patient subgroups with low ICDS index exhibited heightened immune activation features and elevated activity in PCD and stress pathways. Single-cell analysis further revealed that macrophages were the central drivers of immune characteristics underlying prognostic differences within the ICDS prognostic model. Pseudotime analysis and cellular experiments indicated that the model gene GAL3ST4 promotes the transition of macrophages toward an M2 pro-tumor phenotype. Furthermore, cell communication analysis and experimental validation revealed that the cuproptosis in tumor cells suppress GAL3ST4 expression, thereby inhibiting M2-like macrophage polarization. Conclusion In summary, we constructed the ICDS prognostic model and uncovered the mechanism by which tumor cells downregulate GAL3ST4 expression via cuproptosis to inhibit M2-like macrophage polarization, providing new targets and biomarkers for CRC treatment and prognosis evaluation. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-025-06143-9. Keywords: Colorectal cancer, Machine learning, Single cell RNA sequencing, Programmed cell death, Macrophage polarization Introduction Colorectal cancer (CRC) is the third most common cancer worldwide and the second leading cause of cancer-related deaths globally [[38]1]. The treatment for CRC currently includes a variety of methods, such as chemotherapy, targeted therapy and immunotherapy [[39]2, [40]3]. However, many patients fail to benefit from these therapies. Consequently, scientists worldwide are committed to exploring innovative diagnostic and therapeutic strategies and developing novel drugs to expand the range of treatment options [[41]4–[42]6]. Notably, cancer treatments can induce various types of programmed cell death (PCD), which significantly influence therapeutic outcomes through distinct mechanisms [[43]7, [44]8]. However, these mechanisms remain poorly studied and insufficiently understood. Tumor microenvironment (TME) is composed of various immune cells, including macrophages, fibroblasts, T cells, and neutrophils, along with cytokines and biochemical components such as the extracellular matrix, providing critical support for tumor initiation and progression [[45]9–[46]11]. The remodeling of TME is often accompanied by significant alterations in the states of immune cells, profoundly impacting cancer treatment outcomes, particularly immunotherapy [[47]9, [48]12].The remodeling of TME is regulated by multiple factors, with recent studies revealing extensive crosstalk between the tumor immune microenvironment, PCD and cellular stress [[49]13–[50]15]. Therefore, elucidating the mechanisms by which PCD and cellular stress contribute to immune microenvironment remodeling will provide fresh perspectives into how PCD influences cancer therapy. Certain forms of PCD, such as pyroptosis and necroptosis, are referred to as immunogenic cell death due to their ability to significantly impact immune cell states and reshape TME [[51]16]. In CRC, TME remodeling induced by immunogenic cell death not only enhances the efficacy of chemotherapy and radiotherapy but also significantly increases the response rate of microsatellite-stable (MSS) colorectal cancer patients to immunotherapy [[52]17–[53]19]. Despite recent advances in understanding the role of immunogenic cell death in regulating TME remodeling in CRC, such as immunogenic cell death triggering immune responses and reshaping the immune microenvironment through the release of damage-associated molecular patterns (DAMPs) interacting with pattern recognition receptors (PRRs) on immune cells [[54]20–[55]23], the specific mechanisms of different forms of immunogenic cell death in this process and their precise interactions with immune cell subsets in CRC remain largely unknown. In this study, we developed a prognostic model that integrates immunogenic cell death features and cellular stress characteristics using immune cell infiltration data. The model demonstrated excellent predictive performance across multiple independent cohorts, with model scores significantly negatively correlated with patient prognosis while revealing differences in immune responses and immunogenic cell death-related biological processes among patients with varying prognoses. Single-cell transcriptomics analysis revealed a strong association between model grouping and differences in GAL3ST4 expression levels as well as macrophage polarization states. Further experimental validation demonstrated that GAL3ST4 promoted the polarization of macrophages toward an M2 pro-tumor phenotype, a process regulated by tumor cell-mediated cuproptosis. By integrating multi-omics data and experimental analyses, we discovered that tumor cells suppress GAL3ST4 expression through the process of cuproptosis, thereby inhibiting the polarization of M2-like macrophages. These findings provided critical theoretical foundations and potential applications for understanding the mechanisms of CRC immune microenvironment remodeling and related therapeutic strategies. Methods Processing of Bulk RNA-seq data CRC RNA-seq data and clinical information from the TCGA database were retrieved via the UCSC Xena browser ([56]http://xena.ucsc.edu/). Gene expression data obtained from the GEO database ([57]GSE29621, [58]GSE39582, [59]GSE17536, [60]GSE17537, [61]GSE33382, and [62]GSE72970) were processed with log2(x + 1) transformation to acquire samples with normalized RNA-seq or microarray data. Single-cell RNA-seq data collection and analysis A single-cell RNA sequencing (scRNA-seq) dataset comprising 23 patients was obtained from the GEO database ([63]GSE132465). Cells with > 1,000 unique molecular identifier (UMI) counts > 200 genes < 6,000 genes, and < 10% of mitochondrial gene expression in UMI counts were filtered. Following the standard dimensionality reduction and clustering procedures, clusters were annotated with classical cell markers. For the sub-clustering analysis of myeloid cells and macrophages, we extracted cell populations initially annotated as myeloid cells and macrophages and repeated the aforementioned workflow. Myeloid cells were annotated based on classical myeloid markers, while macrophages were further annotated and classified according to their function. To differentiate between malignant and non-malignant epithelial cells, the ‘InferCNV’ R package was used to deduce large-scale chromosomal copy-number variations (CNVs) from single-cell transcriptome profiles. The R package ‘CellChat’ with default parameters was employed to compute the communication intensity between cuproptosis malignant epithelial cells and other cell types, while the ‘Monocle2’ package, also with default settings, was used to delineate the phenotypic transition trajectory of macrophages. Construction of the ICDS model gene set Selection of Immune-Related Genes: single-sample gene set enrichment analysis (ssGSEA) via the R package GSVA was employed to quantify the relative infiltration of 28 immune cell types in the TCGA-CRC cohort (97 genes for 28 immune cells downloaded from the Molecular Signatures Database (MSigDB)). Subsequently, clustering was performed using the ‘Consensus Cluster Plus’R package based on the infiltration profiles of these immune cells. Thereafter, a comprehensive assessment utilizing the consensus score matrix and CDF curve was conducted to determine the optimal number of clusters. The ‘WGCNA’R package was then employed to perform unsigned weighted gene co-expression network analysis (WGCNA) to identify immune-related gene modules (n = 4295). Selection of programmed cell death and stress genes: The R package ‘GSVA’ was used for ssGSEA to evaluate the enrichment of four pathways associated with programmed cell death and stress responses in the TCGA-CRC cohort. Clustering was then carried out with the ‘Consensus Cluster Plus’ R package based on the enrichment data of these pathways. The method for determining the optimal number of clusters is the same as described above. Finally, differential analysis was performed using the R package ‘limma’ (|logFC|> 0.2), resulting in the identification of 2,769 differentially expressed genes associated with programmed cell death and stress. Construction of the ICDS Model Gene Set: The final ICDS model gene set (n = 579) was established as the intersection of programmed cell death and stress genes with immune-related genes, deemed a reliable representation. Construction of the ICDS machine learning model We employed a three-step strategy to construct the ICDS machine learning model, with each patient's final score considered the ICDS index. Step 1: Univariate Cox Regression: Prognostic genes were identified through univariate Cox regression analysis, defining genes with hazard ratios (HR) and 95% confidence intervals (CI) with p < 0.05. Step 2: LASSO Cox Regression: To prevent overfitting, LASSO Cox regression was applied using the ‘glmnet’ R package to narrow down the candidate genes. Genes with the minimal lambda value were retained in the model. Step 3: Stepwise Cox Regression: Stepwise Cox regression was ultimately used to select the genes most suitable for model construction. The model with the minimal Akaike Information Criterion (AIC) for survival genes was designated as the final: [MATH: ICDS=∑< /mo>i=1n Coefi∗Expi :MATH] where n denotes the number of model genes, Coefi denotes each model gene risk coefficient, and Expi denotes each model gene mRNA level. Prognostic evaluation of ICDS index Kaplan–Meier survival curves were generated using the ‘survival’ package in R, followed by univariate and multivariate Cox regression analyses. Forest plots were created using the ‘forestplot’ package. Clinical nomograms were constructed utilizing the ‘survival’ and ‘rms’ packages. Analysis of immune cell and cytokine infiltration Utilizing the R package ‘GSVA’, ssGSEA was employed to quantify the relative abundance of various TME-infiltrating cells and cytokines within ICDS subgroups [[64]24]. Predictive analysis of immunotherapy Utilizing the R package ‘easier’, we computed immunotherapy response scores for ICDS subgroups based on cellular components, pathway activity, transcription factor activity, ligand-receptor weights, and intercellular interactions. Pathway enrichment analysis Through the R package ‘clusterProfiler’, pathways from the Gene Ontology (GO) Biological Process Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) were subjected to over-representation analysis (ORA) and gene set enrichment analysis (GSEA). Differential genes in the single-cell dataset were identified using the FindAllMarkers function (p < 0.01, avg_Log2FC ≥ 0.25), while bulk transcriptomic differential genes were analyzed with the R package ‘limma’ (p < 0.05, |logFC|< 0.5). Single-cell ICDS differential gene scoring Quantitative analysis of differential gene expression between ICDS subgroups in the TCGA-CRC cohort was conducted using the R package ‘AUCell’. Visualization of cell type-specific scores was achieved through the use of violin plots and ridge plots, produced by the R packages ‘ggplot2’ and ‘ggridge’, respectively. Cell culture The human colon cancer cell lines DLD1 and HCT116 were obtained from the American Type Culture Collection (ATCC). The HCT116 cell line was cultured in DMEM medium (Gibco, Carlsbad, USA) supplemented with 10% fetal bovine serum, while the DLD1 cell line was maintained in RPMI-1640 medium (Gibco, Carlsbad, USA) with 10% fetal bovine serum. THP-1 cells were purchased from Procell (Wuhan, China). THP-1 cells were maintained in RPMI 1640 medium (Gibco, Carlsbad, USA) augmented with 0.05 mM 2-mercaptoethanol (Sigma-Aldrich) and 10% fetal bovine serum. All cell lines were kept in a cell incubator at 37 °C with a 5% CO2 atmosphere. Macrophage induction and polarization Initially, THP-1 cells were exposed to 100 nM PMA (Sigma-Aldrich) fo 24 h to facilitate their differentiation into M0 phase macrophages. Subsequently, the M0 macrophage culture medium was supplemented with 20 ng/ml IL-4 (PeproTech) and 20 ng/ml IL-13 (PeproTech), followed by a 48-h incubation to induce polarization into M2 macrophages. Collection of cuproptosis cell fragments The HCT116 and DLD1 colon cancer cell lines were uniformly plated in six-well plates and cultured until they reached 80% confluence. Then, 0.5 µM Elesclomol (Selleck, S1052) was added, and the cells were incubated for 48 h. After 48 h, extensive cell death was observed. The cell culture medium was then subjected to centrifugation at 1000 g for 10 min, and the supernatant was collected for co-culturing with M2 macrophages. Real-time quantitative PCR (qRT-PCR) assay The reverse transcription of total RNA was performed utilizing Hifair® III Reverse Transcriptase (Yeasen Biotechnology Co., Ltd., Shanghai, China) according to the manufacturer's specifications. Subsequently, the resultant product was immediately stored at −80 °C for future use. Quantitative real-time polymerase chain reaction (qRT-PCR) was conducted using a QuantStudio 3 Real-Time PCR System (Thermo Fisher) in conjunction with Hieff® qPCR SYBR Green Master Mix (Yeasen Biotechnology Co., Ltd., Shanghai, China). The qRT-PCR protocol entailed an initial denaturation step at 95 °C for 5 min, followed by 40 amplification cycles, each comprising denaturation at 95 °C for 10 s, primer-specific annealing at 60 °C for 20 s, and extension at 72 °C for 20 s. Primer sequences are provided in Table S1. The relative expression levels of the target genes were determined using the 2 − ΔΔCt method, with GAPDH employed as the internal reference gene. Western blot Proteins were extracted from cell lines using RIPA buffer, and equivalent protein samples (30 μg, quantified by the BCA Protein Assay Kit) were separated on SDS-PAGE gels. The proteins were subsequently transferred onto PVDF membranes. Post-transfer, the membranes were blocked with TBST containing 5% BSA, followed by overnight incubation at 4 °C with primary antibodies targeting GAL3ST4 (1:3000). The membranes were then hybridized with secondary antibodies at 37 °C for 1 h, after which band intensities were detected using the ECL-plus kit. Immunofluorescence of cells Inoculate the prepared cells into a 24-well plate with coverslips, cultivating until a suitable cell density is attained; aspirate the culture medium and wash with pre-cooled PBS for 5 min, repeated three times; fix the cells with 4% paraformaldehyde, followed by washing with PBS for 5 min, repeated three times; treat with a permeabilization solution containing 0.1% Triton X-100 for 15 min, followed by washing with PBS for 5 min, repeated three times; treat with a permeabilization solution containing 0.1% Triton X-100 for 15 min, followed by washing with PBS for 5 min, repeated three times; aspirate the blocking solution and add the primary antibody (1:200), incubating overnight at 4 °C; the next day, aspirate the primary antibody, wash with PBS for 5 min, repeated three times, then add the fluorescent secondary antibody, incubating in the dark at room temperature for 1 h, followed by PBS washing for 5 min, repeated three times; use a pipette tip to apply an antifade mounting medium containing DAPI onto a glass slide, and gently place the coverslip on top; observe the cells under a fluorescence microscope. Immunofluorescence Assay for Tissue Sections Subject the paraffin sections to a deparaffinization agent (Wuhan, Seville) until they achieve complete hydration; apply an antigen retrieval buffer (Wuhan, Seville) to the prepared tissue sections to restore antigenicity, subsequently rinsing with PBS for three repetitions of 5 min; administer a 3% BSA solution onto the sections for a blocking period of 30 min; prepare a mixed solution containing two primary antibodies targeting distinct protein markers, dispense the solution onto the sections, and incubate overnight at 4 °C in a humidified chamber, followed by washing with PBS for three cycles of 5 min each; introduce the appropriate secondary antibodies, incubating for 1 h at room temperature while shielded from light, then perform PBS washes for three repetitions of 5 min each; apply DAPI staining solution to the sections and incubate at room temperature for 10 min, subsequently washing with PBS for three repetitions of 5 min; apply DAPI staining solution to the sections and incubate at room temperature for 10 min, subsequently washing with PBS for three repetitions of 5 min; utilize an antifade mounting agent for the sections, and subsequently acquire images with a fluorescence scanner (Pannoramic MIDI, 3DHISTECH). Immunohistochemistry Immerse paraffin sections in a deparaffinization solution (Wuhan, Seville) until fully rehydrated; apply an antigen retrieval solution (Wuhan, Seville) to the processed tissue sections for antigen restoration, followed by three washes with PBS for 5 min each; incubate the sections in a 3% hydrogen peroxide solution at room temperature, shielded from light, for 25 min to inhibit endogenous peroxidase activity, followed by three washes with PBS for 5 min each; incubate the sections in a 3% hydrogen peroxide solution at room temperature, shielded from light, for 25 min to inhibit endogenous peroxidase activity, followed by three washes with PBS for 5 min each; incubate the sections in a 3% hydrogen peroxide solution at room temperature, shielded from light, for 25 min to inhibit endogenous peroxidase activity, followed by three washes with PBS for 5 min each; introduce the corresponding secondary antibody and incubate at room temperature for 1 h, followed by three washes with PBS for 5 min each; dispense freshly prepared DAB chromogen, controlling the color development time under a microscope, with positive staining appearing brown, and subsequently rinse with tap water to halt the reaction; subject the sections to hematoxylin counterstaining for approximately 3 min, followed by rinsing under running water; then apply a differentiation solution for a few seconds, rinse again, and use a bluing solution followed by another rinse; sequentially immerse the sections in 75% ethanol for 5 min, 85% ethanol for 5 min, absolute ethanol I for 5 min, absolute ethanol II for 5 min, n-butanol for 5 min, and xylene I for 5 min for dehydration and clearing; then remove the sections from xylene, allow them to air dry briefly, and mount with adhesive; observe the sections under a bright-field microscope. Cell transfection Transfection is conducted when the cell density reaches 70%. siRNA is directly mixed with the transfection reagent (Advanced DNA RNA Transfection Reagent, AD600150) in a 1:1 ratio. The mixture is gently pipetted 10–15 times to ensure thorough homogenization and left to stand at room temperature for 10–15 min. The nucleic acid complexes are then added to the cell culture medium for transfection. After 48 h, the mRNA expression levels are assessed using PCR. The si-RNA (si-GAL3ST4) sequences are detailed in Table S2. Cell coculture and cell metastasis assays Thereafter, DLD1 and HCT116 cells were introduced into Transwell chambers (8-μm pore membrane, Corning, NY, USA) at a density of 2 × 10^5 cells per chamber, with varying treatments applied according to experimental requirements. Following a 48-h incubation, the chambers were fixed using 4% paraformaldehyde, and cells adhering to the upper side of the membrane were carefully removed by gentle wiping. The membranes were subsequently stained with crystal violet and examined under an optical microscope. The CCK8 assay The cells were initially seeded into 96-well culture plates at a density of 2000 cells per well. After a 24-h incubation, the cells underwent various treatments. Cell viability was evaluated at 0, 24, 48, 72, and 96 h post-treatment using the CCK8 assay, following the manufacturer's protocol. Human samples For the internal CRC sample cohort, three patients with a CRC diagnosis were recruited, and their clinical data were presented in Table S3. Tissue samples were obtained from the First Affiliated Hospital of Zhengzhou University, and all procedures were performed in accordance with the approval granted by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (Approval No: 2023-KY-0353-003). All participating patients provided informed consent. Statistical Analysis All data analyses were conducted using R software version 4.2.2. Pearson's correlation test was used for correlation analysis. Comparisons between the two groups were calculated using unpaired Student’s t-test or unpaired nonparametric Mann–Whitney test. Each assay was independently replicated more than three times, with results displayed as means ± standard error of the mean (SEM) in the figures. Multiple time point comparisons between two groups, as well as comparisons among more than two groups, were computed using analysis of variance (ANOVA). Results Construction of the immune-related cell death and stress (ICDS) prognosis model To identify the immune-related cell death and stress (ICDS) genes, we initially conducted enrichment scoring of immune cells, programmed cell death (PCD) and stress pathways using the Gene Set Variation Analysis (GSVA) algorithm. Subsequently, the Consensus Clustering method was used to distinguish between high-expression and low-expression groups (Fig. S1A-F). Following this, the Weighted Correlation Network Analysis (WGCNA) (Fig. S1G-J) was utilized to screen for cell death and stress related genes, while Linear Models for Microarray Data (Limma) was applied to identify immune-related genes (Fig. S1K). This analysis revealed a total of 4,295 immune-related genes and 2,769 cell death and stress related genes. By filtering for genes common to both categories, we identified 579 the ICDS genes (Fig. [65]1A). To further screen clinically relevant ICDS genes, we performed univariate Cox analysis and Lasso regression analysis (Fig. [66]1B-C), resulting in the identification of 29 ICDS genes. Then, based on Stepwise Cox Proportional Hazards Regression, the ICDS model composed of 9 key genes was finally determined, and the ICDS index was calculated (Fig. [67]1D). Univariate Cox analysis of overall survival time in the TCGA-CRC cohort indicated that CXXC5, ATP6V1B2, FAS, and MMP1 were associated with improved survival for colorectal cancer patients. Conversely, BST2, TRIB2, ENO2, GAL3ST4, and FLT1 were identified as adverse factors affecting overall survival (Fig. [68]1E). Fig. 1. [69]Fig. 1 [70]Open in a new tab Construction of the immune-related cell death and stress (ICDS) prognostic model and its clinical features. A Gene selection process for the ICDS model. B LASSO coefficient curve for 29 predictor variables. C Identification of the optimal lambda penalty in the LASSO regression model. D Nine model genes and their coefficients selected through LASSO Cox regression analysis. E Forest plot showing the univariate Cox values of the model genes. Red dots indicate negative effects, blue dots indicate positive effects. The points represent estimates, and the error bars indicate 95% confidence intervals. F–H Boxplots showing the association between the ICDS index and variations in clinical indicators within the TCGA-CRC cohort. I Display the differential expression of signature genes across ICDS groups in the TCGA-CRC cohort, with the ICDS index and clinical factors presented as figure legends. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 The clinical staging results from the TCGA-CRC cohort demonstrated a significant correlation between the ICDS index and the T, N and the clinical stage. Specifically, in the T stage, patients in the T3 and T4 categories exhibited higher ICDS index (Fig. [71]1F). In the N stage, individuals with N1 and N2 classifications also showed elevated ICDS index (Fig. [72]1G). Furthermore, in the clinical stage, patients in stages III and IV had higher ICDS index compared to those in stages I and II (Fig. [73]1H), which were linked to poorer clinical outcomes. The relationship between these clinical traits and the ICDS index had also been corroborated by other datasets (Fig. S2B-E). Among these cohorts, no significant correlation between the ICDS index and the M stage was detected (Fig. S2A). A further analysis was conducted to investigate the differential expression of modeling genes between the high and low ICDS index groups within the TCGA-CRC cohort. The results revealed that, except for CXXC5, all other genes exhibited significant differences in expression (F[74]ig. [75]1I). Moreover, our findings reveal correlations between the ICDS index and model genes with T stage, N stage, and clinical stage, suggesting their potential utility in clinical prognosis. (F[76]ig. [77]1I and Fig. S2F-G). Prognostic value of the ICDS index in multiple colorectal cancer cohorts To assess the applicability of the ICDS index in predicting the survival probability of CRC patients, we categorized patients into high ICDS index and low ICDS index subgroups based on quartiles across seven datasets. Kaplan–Meier survival analysis revealed that CRC patients in the high ICDS index subgroup exhibited a significantly elevated risk of death across all cohorts (p < 0.05) (Fig, 2A-F right and Fig. S2H). The comparison of receiver operating characteristic (ROC) curves revealed that the ICDS index showed strong predictive accuracy for overall survival in CRC cohorts at 1, 3, and 5 years. (Fig, 2A-F left and Fig. S2I). We performed both univariate and multivariate Cox analyses to thoroughly evaluate the prognostic significance of the ICDS index within the TCGA-CRC cohort. The results of both indicate that the ICDS index is an independent risk factor in the TCGA-CRC cohort (p < 0.05). (Fig. [78]2G and Fig. S2J). In addition, we constructed a predictive nomogram model that integrates multiple clinical factors, including the ICDS index, to predict overall survival (OS) at 1, 3, and 5 years (Fig. [79]2H). The C-index of the nomogram models was 0.820 (95% CI 0.70–0.872). The bootstrap-corrected calibration curve showed that the nomogram model for overall survival exhibits good predictive efficiency (Fig. S2K). In summary, our findings indicated that the ICDS index could robustly predict the prognosis of patients with CRC. Figure2. [80]Figure2 [81]Open in a new tab Prognostic validation of the ICDS model. A-F ROC curves (1-, 3-, and 5-year) and overall survival Kaplan–Meier curves for high and low ICDS index groups across datasets (TCGA-CRC, [82]GSE29621, [83]GSE38832, [84]GSE39582, [85]GSE17536, [86]GSE17537). G Forest plot showing the univariate Cox regression hazard ratio (HR) and 95% confidence interval (CI) for overall survival, based on clinical characteristics and ICDS index. H Nomogram developed to predict overall survival probability based on clinical characteristics and ICDS index. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 Immune cell infiltration, programmed cell death and cell stress characteristics of the ICDS index We evaluated the impact of the ICDS index on immune cell infiltration based on our model features. In comparison to the high-ICDS index group, the low-ICDS index group exhibited increased levels of activated CD4 + T cells, activated dendritic cells (DCs), central memory CD8 T cells, effector memory CD4 T cells, cytotoxic T cells, macrophages, mast cells, myeloid-derived suppressor cells (MDSCs), as well as helper T cells of types 1, 2, and 17 (Fig. [87]3A). For these cells exhibiting significant differences, aside from MDSCs, which might play a role in immunosuppression, the remaining cells primarily contributed to immune antigen presentation, phagocytosis, or direct tumor cell killing, while also promoting immunity through the secretion of various cytokines. This phenomenon might suggest that the low-ICDS index subgroup exhibited a highly activated immune state, which could be associated with its favorable prognosis. We further analyzed the cytokine level differences between the two subgroups, and the results showed significantly elevated expression levels of pro-inflammatory cytokines, including IL1B, IL6, IL18, TNF, OSM, CSF2, CCL2, CCL3, CCL4, and IL10, in the subgroup with a lower ICDS index (Fig. [88]3B). These findings further indicated that the subgroup with a lower ICDS index displayed a heightened immune-promoting state, consistent with the observed differences in immune cell profiles. In the comparison of immune checkpoints, ICOSLG (CD275) and SIGLEC15 (CD33) were found to be up-regulated in the high ICDS index group. Conversely, the low ICDS index group demonstrated high expression of immune-suppressive checkpoints such as CD274, CD70, CTLA4, PDCD1LG2 (PD-L2), and IDO1, as well as immune-promoting checkpoints including ICOS and TNFRSF9 (CD137) (Fig. [89]3C). These findings indicated that the two patient groups had disparate responses to immunotherapy, with the low-ICDS index group likely exhibiting superior therapeutic effectiveness. This conclusion was corroborated by immunotherapy predictions derived from the Easier package (Fig. [90]3D). According to the Estimate algorithm, patients with low ICDS index exhibited higher immune scores and lower tumor purity scores (Fig. [91]3E-F), suggesting that low-ICDS index patients have greater immune infiltration, which aligns with the scoring of immune cells and cytokines. Interestingly, the stromal scores displayed no significant differences among the various ICDS index subgroups (Fig. [92]3G), suggesting that variations between ICDS subgroups were largely linked to immune cell differences. To investigate the functional differences between the high and low-ICDS index groups, we conducted GO functional pathway enrichment analysis on differential genes. The low ICDS index group was mainly enriched in immune-related pathways, while the high ICDS index group was predominantly enriched in pathways associated with stromal growth (Fig. [93]3H). The gene mutation waterfall plot revealed that the most common mutated genes in both subgroups were APC, TP53, TTN, and KRAS (F[94]ig. [95]3I), with only TP53 showing a significant difference (Fig. S3A). Meanwhile, genomic changes (FGA), genomic amplifications (FGG), and genomic deletions (FGL) (Fig. S3B-D), along with tumor mutation burden (TMB), single nucleotide polymorphisms (SNPs), and insertion and deletion mutations (INDELs) between the two groups showed no significant differences (Fig. S3E-G). Genomic findings indicated that the immune disparities between the two patient groups did not correlate with tumor mutational burden, suggesting that changes in adjuvant-like factors were the primary drivers of the observed differences in immune infiltration. Fig. 3. [96]Fig. 3 [97]Open in a new tab Characteristics and differences between the ICDS subgroups. A-C Differences between the ICDS subgroups in 28 immune and stromal cells (A), immune factors (B), and immune checkpoints (C). D Boxplot showing the differences in immunotherapy response scores (Easier Score) between the ICDS subgroups. E–G Boxplots displaying differences in Immune Score (E), Tumor Purity (F), and Stromal Score (G) between ICDS subgroups. H Gene Ontology (GO) pathway enrichment analysis showing significant differences in differential genes between the ICDS subgroups. I Tumor mutation burden (TMB) distribution and mutation waterfall plot for the ICDS subgroups. J-K Gene set enrichment analysis (GSEA) of cell death (J) and cell stress (K) showing relative enrichment differences between the ICDS subgroups. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 Further investigation into the differences in cell death and stress between the two groups using GSEA analysis showed that the low-ICDS index group exhibited significant enrichment in multiple cell death pathways, especially cuproptosis (Fig. [98]3J), along with endoplasmic reticulum stress and oxidative stress (Fig. [99]3K). In conclusion, these findings validated that ICDS subgroups exhibited marked differences in immune cell infiltration, programmed cell death and cell stress. Single-cell analysis revealed the key differential cell subpopulations and regulatory genes of the ICDS model To investigate the critical cellular subpopulations contributing to ICDS subgroup differences, we performed single-cell analysis. After stringent quality control and filtering of the [100]GSE132465 single-cell dataset, we retained cells originating from tumor tissues, resulting in 27,292 single cells. Using unsupervised graph-based clustering on all cells and identifying them based on classical cell markers, we identified six main cell types: B cell, Endothelial cell, Epithelial cell, Fibroblast, T/NK cell, and Myeloid cell (Fig. [101]4A-B). Then, AUCell was used to calculate the enrichment scores of ICDS subgroup differential genes across each cell type, revealing that myeloid cells were highly enriched with ICDS subgroup differential genes (Fig. [102]4C-D). To further identify the key subpopulations contributing to the differences, we conducted additional unsupervised graph-based clustering on myeloid cells, classifying them into three cell types (Dendritic cell, Macrophage, Monocyte) using classical markers (Fig. [103]4E and Fig. S4A). Subsequent AUCell calculations revealed that the differential gene scores were highest in macrophages (Fig. [104]4F). Macrophages are broadly classified into pro-inflammatory M1-type and anti-inflammatory M2-type based on their functions. We classified the three macrophage clusters (Fig. [105]4G-H) as M1-like and M2-like types based on their functional roles. Cluster 1 was primarily enriched in functions related to synapse modification, lipoprotein catabolism, collagen degradation, and regulation of lipid transport, which are associated with wound healing and tissue repair (Fig. S4B). Thus, we defined this cluster as M2-like macrophages. Conversely, Cluster 2 (Fig. S4C) and Cluster 3 (Fig. S4D) were enriched in functions related to leukocyte aggregation, neutrophil chemotaxis, inflammatory response and MHC-II complex assembly, associated with pro-inflammatory and antigen presentation roles. Therefore, we defined these cells as M1-like macrophages. Interestingly, we found significant differences between M1-like and M2-like macrophages when we recalculated enrichment scores using differential genes (F[106]ig. [107]4I). These findings suggest that variations in macrophage types could play a major role in driving the immune differences observed between the ICDS subgroups. Fig. 4. [108]Fig. 4 [109]Open in a new tab Single-cell analysis revealed the key differential cell subpopulations and regulatory genes of ICDS. A The t-SNE plot of unsupervised clustering of cells from the single-cell dataset. B Heatmap of characteristic genes for the cell types. The dot size reflects the percentage of cells expressing these genes, and dot color intensity indicates the average expression levels of the genes. C Ridge plot displaying differences in Area Under the Curve (AUCell) scores for differential genes derived from the ICDS subgroup across cell types. D The t-SNE plot reveals differences in the density distribution of enrichment scores for differential genes from the ICDS subgroup across cells. E Unsupervised clustering of myeloid cells, represented as a t-SNE plot. F The box plot illustrates differences in AUCell scores for the ICDS differential genes among three myeloid cell types. G-H Unsupervised clustering of macrophages resulted in three clusters (G), which were further categorized into two functional subtypes (H), as shown in the t-SNE plot. I The box plot illustrates differences in AUCell scores for ICDS differential genes between two macrophage populations. J The heatmap shows the relative expression values of model genes across different cell subpopulations. K-L The t-SNE plots display the differences in density distribution of GAL3ST4 in myeloid cells (K) and macrophages (L). M Gene preferences of different