Abstract Introduction This study leverages bioinformatics and medical big data to integrate datasets from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA), providing a comprehensive overview of immunogenic cell death (ICD)-related gene expression in colorectal cancer (CRC). The research aims to elucidate the molecular pathways and gene networks associated with ICD in CRC, with a focus on the therapeutic potential of cell death inducers, including ferroptosis agents, and their implications for precision medicine. Methods We conducted differential expression analysis and utilized advanced bioinformatic techniques to analyze ICD-related gene expression in CRC tissues. Unsupervised consensus clustering was applied to categorize CRC patients into distinct ICD-associated subtypes, followed by an in-depth immune microenvironment analysis and single-cell RNA sequencing to investigate immune responses and cell infiltration patterns. Experimental validation was performed to assess the impact of cell death inducers on ICD gene expression and their interaction with ferroptosis inducers in combination with other clinical drugs. Results Distinct ICD gene expression profiles were identified in CRC tissues, revealing molecular pathways and intricate gene networks. Unsupervised consensus clustering refined the CRC cohort into unique ICD-associated subtypes, each characterized by distinct clinical and immunological features. Immune microenvironment analysis and single-cell RNA sequencing revealed significant variations in immune responses and cell infiltration patterns across these subtypes. Experimental validation confirmed that cell death inducers directly affect ICD gene expression, highlighting their therapeutic potential. Additionally, combinatorial therapies with ferroptosis inducers and clinical drugs were shown to influence drug sensitivity and resistance in CRC. Discussion Our findings underscore the importance of ICD-related genes in CRC prognosis and therapeutic targeting. The study provides actionable insights into the efficacy of cell death-inducing therapies, particularly ferroptosis inducers, and their regulatory mechanisms in CRC. These discoveries support the development of precision medicine strategies targeting ICD genes and offer valuable guidance for translating these therapies into clinical practice, with the potential to enhance CRC treatment outcomes and patient survival. Keywords: immunogenic cell death (ICD), colorectal cancer (CRC), ferroptosis inducers, combinatorial therapies, DAMPs, damage-associated molecular patterns Graphical Abstract [33]graphic file with name fimmu-15-1458270-g010.jpg Introduction Unlike proinflammatory necrosis, apoptosis has traditionally been regarded as a form of acceptable cell death. It is currently established that damage-associated molecular patterns (DAMPs), which are released after the apoptosis of cancerous or diseased cells, can trigger both innate and adaptive immune responses. Immunogenic cell death (ICD), which is triggered by the release of DAMPs, can result in the generation of tumor-specific CD8+ T cells and immunological memory ([34]1). ICD seems like a potential way to treat cancer, but as of right now, only a small number of tumor therapeutics—including radiation, oncolytic virotherapy, oxaliplatin, cetuximab, bortezomib, photodynamic therapy, and extracorporeal photochemotherapy—are known to cause ICD. These medications encourage alterations in the extracellular milieu or cell surface that are consistent with the induction of ICD ([35]2, [36]3). The avoidance of immunological destruction, which can be accomplished by secreting immunosuppressive substances, enlisting immunosuppressive cells as T regulatory cells, and downregulating the activation of cytotoxic lymphocytes, is one of the primary characteristics of cancer as outlined by Hanahan and Weinberg in 2011 ([37]4). In this context, the natural history of a number of solid malignancies, such as melanoma, lung cancer, and urological tumors, has changed due to the use of immune checkpoint inhibitors (ICIs) ([38]5, [39]6). Antigen-presenting cells (APCs) and cancer cells, in fact, express ligands on their plasma membranes that bind to immunological checkpoints on lymphocytes and prevent their activation; ICIs disrupt this process and allow T-cell function to be restored ([40]6). The relationship between ICIs and ICD inducers and their possible function in inducing long-lasting immune responses is poorly understood. This could ultimately lead to cancer vaccination, particularly in diseases like gastrointestinal (GI) cancer where the benefits of immunotherapy alone are not very great ([41]7). With 3.4 million cancer deaths from the disease in 2018 and nearly 26% of all cancer cases worldwide, gastrointestinal cancer is a serious problem ([42]8). Zheng’s study highlights the ways in which ferroptotic colon cancer cells and antitumor immunity work in concert by boosting immunogenicity, specifically the production of HMGB1 and the exposure of CRT ([43]9, [44]10). In addition to accelerating DC maturation, these cells also enhance the immunological milieu by raising the frequency of antitumor immune cells and T cell release of IFN-γ ([45]11). Crucially, immune checkpoint inhibitors (ICIs) become much more effective when ferroptosis is induced, and this improvement is reliant on antitumor immunity being activated ([46]12). The study also discovered that while certain cells, including MDSCs, are resistant to ferroptosis, neutral ceramidase N-acylsphingosine amidohydrolase can be inhibited to overcome this resistance ([47]13). Immune cells, such as CD8+ T cells, collaborate with arachidonic acid to induce iron poisoning in cancer cells ([48]14). These results demonstrate the potential of ferroptosis in colon cancer immunotherapy and offer novel therapeutic approaches for future ferroptosis-targeted immunotherapy combinations ([49]14). The identification of treatment-related targets and the investigation of the activating effects of ferroptosis inducers on immunogenic cell death and antitumor immunity are necessary because the therapeutic mechanism of targeting immunogenic cell death-related genes in tumors through the ferroptosis mechanism of tumor cells remains unknown ([50]11). Materials and methods Data download and processing RNA-seq data, related clinical data, and survival data for colorectal cancer were downloaded from the UCSC-Xena platform. Included in this were count numbers from the TCGA-COAD and TCGA-READ data sets, which the TCGAbiolinks package was used to standardize and process further ([51]15). To be more precise, the sample bearing the ID “-11A” was classified as neighboring normal tissue, whilst the remaining samples were identified as malignant tissue. Of them, 689 samples containing prognostic data—51 normal tissue samples and 638 malignant tissue samples—were used in the creation and examination of the models that followed. Gene-expression data processing Initially, we concentrated on creating a list of genes that were closely linked to immunogenic cell death (ICD) in order to better understand the molecular mechanisms underlying this phenomenon. These genes, which have been chosen from previous research, serve as the cornerstone of our ICD investigation. We then successfully created a gene expression matrix specifically for ICD by using differential analysis techniques to match these genes with expression profiling data from tumor and normal tissues. The primary dataset for our ensuing analyses is this matrix. The exact roles of ICD-related genes in tumor development were determined by using the Wilcoxon rank-sum test to the ICD gene expression data from tumor and normal tissue samples in the TCGA database. The study’s findings demonstrated that there were differences in the patterns of these genes’ expression between the two sets of samples. Furthermore, Pearson correlation analysis was used to explore potential interactions between these genes, revealing expression connections that improved our understanding of the ICD gene network. We used the STRINGdb program to build a protein-protein interaction (PPI) network based on known protein interactions for ICD-related genes in order to reveal the molecular processes of ICD at the protein level ([52]16). This network offers insights into the intricate processes involved in ICD by visualizing important proteins and their connections. Lastly, we used the clusterProfiler tool to conduct pathway enrichment analysis in order to obtain a better understanding of the possible roles that these genes may play in biological processes. Significant enrichments of ICD-related genes were found in particular biological pathways by this study, suggesting that these genes may play roles in immune response, tumor growth, and other biological processes. These thorough findings set a strong basis for further study in this area and offer vital insights on the function of ICD in tumor development. Unsupervised clustering analysis identified ICDRGs expression patterns To gain a deeper understanding of the molecular characteristics and subtype prediction of CRC, we employed an approach focused on the differential expression of ICD-related genes. Utilizing the “ConsensusClusterPlus” package (version 1.62.0), we performed unsupervised clustering analysis on CRC patients, setting the parameters maxK to “6”, clusterAlg to “hc” (hierarchical clustering), and distance to “Pearson”. During this process ([53]17), we identified the optimal number of immune subtypes (K) by tracking cells, analyzing cumulative distribution functions, and their relative changes. Subsequently, we collected and analyzed sample information from different clinical subgroups (such as age, stage, TNM classification, and survival status) within each subtype. To visually present this information, we leveraged the tidyverse (version 1.3.2) and ggplot2 (version 3.3.6) packages to visualize the classification information of different subtypes across various subgroups ([54]18, [55]19). In the subtype analysis phase, we utilized the TCGA dataset and employed the Kaplan–Meier curve method to assess the correlation between ICD-related gene consensus clustering subtypes and survival outcomes, based on subtype classification, patient survival time, and status. Additionally, we generated a heatmap of gene expression, revealing the relationship between gene expression, subtype classification, and clinical characteristics by observing the expression values of genes in TCGA samples and the clinical features of each sample. This integrated analysis not only provided us with a profound understanding of CRC molecular subtypes but also highlighted the potential value of ICD-related genes in predicting patient prognosis and guiding treatment. These findings offer new strategies and insights for precision medicine in CRC, paving the way for further research and clinical applications. Analysis of immune microenvironment differences between subtypes In this study, we employed two advanced bioinformatics algorithms to extensively profile the immune microenvironment of colorectal cancer samples from The Cancer Genome Atlas (TCGA). The application of these methods not only sheds light on the intricate mechanisms of tumor immune response but also provides crucial scientific evidence for future clinical diagnosis and therapeutic strategies. Firstly, we leveraged the CIBERSORT algorithm ([56]http://cibersort.stanford.edu/index.php), a robust tool based on linear support vector machine (SVM) principles, to accurately deconvolve the expression matrix of immune cell subtypes in tumor samples ([57]20). Through both relative and absolute modes of computation, we determined the composition of 22 immune cell types within the samples. The relative mode ensures that the proportions of different immune cell subtypes sum to one, reflecting their relative contributions to the overall immune response. Meanwhile, the absolute mode directly provides the absolute counts of each immune cell type, offering a more intuitive understanding of their actual roles in immune responses. Furthermore, we utilized the Wilcoxon rank-sum test to analyze the differences in the proportions of various immune cell subtypes and calculated the corresponding statistical significance (p-values). Additionally, we plotted boxplots based on the subgrouping of immune cell types to visually demonstrate these differences. Secondly, to assess the enrichment of different immune gene sets in tumor samples, we employed the ssGSEA (2.0) algorithm to evaluate 17 immune gene sets from the import database ([58]21). This approach comprehensively reflects the overall expression levels of gene sets in samples, thus uncovering their potential roles in immune responses. Specifically, we focused on the expression differences of HLA family genes among different immune cell subtypes, as these genes play a pivotal role in immune responses. By extracting the expression data of HLA family genes and combining it with the Wilcoxon rank-sum test, we compared the expression differences of immune checkpoint genes and HLA family genes between different immune cell subtypes and calculated the statistical significance (p-values) between high-risk and low-risk groups. Finally, we utilized the ggpubr package in R to plot boxplots, graphically representing these differences ([59]22). In summary, by integrating the CIBERSORT algorithm and the ssGSEA algorithm, we comprehensively and in-depth profiled the immune microenvironment of TCGA colorectal cancer samples. These findings not only enhance our understanding of the complex mechanisms of tumor immune response but also provide robust data support and scientific evidence for future clinical diagnosis and therapeutic strategies. Functional enrichment analysis of genes Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment studies were conducted using several R packages, such as “org.Hs.Eg.Db,” “ggplot2,” and “clusterProfiler,” to investigate the biological processes enriched by the genes in key modules ([60]23). Any GO or KEGG terms that showed significant enrichment were examined further with a p-value threshold of 0.05. Establishing a prognosis signature related to ICDGs Further, utilizing the limma package in R language, which employs linear regression and empirical Bayes methods, we conducted differential expression analysis comparing tumor samples with normal samples across all TCGA colorectal cancer datasets. This analysis yielded gene-specific information such as P-values and logFC. Additionally, we applied the Benjamini & Hochberg method for multiple testing correction, resulting in adjusted P-values (adj.P.Value) ([61]24). We evaluated the genes based on both fold-change and significance, setting thresholds as adj.P.Value < 0.05 and |logFC| > 1. Incorporating the clinical survival and prognostic data from the TCGA colorectal cancer samples, we screened for genes significantly associated with overall survival prognosis using univariate Cox regression analysis from the survival package in R. Genes with a P-value less than 0.05 were considered significantly correlated. LASSO Cox analysis identified the genes most relevant to overall survival. To prevent overfitting, we performed 10-fold cross-validation using the glmnet package in R. Based on RNA expression levels, we then calculated a risk score for each patient using the formula: [MATH: Risk sc ore= ∑i=1nc oefi X < /mo>id :MATH] (1) where coefi represents the coefficient and Xi is the normalized count for each gene. Using clinical follow-up data from the [62]GSE17536 and [63]GSE39582 colorectal cancer datasets, we computed the Risk score for each sample in both the TCGA training set and the GEO validation set. Samples were then divided into High (Risk score above the median) and Low (Risk score below or equal to the median) groups. The Kaplan-Meier curve method from the survival package was employed to assess the correlation between these groups and actual survival outcomes. This risk model was further validated for robustness by evaluating risk scores on an external independent dataset. Quality control and cell-type identification For quality control, Seurat (version 4.3.0) was used to count unique molecular identifiers (UMIs) and mitochondrial genes. Cells with more than 100 UMIs and less than 15% mitochondrion-derived UMI counts were selected ([64]25). This study selected the top 20 components and first 2000 variable genes. The “ScaleData” function was used to regress the inflow of UMIs and the percentage of mitochondrion-derived UMI counts. Subsequently, the main cell clusters were identified by Seurat’s “FindClusters” function. Unbiased cell type recognition was visualized by umap. Cell clusters annotation To delineate the specific marker genes for individual cell clusters, we employed the ‘FindMarkers’ function in Seurat to contrast the gene expression profiles of cells within a given cluster against all other cells in the dataset. This function utilizes a Wilcoxon rank-sum test to identify differentially expressed genes (DEGs) between the two groups. Subsequently, the obtained P values were adjusted for multiple testing using Bonferroni correction, considering the total number of genes analyzed. Marker genes were defined as those with an adjusted P value below 0.05 and exhibiting at least a twofold higher average expression level in the cluster of interest compared to all other clusters. To validate the annotation of each cell cluster, we relied on the expression of canonical marker genes. Specifically, T/NK cells were identified by the expression of CD3D, CD3E, and NKG7; memory B cells were characterized by CD79A and MS4A1; plasma cells were defined by the presence of IGHG1, JCHAIN, and MZB1. Monocytes and macrophages were distinguished by CD68, CD163, CD14, and LYZ, while dendritic cells were identified through CD74, CLEC9A, and CD1C. Fibroblasts exhibited high expression of COL1A1, ACTA2, and TAGLN, while endothelial cells were marked by VWF, PLVAP, and CLDN5. Epithelial cells, on the other hand, were distinguished by EPCAM, KRT8, and KRT18. This comprehensive analysis, leveraging both statistical testing and the expression of established marker genes, allowed us to accurately annotate and characterize each cell cluster. Chemotherapy response analysis Information on drug responsiveness and drug targeting pathways was gathered using the Cancer Drug Sensitivity Genomics (GDSC) platform. Subsequently, the pRRophetic package in R was employed to predict the drug sensitivity of various phenotypes from gene expression data, and to obtain the pharmaceuticals associated with different classes after tabulating the sensitivity data of the two GBM classes to various drugs, thereby setting the stage for the selection of clinical drugs ([65]26). Quantitative reverse transcriptase PCR The HCT-116(human colon cancer cell) lines were procured from ATCC. Following incubation for a duration of 12 to 24 hours, these cells were subjected to treatment with various compounds, including 2uM RSL3 (Selleck, catalog number #S8155), 0.25 μM disulfiram (Selleck, catalog number #S1680), and 20uM oxaliplatin (Selleck, catalog number S1224), for specified durations. Additionally, control groups were maintained without any drug exposure to serve as references. This experimental