Skip to content
BY 4.0 license Open Access Published by De Gruyter Open Access November 9, 2023

Construction of immunogenic cell death-related molecular subtypes and prognostic signature in colorectal cancer

  • Chun Yu , Weixuan Yang , Li Tian , Yue Qin , Yaoyao Gong EMAIL logo and Wenfang Cheng EMAIL logo
From the journal Open Medicine

Abstract

Immunotherapy is a promising treatment for advanced colorectal cancers (CRCs). However, immunotherapy resistance remains a common problem. Immunogenic cell death (ICD), a form of regulated cell death, induces adaptive immunity, thereby enhancing anti-tumor immunity. Research increasingly suggests that inducing ICD is a promising avenue for cancer immunotherapy and identifying ICD-related biomarkers for CRCs would create a new direction for targeted therapies. Thus, this study used bioinformatics to address these questions and create a prognostic signature, aiming to improve individualized CRC treatment. We identified two ICD -related molecular subtypes of CRCs. The high subtype showed pronounced immune cell infiltration, high immune activity, and high expression of human leukocyte antigen and immune checkpoints genes. Subsequently, we constructed and validated a prognostic signature comprising six genes (CD1A, TSLP, CD36, TIMP1, MC1R, and NRG1) using random survival forest analyses. Further analysis using this prediction model indicated that patients with CRCs in the low-risk group exhibited favorable clinical outcomes and better immunotherapy responses than those in the high-risk group. Our findings provide novel insights into determining the prognosis and design of personalized immunotherapeutic strategies for patients with CRCs.

Graphical abstract

1 Introduction

Colorectal cancer (CRC) is currently the third most commonly diagnosed and the second most deadly cancer worldwide [1]. According to the Global Cancer Observatory, approximately 1.93 million (10%) new cases and 0.94 million (9.4%) cancer deaths due to CRC were recorded in 2020 [1]. The most common pathologic staging system for CRC is the 8th edition of the American Joint Committee Cancer Tumor Node Metastasis system, which ranges from stage I to IV, with stage 0 being the earliest stage of CRC [2]. Surgery and postoperative chemotherapy/radiotherapy have improved the prognosis of patients with early-stage CRC; however, treating advanced CRC remains a challenge [3,4]. The mechanisms underlying the development and progression of CRC are complex and have not been fully elucidated. Known factors include patient-intrinsic factors (e.g., age, genetics, and the microbiome), environmental factors (e.g., diet and lifestyle, infections and chronic inflammation, tobacco, alcohol, and pollution), and cancer cell-intrinsic mechanisms (e.g., genetic alterations, oncogenic cell signaling, and the epithelial–mesenchymal transition) [1,5]. Immunotherapy is a valuable treatment option for those with late-stage CRC. For example, immune checkpoint inhibitors received regulatory approval in 2017 for treating mismatch-repair-deficient (dMMR) or microsatellite instability-high metastatic CRC [6,7,8]. However, immunotherapy resistance remains common. Therefore, identifying biomarkers indicative of the immunotherapy response and prognosis of patients with CRC is vital.

According to the 2018 Nomenclature Committee on Cell Death classifications, immunogenic cell death (ICD) is a particular form of regulated cell death that can trigger an adaptive immune response [9,10]. Cells undergoing ICD are characterized by the release of damage-associated molecular patterns, including heat shock proteins 70 and 90, calreticulin, high-mobility group protein B1, and ATP [11,12]. Mounting studies have shown that treatment-driven ICD enhances the treatment effects of conventional radiotherapy and chemotherapy by inducing effective immune responses [13,14]. Moreover, studies have found that ICD combined with immune checkpoint inhibitors improves the anti-tumor immune effect [15,16]. Since ICD is inducible with cytotoxic drugs, it offers a potential approach for cancer immunotherapy [17].

In recent years, many preclinical studies have investigated the molecular mechanisms of ICD, but few studies have evaluated the possibility of ICD in a clinical context. Therefore, this study recognized ICD-associated subtypes of CRC and differences between subtypes in terms of prognosis, immune landscape, and somatic mutation. And then established a prognostic signature based on ICD-related gene expression. Overall, we aimed to improve immunotherapy response predictions for patients with CRC and provide new strategies for individualized therapy selection.

2 Materials and methods

2.1 Data source

The Cancer Genome Atlas (TCGA)-colon adenocarcinoma (COAD) and TCGA-rectal adenocarcinoma (READ) RNA sequencing data and matching clinical information were downloaded from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/) to create the training set. Principal component analysis of the rectal and colon samples revealed no significant differences, and merging the samples did not require adjustments [18]. Thus, the COAD and READ samples were combined into a TCGA-CRC cohort. The GSE17538 dataset was downloaded from the Gene Expression Omnibus (GEO) database as a validation set [19]. Additionally, we obtained the list of ICD-related genes from the GeneCards database (http://www.genecards.org/).

2.2 Consensus clustering

We performed a cluster analysis to identify ICD-related molecule subtypes using the “ConsensusClusterPlus” package in R software (R Core Team, Vienna, Austria), as previously reported [20]. Subsequently, to ensure the stability of the results, we examined an optimal cluster size between k = 2–9 and repeated the analysis 1,000 times. The cluster map was created using the “pheatmap” tool in R.

2.3 Differentially expressed (DE) gene and functional enrichment analyses

DE genes were identified utilizing the “Limma” package in R, and the false discovery rate (FDR) was used to correct for false-positive results. The filter criteria were FDR < 0.05 and |logFoldchange (FC)| > 1. The DE genes were then subjected to Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses using the R package “clusterProfiler” [21]. GO terms and KEGG pathways were considered significantly enriched when the corrected p-value (i.e., the q-value) was <0.05.

2.4 Gene set enrichment analysis (GSEA)

We conducted a GSEA to assess functional variations between ICD-high and ICD -low clusters; c2kegg and c5go were used as the reference gene sets. The screening conditions were p-values of <0.05, FDR q-values of  <0.1, and |normalized enrichment scores (NESs)| of >1.5.

2.5 Somatic mutation analysis

We downloaded the COAD and READ somatic mutation data from the GDC data portal. The “Maftools” package in R was used to create waterfall plots and visualize the mutated genes, as previously reported [22].

2.6 Immune microenvironment differences between ICD high and low cohorts

The estimate score, immune score, stromal score, and tumor purity were calculated for each CRC sample using the ESTIMATE algorithm through the “estimate” package in R [18]. We then utilized the CIBERSORT algorithm to evaluate the relative proportions of 22 kinds of infiltrating immune cells in each sample, as previously reported [23]. Subsequently, we compared the expression levels of the immune checkpoint and human leukocyte antigen (HLA) genes between the two clusters.

2.7 Establishing and validating an ICD-related prognostic risk signature

Univariate Cox regression analysis and random survival forest-variable hunting algorithm (a powerful ensemble algorithm based on machine learning) were used for gene selection and model construction. First, a univariate Cox analysis was used to identify prognosis-related ICD genes. Then, using the R package “randomForestSRC,” we set the parameter ntree at 1,000 to predict the significant ICD genes from preliminary screened candidates. Then, based on the multivariate Cox analysis, Kaplan–Meier (KM) tests were performed for several gene combination signatures, and p-values calculated by KM analyses were sorted for comparisons and determining the optimal combination. Receiver operating characteristic (ROC) curves were used to assess the prognostic accuracy of the risk signature. The independent GSE17538 data set was used to verify the model’s stability. Finally, after identifying independent prognostic factors of CRC, we constructed a nomogram by using the “rms” package in R and evaluated its performance.

2.8 Immunotherapy response predictions

The tumor immune dysfunction and exclusion (TIDE) algorithm is a computing architecture that integrates data on two tumor immune escape mechanisms; we used this algorithm to predict the response to immune checkpoint therapy [24]. The TIDE database (http://tide.dfci.harvard.edu/) was used to calculate the TIDE score for each TCGA-CRC sample; then, we compared the relevance between the risk score and immunotherapy response.

2.9 Statistical analysis

R v4.1.2 (https://cran.r-project.org/) or GraphPad Prism software was used to analyze data. Continuous variables were reported as means and standard deviations and compared using a Student’s t-tests with p-values. Data were visualized using R v4.1.2, GraphPad Prism, SangerBox (http://vip.sangerbox.com), Figdraw (https://www.figdraw.com), and Hiplot (https://hiplot.com.cn) software.

3 Results

3.1 Differential ICD-associated genes and biological functions

Figure 1 presents a workflow overview. First, 859 ICD-related genes were obtained from the GeneCards database; they were screened based on protein-coding genes and relevance scores larger than the median. Then, we identified 298 DE-ICD genes; 154 were up-regulated, and 144 were down-regulated based on a differential analysis between tumor and normal samples in TCGA-CRC (Figure 2a and b). GO and KEGG analyses indicated that the DE-ICD genes were primarily enriched in the activation and regulation of inflammatory cells, activities associated with immunity, and cancer-related biological processes and signaling pathways (Figure 2c and d). To further understand the connections among these DE-ICD genes, we constructed a protein-protein interaction network using the STRING database (https://string-db.org), which was visualized usingCytoscape (Figure A1). Figure 2e–g presents the top ten hub genes and the two most significant modules identified using the CytoHubba and MCODE Cytoscape plugins.

Figure 1 
                  Study workflow diagram. CRC: colorectal cancer; DE: differentially expressed; GO: Gene Ontology; GSEA: gene set enrichment analysis; HLA: human leukocyte antigen; ICD: immunogenic cell death; KEGG: Kyoto Encyclopedia of Genes and Genomes; KM: Kaplan–Meier; PPI: protein–protein interaction; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas; TIDE: tumor immune dysfunction and exclusion.
Figure 1

Study workflow diagram. CRC: colorectal cancer; DE: differentially expressed; GO: Gene Ontology; GSEA: gene set enrichment analysis; HLA: human leukocyte antigen; ICD: immunogenic cell death; KEGG: Kyoto Encyclopedia of Genes and Genomes; KM: Kaplan–Meier; PPI: protein–protein interaction; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas; TIDE: tumor immune dysfunction and exclusion.

Figure 2 
                  DE-ICD genes and a functional enrichment analysis of the DE-ICD genes. (a) Volcano plot of 154 up-regulated and 144 down-regulated ICD genes in the TCGA-CRC dataset (FDR < 0.05 and |logFoldchange| > 1). (b) Heatmap illustrating the DE-ICD genes in normal and CRC samples. (c) The top 30 enriched terms in GO enrichment analysis of the DE-ICD genes. (d) The top 30 enriched terms in the KEGG analysis. (e) The top 10 hub genes. (f and g) Modules obtained from the “MCODE” plugin. CRC: colorectal cancer; DE: differentially expressed; FDR: false discovery rate; ICD: immunogenic cell death; TCGA: The Cancer Genome Atlas.
Figure 2

DE-ICD genes and a functional enrichment analysis of the DE-ICD genes. (a) Volcano plot of 154 up-regulated and 144 down-regulated ICD genes in the TCGA-CRC dataset (FDR < 0.05 and |logFoldchange| > 1). (b) Heatmap illustrating the DE-ICD genes in normal and CRC samples. (c) The top 30 enriched terms in GO enrichment analysis of the DE-ICD genes. (d) The top 30 enriched terms in the KEGG analysis. (e) The top 10 hub genes. (f and g) Modules obtained from the “MCODE” plugin. CRC: colorectal cancer; DE: differentially expressed; FDR: false discovery rate; ICD: immunogenic cell death; TCGA: The Cancer Genome Atlas.

3.2 ICD-related subtypes

We determined the ICD-related clusters of TCGA-CRC through consensus clustering. Two subgroups were identified from the univariate analysis based on the prognosis-related genes (Figure 3a and b). DE-ICD gene expression levels were compared between the two clusters; cluster C1 was defined as the ICD-high group, and cluster C2 as the ICD-low group (Figure 3c). Then, a survival analysis between the two groups was performed, resulting in different clinical outcomes. The ICD-low and ICD-high subtypes were linked to favorable and unfavorable prognoses, respectively (Figure 3d).

Figure 3 
                  Identification of ICD-associated subtypes by consistent clustering. (a) Consensus clustering matrix when k = 2. (b) Relative change in the area under the cumulative distribution function curve for k = 2–9. (c) Heatmap of ICD-related gene expressions between the ICD-high and ICD-low subtypes. (d) KM curves of OS in ICD-high and ICD-low subtypes. CDF: cumulative distribution function; ICD: immunogenic cell death.
Figure 3

Identification of ICD-associated subtypes by consistent clustering. (a) Consensus clustering matrix when k = 2. (b) Relative change in the area under the cumulative distribution function curve for k = 2–9. (c) Heatmap of ICD-related gene expressions between the ICD-high and ICD-low subtypes. (d) KM curves of OS in ICD-high and ICD-low subtypes. CDF: cumulative distribution function; ICD: immunogenic cell death.

3.3 Signal pathways, somatic mutations, and tumor microenvironment differences between the ICD subtypes

Further comparison of the low and high ICD subtypes identified 1,302 dysregulated genes (Figure 4a); Figure 4b presents the top 50 genes ordered by the logFC values. GO enrichment and KEGG pathway analyses were performed using the genes up-regulated in the ICD-high subtype. These genes were enriched in biological functions and signaling pathways, including leukocyte migration, signaling receptor activator activity, cytokine activity, cytokine–cytokine receptor interaction, B-cell receptor, and the PI3K–Akt signaling pathway (Figure 4c and d). Several immune-related pathways were also associated with the ICD-high subtype, including natural killer cell-mediated cytotoxicity (NES = 2.3681, FDR < 0.0001), immune receptor activity (NES = 2.4150, FDR < 0.0001), T-cell receptor signaling pathway (NES = 2.1195, FDR = 0.0011), and B-cell receptor signaling pathway (NES = 2.2384, FDR = 0.0003), identified by GSEA (Figure 4e and f).

Figure 4 
                  Analyses of DE genes, signaling pathways, and somatic mutations between the high and low ICD subtypes. (a) Volcano plot of 1,302 dysregulated genes and (b) a heatmap presenting the top 50 genes ordered based on the log foldchange values between the two subtypes. (c) The top 10 enriched terms in GO enrichment analysis and (d) the top 30 enriched terms in KEGG analysis of the genes up-regulated in the ICD-high subtype. GSEA between the ICD- (e) high and (f) low subtypes. The most frequently mutated genes in the ICD- (g) high and (h) low subtypes. DE: differentially expressed; ICD: immunogenic cell death.
Figure 4

Analyses of DE genes, signaling pathways, and somatic mutations between the high and low ICD subtypes. (a) Volcano plot of 1,302 dysregulated genes and (b) a heatmap presenting the top 50 genes ordered based on the log foldchange values between the two subtypes. (c) The top 10 enriched terms in GO enrichment analysis and (d) the top 30 enriched terms in KEGG analysis of the genes up-regulated in the ICD-high subtype. GSEA between the ICD- (e) high and (f) low subtypes. The most frequently mutated genes in the ICD- (g) high and (h) low subtypes. DE: differentially expressed; ICD: immunogenic cell death.

Furthermore, we explored somatic mutations between the two subtypes. APC, TP53, TTN, KRAS, and SYNE1 were the top five genes with the highest mutation frequencies in both the high and low ICD subtypes, but their relative frequencies differed in each group (Figure 4g and h).

Increasing evidence has demonstrated that ICD is important in activating anti-tumor immune responses. Thus, we dissected the differences between the tumor microenvironment in the high and low ICD subtypes. The ESTIMATE algorithm showed that the ICD-high subtype had a higher immune score and lower tumor purity than the ICD-low subtype (Figure 5a). Subsequently, the CIBERSORT algorithm was used to analyze the differentiation in infiltration of immune cells between the two subtypes; Figure 5b presents the results of each TCGA-CRC sample. Patients in the ICD-high subtype had higher percentages of CD8 T cells, naive B cells, activated memory CD4 T cells, macrophages, and neutrophils than those in the ICD-low subtype (Figure 5c). In addition, most HLA and immune checkpoint genes were up-regulated in the ICD-high subtype (Figure 5d and e). Therefore, the ICD-high subtype may be related to the immune-hot phenotype.

Figure 5 
                  Tumor microenvironment analyses between the high and low ICD subtypes. (a) Violin plots of immune and tumor purity scores. (b) Relative proportion of immune infiltration in each CRC sample obtained from TCGA. (c) Immune cell infiltration differences between the high and low ICD subtypes. DE of (d) multiple HLA genes and (e) immune checkpoint genes between the high and low ICD subtypes. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001. CD: cluster of differentiation; DE: differential expression; ICD: immunogenic cell death; NK: natural killer.
Figure 5

Tumor microenvironment analyses between the high and low ICD subtypes. (a) Violin plots of immune and tumor purity scores. (b) Relative proportion of immune infiltration in each CRC sample obtained from TCGA. (c) Immune cell infiltration differences between the high and low ICD subtypes. DE of (d) multiple HLA genes and (e) immune checkpoint genes between the high and low ICD subtypes. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001. CD: cluster of differentiation; DE: differential expression; ICD: immunogenic cell death; NK: natural killer.

3.4 Establishing and validating an ICD-related prognostic risk model

To better predict the prognosis of patients with CRC, we built a prognostic model comprising six ICD genes. In the univariate Cox analysis of 298 DE-ICD genes, 41 were associated with overall survival (OS) in the TCGA-CRC cohort (Figure 6a). These genes were reserved for further analysis, and a six-gene prognostic risk model was established in the random survival forest analysis (Figure 6b and c). The risk score formula was: risk score = (−1.2621 × CD1A) + (0.7546 × TSLP) + (0.4412 × CD36) + (0.4615 × TIMP1) + (0.3879 × MC1R) + (−0.9271 × NRG1).

Figure 6 
                  Construction and validation of the ICD risk signature. (a) Volcano plot of 41 genes associated with OS identified by univariate Cox analysis. (b) Top 10 genes ranked by importance in the random forest analysis. (c) Combinations of genes as a prognostic model were ordered based on the −log10 (p-value). Distribution of individual risk scores and survival statuses and heatmaps of the prognostic six-gene signature demonstrate significant differences between the high-risk and low-risk groups in (d) TCGA-CRC and (f) GSE17538 cohorts. KM and ROC analysis in (e) TCGA-CRC and (g) GSE17538 cohorts. AUC: area under the curve; CI: confidence interval; CRC: colorectal cancer; HR: hazard ratio; ICD: immunogenic cell death; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas.
Figure 6

Construction and validation of the ICD risk signature. (a) Volcano plot of 41 genes associated with OS identified by univariate Cox analysis. (b) Top 10 genes ranked by importance in the random forest analysis. (c) Combinations of genes as a prognostic model were ordered based on the −log10 (p-value). Distribution of individual risk scores and survival statuses and heatmaps of the prognostic six-gene signature demonstrate significant differences between the high-risk and low-risk groups in (d) TCGA-CRC and (f) GSE17538 cohorts. KM and ROC analysis in (e) TCGA-CRC and (g) GSE17538 cohorts. AUC: area under the curve; CI: confidence interval; CRC: colorectal cancer; HR: hazard ratio; ICD: immunogenic cell death; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas.

Patients in the TCGA-CRC cohort were stratified into high- and low-risk groups based on the median risk score value, and the association between survival status and risk score was visualized. The number of deceased patients in the high-risk group was considerably higher than that in the low-risk group. Figure 6d illustrates the expression level differences of the six genes between the two groups.

KM analysis indicated that patients with CRC and low-risk scores in the TCGA cohort had better OS times than those with high-risk scores. The predictability and sensitivity of the risk score model were validated by ROC analysis (Figure 6e). The areas under the curve of the prognostic model at 1, 2, and 3 years in the TCGA dataset were 0.70, 0.72, and 0.74, respectively, demonstrating its excellent predictive model performance. As an independent external validation set, GSE17538 was used to further evaluate the efficiency of the constructed prognostic model, which also confirmed the model’s good performance (Figure 6f and g).

In addition, we performed univariate and multivariate Cox regression analyses, finding that age, tumor stage, and risk score were independent prognostic factors in patients with CRC (Figure 7a and b). Subsequently, a nomogram including the risk score, stage, and age was created (Figure 7c). Calibration plots revealed that the predicted probability of OS probabilities at 3, 4, and 5 years were highly concordant with the actual OS (Figure 7d–f).

Figure 7 
                  Construction and internal validation of a predictive nomogram of OS for patients with CRCs. (a) Univariate and (b) multivariate Cox analyses evaluating the independent prognostic value of the risk signature in patients with CRC. (c) Nomogram based on age, stage, and signature risk score. Calibration plots of the nomogram for predicting the probability of (d) 3-, (e) 4-, and (f) 5-year survival. CRC: colorectal cancer; OS: overall survival.
Figure 7

Construction and internal validation of a predictive nomogram of OS for patients with CRCs. (a) Univariate and (b) multivariate Cox analyses evaluating the independent prognostic value of the risk signature in patients with CRC. (c) Nomogram based on age, stage, and signature risk score. Calibration plots of the nomogram for predicting the probability of (d) 3-, (e) 4-, and (f) 5-year survival. CRC: colorectal cancer; OS: overall survival.

3.5 Correlations between tumor immune cell infiltration and immunotherapy response and the ICD risk signature

ICD plays a considerable role in anti-tumor immune response. Thus, we assessed the connections between the risk scores and infiltrating immune cells. The risk scores of samples in the TCGA-CRC cohort samples negatively correlated with activated memory CD4 T cells, CD8 T cells, and plasma cells (Figure 8a). Similar results were obtained using the GSE17538 dataset (Figure 8b).

Figure 8 
                  Correlations between the risk signature and tumor microenvironment. Scatter plots present the relationships between the risk score and the infiltration of plasma cells, activated memory CD4 T cells, and CD8 T cells in the (a) TCGA-CRC and (b) GSE17538 datasets. (c) The TIDE scores in the high- and low-risk groups. CRC: colorectal cancer; TCGA: The Cancer Genome Atlas; TIDE: tumor immune dysfunction and exclusion.
Figure 8

Correlations between the risk signature and tumor microenvironment. Scatter plots present the relationships between the risk score and the infiltration of plasma cells, activated memory CD4 T cells, and CD8 T cells in the (a) TCGA-CRC and (b) GSE17538 datasets. (c) The TIDE scores in the high- and low-risk groups. CRC: colorectal cancer; TCGA: The Cancer Genome Atlas; TIDE: tumor immune dysfunction and exclusion.

We then used the TIDE algorithm to evaluate the predictive power of the ICD-related risk score in response to immunotherapy. Compared with those in the high-risk group, patients in the low-risk group were associated with low TIDE scores in the TCGA cohort (Figure 8c), demonstrating that patients with CRC and low-risk scores should benefit from immunotherapy more than those with high-risk scores.

4 Discussion

In recent years, immunotherapy has become an important therapeutic method for patients with metastatic CRC. However, owing to the low response rate and side effects, cancer immunotherapy is effective in only a small subset of patients; however, inducing ICDs may mitigate these challenges [25]. As a new form of regulatory cell death, ICD has been found to induce adaptive immunity, thereby enhancing anti-tumor immunity. Thus, identifying ICD-related biomarkers may help distinguish CRC patients who could benefit from immunotherapy, similar to ICD markers in melanoma [26], head and neck squamous cell carcinoma [27], and neuroblastoma [28]. Furthermore, a recent study reported a breast cancer prognostic risk model comprising seven ICD genes for estimating the prognosis of patients with breast cancer and their response to immunotherapy [29], similar to our research aim.

To our knowledge, this study is the first to describe the role of ICD genes in tumor microenvironment differentiation and CRC prognosis. We identified two distinct subtypes of patients with CRC based on the expression of ICD-related genes, which correlated with diverse infiltration levels of various immune cells and different survival prognoses. The ESTIMATE algorithm showed that patients with the ICD-high subtype had high immune scores, low tumor purity, and increased immune cell infiltration. Furthermore, the up-regulated genes in the ICD-high subtype were enriched in leukocyte migration, signaling receptor activator activity, cytokine activity, and cytokine–cytokine receptor interaction; immune-related signaling pathways were also significantly enriched. Most HLA genes and immune checkpoint genes were also up-regulated in the ICD-high group. Some have reported that tumors can escape immune surveillance by reducing the expression of major histocompatibility complex (MHC) molecules [30]. Immune checkpoint inhibitors have also become a promising treatment method in cancer immunotherapy [31]. Our results demonstrated that the ICD-high subgroup was related to the immune-hot phenotype.

Moreover, we constructed and validated a six-gene (CD1A, TSLP, CD36, TIMP1, MC1R, and NRG1) prognostic risk signature. The risk score was associated with the immune context of the tumor microenvironment, and it has the potential ability to predict an immunotherapy response. TSLP is a pro-Th2 cytokine primarily expressed by epithelium cells and is thought to be a key cytokine linking innate and adaptive immune systems [32,33]. Many reports have noted that TSLP-dependent inflammatory Th2-type response and pro-tumorigenic functions play roles in pancreatic, breast, and cervical cancers [34,35,36]. TSLP also affects CRC progression by regulating the function of tumor-specific regulatory T cells [37]. CD36, a transmembrane glycoprotein receptor, is expressed in tumor, stromal, and immune cells and plays a significant role in regulating cell adhesion, immune response, metastasis, and angiogenesis in tumors [38,39]. Decreased CD36 expression might contribute to tumor cell evasion from the immune system and targeting CD36 may be a potential therapeutic strategy in cancer immunotherapy [40,41]. TIMP1 suppresses immune responses indirectly by regulating matrix metalloproteinase expression [42,43]; it also has an essential role in the development of gastrointestinal malignancies [44,45]. In patients with CRC, up-regulated TIMP1 expression is related to distant metastasis, vascular invasion, and lymphatic metastasis [46]. MC1R was initially described in cells of melanocytic origin but was later found on fibroblasts and most immune cells, indicating that it can affect innate and adaptive immunities [47,48]. Reports suggest that high MC1R expression might be significantly related to microsatellite instability and poor prognosis in CRC [49]. CD1A is an MHC class I-like molecule expressed by immune cells, including dendritic and Langerhans cells [50]. CD1A-positive dendritic cell infiltration is associated with favorable prognosis in patients with CRC [51], esophageal carcinoma [52], and breast cancer [53].

Our study has some limitations. Above all, the prognostic risk model was only validated on the GEO database, and prospective studies with more extensive sample numbers are needed to demonstrate the clinical value of the model. Furthermore, in vivo and in vitro biological experiments of the six genes in CRC are needed. Finally, this study only suggested a possible relationship between immune status and ICD subtypes or the risk model. Thus, we will collect sufficient samples in future studies to evaluate the value of this typing method and risk model in immunotherapy.

5 Conclusion

In this study, bioinformatics analyses of abnormally expressed ICD genes identified two distinct ICD subtypes in CRC: low and high. The ICD-low and ICD-high subtypes were associated with favorable and unfavorable prognoses, respectively. Moreover, we used ICD-related genes (CD1A, TSLP, CD36, TIMP1, MC1R, and NRG1) to construct and validate a prognostic risk signature of CRC; patients with low-risk scores had longer OS and might benefit from immunotherapy more than those with high-risk scores. This study highlights the predictive value of the ICD-related genes in CRC and immunotherapy response prognoses.


# Chun Yu and Weixuan Yang have contributed equally to this work.

tel: +86 13062589271, tel: +86 13621591027

  1. Funding information: The present study was supported by grants from the First Affiliated Hospital of Nanjing Medical University Clinical Capacity Enhancement Project (No. JSPH-MC-2022-5). The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

  2. Author contributions: WC and YG were involved in the study conception and the study design. All authors were involved in the implementation of the study and interpreting the data. CY, WY, LT, and YQ participated in the data acquisition, and bioinformatics analysis of the study. CY and YG, with the help of WC, wrote the manuscript. All authors have read and approved the final manuscript.

  3. Conflict of interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

  4. Data availability statement: Publicly available datasets were analyzed in this study. All data can be found in TCGA and GEO datasets (accession number: GSE17538).

Appendix

Figure A1 
                  Protein-protein interactions network of the differentially expressed (DE) - immunogenic cell death (ICD) genes.
Figure A1

Protein-protein interactions network of the differentially expressed (DE) - immunogenic cell death (ICD) genes.

References

[1] Hossain MS, Karuniawati H, Jairoun AA, Urbi Z, Ooi J, John A, et al. Colorectal cancer: a review of carcinogenesis, global epidemiology, current challenges, risk factors, preventive and treatment strategies. Cancers (Basel). 2022;14(7):1732.10.3390/cancers14071732Search in Google Scholar PubMed PubMed Central

[2] Weiser MR. AJCC 8th edition: colorectal cancer. Ann Surg Oncol. 2018;25(6):1454–5.10.1245/s10434-018-6462-1Search in Google Scholar PubMed

[3] Zhao W, Jin L, Chen P, Li D, Gao W, Dong G. Colorectal cancer immunotherapy-Recent progress and future directions. Cancer Lett. 2022;545:215816.10.1016/j.canlet.2022.215816Search in Google Scholar PubMed

[4] Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019;16(6):361–75.10.1038/s41575-019-0126-xSearch in Google Scholar PubMed PubMed Central

[5] Tirendi S, Marengo B, Domenicotti C, Bassi AM, Almonti V, Vernazza S. Colorectal cancer and therapy response: a focus on the main mechanisms involved. Front Oncol. 2023;13:1208140.10.3389/fonc.2023.1208140Search in Google Scholar PubMed PubMed Central

[6] Marcus L, Lemery SJ, Keegan P, Pazdur R. FDA approval summary: pembrolizumab for the treatment of microsatellite instability-high solid tumors. Clin Cancer Res. 2019;25(13):3753–8.10.1158/1078-0432.CCR-18-4070Search in Google Scholar PubMed

[7] Overman MJ, McDermott R, Leach JL, Lonardi S, Lenz HJ, Morse MA, et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (CheckMate 142): an open-label, multicentre, phase 2 study. Lancet Oncol. 2017;18(9):1182–91.10.1016/S1470-2045(17)30422-9Search in Google Scholar PubMed PubMed Central

[8] Overman MJ, Lonardi S, Wong KYM, Lenz HJ, Gelsomino F, Aglietta M, et al. Durable clinical benefit with nivolumab plus ipilimumab in DNA mismatch repair-deficient/microsatellite instability-high metastatic colorectal cancer. J Clin Oncol. 2018;36(8):773–9.10.1200/JCO.2017.76.9901Search in Google Scholar PubMed

[9] Galluzzi L, Vitale I, Aaronson SA, Abrams JM, Adam D, Agostinis P, et al. Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ. 2018;25(3):486–541.10.1038/s41418-018-0102-ySearch in Google Scholar PubMed PubMed Central

[10] Galluzzi L, Vitale I, Warren S, Adjemian S, Agostinis P, Martinez AB, et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J Immunother Cancer. 2020;8(1):e000337.10.1136/jitc-2019-000337corr1Search in Google Scholar PubMed PubMed Central

[11] Collett GP, Redman CW, Sargent IL, Vatish M. Endoplasmic reticulum stress stimulates the release of extracellular vesicles carrying danger-associated molecular pattern (DAMP) molecules. Oncotarget. 2018;9(6):6707–17.10.18632/oncotarget.24158Search in Google Scholar PubMed PubMed Central

[12] Krysko DV, Garg AD, Kaczmarek A, Krysko O, Agostinis P, Vandenabeele P. Immunogenic cell death and DAMPs in cancer therapy. Nat Rev Cancer. 2012;12(12):860–75.10.1038/nrc3380Search in Google Scholar PubMed

[13] Deutsch E, Chargari C, Galluzzi L, Kroemer G. Optimising efficacy and reducing toxicity of anticancer radioimmunotherapy. Lancet Oncol. 2019;20(8):e452–63.10.1016/S1470-2045(19)30171-8Search in Google Scholar PubMed

[14] Galluzzi L, Humeau J, Buqué A, Zitvogel L, Kroemer G. Immunostimulation with chemotherapy in the era of immune checkpoint inhibitors. Nat Rev Clin Oncol. 2020;17(12):725–41.10.1038/s41571-020-0413-zSearch in Google Scholar PubMed

[15] Cai J, Hu Y, Ye Z, Ye L, Gao L, Wang Y, et al. Immunogenic cell death-related risk signature predicts prognosis and characterizes the tumour microenvironment in lower-grade glioma. Front Immunol. 2022;13:1011757.10.3389/fimmu.2022.1011757Search in Google Scholar PubMed PubMed Central

[16] Hossain DMS, Javaid S, Cai M, Zhang C, Sawant A, Hinton M, et al. Dinaciclib induces immunogenic cell death and enhances anti-PD1-mediated tumor suppression. J Clin Invest. 2018;128(2):644–54.10.1172/JCI94586Search in Google Scholar PubMed PubMed Central

[17] Rodrigues MC, Morais JAV, Ganassin R, Oliveira GRT, Costa FC, Morais AAC, et al. An overview on immunogenic cell death in cancer biology and therapy. Pharmaceutics. 2022;14(8):1564.10.3390/pharmaceutics14081564Search in Google Scholar PubMed PubMed Central

[18] Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.10.1038/ncomms3612Search in Google Scholar PubMed PubMed Central

[19] Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. Gastroenterology. 2010;138(3):958–68.10.1053/j.gastro.2009.11.005Search in Google Scholar PubMed PubMed Central

[20] Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.10.1093/bioinformatics/btq170Search in Google Scholar PubMed PubMed Central

[21] Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.10.1089/omi.2011.0118Search in Google Scholar PubMed PubMed Central

[22] Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.10.1101/gr.239244.118Search in Google Scholar PubMed PubMed Central

[23] Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.10.1038/nmeth.3337Search in Google Scholar PubMed PubMed Central

[24] Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, et al. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 2020;12(1):21.10.1186/s13073-020-0721-zSearch in Google Scholar PubMed PubMed Central

[25] Li Z, Lai X, Fu S, Ren L, Cai H, Zhang H, et al. Immunogenic cell death activates the tumor immune microenvironment to boost the immunotherapy efficiency. Adv Sci (Weinheim, Baden-Wurttemberg, Ger). 2022;9(22):e2201734.10.1002/advs.202201734Search in Google Scholar PubMed PubMed Central

[26] Ren J, Yang J, Na S, Wang Y, Zhang L, Wang J, et al. Comprehensive characterisation of immunogenic cell death in melanoma revealing the association with prognosis and tumor immune microenvironment. Front Immunol. 2022;13:998653.10.3389/fimmu.2022.998653Search in Google Scholar PubMed PubMed Central

[27] Wang X, Wu S, Liu F, Ke D, Wang X, Pan D, et al. An immunogenic cell death-related classification predicts prognosis and response to immunotherapy in head and neck squamous cell carcinoma. Front Immunol. 2021;12:781466.10.3389/fimmu.2021.781466Search in Google Scholar PubMed PubMed Central

[28] Wang-Bishop L, Wehbe M, Shae D, James J, Hacker BC, Garland K, et al. Potent STING activation stimulates immunogenic cell death to enhance antitumor immunity in neuroblastoma. J Immunother Cancer. 2020;8(1):e000282.10.1136/jitc-2019-000282Search in Google Scholar PubMed PubMed Central

[29] Wang X, Huang H, Liu X, Li J, Wang L, Li L, et al. Immunogenic cell death-related classifications in breast cancer identify precise immunotherapy biomarkers and enable prognostic stratification. Front Genet. 2022;13:1052720.10.3389/fgene.2022.1052720Search in Google Scholar PubMed PubMed Central

[30] Jongsma MLM, de Waard AA, Raaben M, Zhang T, Cabukusta B, Platzer R, et al. The SPPL3-defined glycosphingolipid repertoire orchestrates HLA class I-mediated immune responses. Immunity. 2021;54(2):387.10.1016/j.immuni.2021.01.016Search in Google Scholar PubMed

[31] Huang L, Li Y, Du Y, Zhang Y, Wang X, Ding Y, et al. Mild photothermal therapy potentiates anti-PD-L1 treatment for immunologically cold tumors via an all-in-one and all-in-control strategy. Nat Commun. 2019;10(1):4871.10.1038/s41467-019-12771-9Search in Google Scholar PubMed PubMed Central

[32] Liu YJ, Soumelis V, Watanabe N, Ito T, Wang YH, Malefyt Rde W, et al. TSLP: an epithelial cell cytokine that regulates T cell differentiation by conditioning dendritic cell maturation. Annu Rev Immunol. 2007;25:193–219.10.1146/annurev.immunol.25.022106.141718Search in Google Scholar PubMed

[33] Liu YJ. TSLP in epithelial cell and dendritic cell cross talk. Adv Immunol. 2009;101:1–25.10.1016/S0065-2776(08)01001-8Search in Google Scholar PubMed PubMed Central

[34] De Monte L, Reni M, Tassi E, Clavenna D, Papa I, Recalde H, et al. Intratumor T helper type 2 cell infiltrate correlates with cancer-associated fibroblast thymic stromal lymphopoietin production and reduced survival in pancreatic cancer. J Exp Med. 2011;208(3):469–78.10.1084/jem.20101876Search in Google Scholar PubMed PubMed Central

[35] Olkhanud PB, Rochman Y, Bodogai M, Malchinkhuu E, Wejksza K, Xu M, et al. Thymic stromal lymphopoietin is a key mediator of breast cancer progression. J Immunol (Baltimore, Md: 1950). 2011;186(10):5656–62.10.4049/jimmunol.1100463Search in Google Scholar PubMed PubMed Central

[36] Xie F, Meng YH, Liu LB, Chang KK, Li H, Li MQ, et al. Cervical carcinoma cells stimulate the angiogenesis through TSLP promoting growth and activation of vascular endothelial cells. Am J Reprod Immunol (New York, NY: 1989). 2013;70(1):69–79.10.1111/aji.12104Search in Google Scholar PubMed

[37] Obata-Ninomiya K, de Jesus Carrion S, Hu A, Ziegler SF. Emerging role for thymic stromal lymphopoietin-responsive regulatory T cells in colorectal cancer progression in humans and mice. Sci Transl Med. 2022;14(645):eabl6960.10.1126/scitranslmed.abl6960Search in Google Scholar PubMed

[38] Pascual G, Avgustinova A, Mejetta S, Martín M, Castellanos A, Attolini CS, et al. Targeting metastasis-initiating cells through the fatty acid receptor CD36. Nature. 2017;541(7635):41–5.10.1038/nature20791Search in Google Scholar PubMed

[39] DeFilippis RA, Chang H, Dumont N, Rabban JT, Chen YY, Fontenay GV, et al. CD36 repression activates a multicellular stromal program shared by high mammographic density and tumor tissues. Cancer Discovery. 2012;2(9):826–39.10.1158/2159-8290.CD-12-0107Search in Google Scholar PubMed PubMed Central

[40] Wang J, Li Y. CD36 tango in cancer: signaling pathways and functions. Theranostics. 2019;9(17):4893–908.10.7150/thno.36037Search in Google Scholar PubMed PubMed Central

[41] Ruan C, Meng Y, Song H. CD36: an emerging therapeutic target for cancer and its molecular mechanisms. J Cancer Res Clin Oncol. 2022;148(7):1551–8.10.1007/s00432-022-03957-8Search in Google Scholar PubMed

[42] Ulug U, Goldman S, Ben-Shlomo I, Shalev E. Matrix metalloproteinase (MMP)-2 and MMP-9 and their inhibitor, TIMP-1, in human term decidua and fetal membranes: the effect of prostaglandin F(2alpha) and indomethacin. Mol Hum Reprod. 2001;7(12):1187–93.10.1093/molehr/7.12.1187Search in Google Scholar PubMed

[43] Schoeps B, Frädrich J, Krüger A. Cut loose TIMP-1: an emerging cytokine in inflammation. Trends Cell Biol. 2022;33(5):413–26.10.1016/j.tcb.2022.08.005Search in Google Scholar PubMed

[44] Vihinen P, Kähäri VM. Matrix metalloproteinases in cancer: prognostic markers and therapeutic targets. Int J Cancer. 2002;99(2):157–66.10.1002/ijc.10329Search in Google Scholar PubMed

[45] Brew K, Dinakarpandian D, Nagase H. Tissue inhibitors of metalloproteinases: evolution, structure and function. Biochim Biophys Acta. 2000;1477(1–2):267–83.10.1016/S0167-4838(99)00279-4Search in Google Scholar PubMed

[46] Vočka M, Langer D, Fryba V, Petrtyl J, Hanus T, Kalousova M, et al. Serum levels of TIMP-1 and MMP-7 as potential biomarkers in patients with metastatic colorectal cancer. Int J Biol Markers. 2019;34(3):292–301.10.1177/1724600819866202Search in Google Scholar PubMed

[47] Maaser C, Kannengiesser K, Specht C, Lügering A, Brzoska T, Luger TA, et al. Crucial role of the melanocortin receptor MC1R in experimental colitis. Gut. 2006;55(10):1415–22.10.1136/gut.2005.083634Search in Google Scholar PubMed PubMed Central

[48] Gruis NA, van Doorn R. Melanocortin 1 receptor function: shifting gears from determining skin and nevus phenotype to fetal growth. J Investig Dermatol. 2012;132(8):1953–5.10.1038/jid.2012.216Search in Google Scholar PubMed

[49] Peng L, Chang J, Liu X, Lu S, Ren H, Zhou X, et al. MC1R is a prognostic marker and its expression is correlated with MSI in colorectal cancer. Curr Issues Mol Biol. 2021;43(3):1529–47.10.3390/cimb43030108Search in Google Scholar PubMed PubMed Central

[50] Pereira CS, Macedo MF. CD1-restricted T cells at the crossroad of innate and adaptive immunity. J Immunol Res. 2016;2016:2876275.10.1155/2016/2876275Search in Google Scholar PubMed PubMed Central

[51] Sandel MH, Dadabayev AR, Menon AG, Morreau H, Melief CJ, Offringa R, et al. Prognostic value of tumor-infiltrating dendritic cells in colorectal cancer: role of maturation status and intratumoral localization. Clin Cancer Res. 2005;11(7):2576–82.10.1158/1078-0432.CCR-04-1448Search in Google Scholar PubMed

[52] Ikeguchi M, Ikeda M, Tatebe S, Maeta M, Kaibara N. Clinical significance of dendritic cell infiltration in esophageal squamous cell carcinoma. Oncol Rep. 1998;5(5):1185–9.10.3892/or.5.5.1185Search in Google Scholar PubMed

[53] La Rocca G, Anzalone R, Corrao S, Magno F, Rappa F, Marasà S, et al. CD1a down-regulation in primary invasive ductal breast carcinoma may predict regional lymph node invasion and patient outcome. Histopathology. 2008;52(2):203–12.10.1111/j.1365-2559.2007.02919.xSearch in Google Scholar PubMed

Received: 2023-04-20
Revised: 2023-09-22
Accepted: 2023-10-10
Published Online: 2023-11-09

© 2023 the author(s), published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 3.3.2024 from https://www.degruyter.com/document/doi/10.1515/med-2023-0836/html
Scroll to top button