Identification and validation of endoplasmic reticulum stress-related genes that enhance immunotherapy in colon cancer

Introduction

Colon cancer is a complex disease characterized by a wide range of genetic and epigenetic alterations (1). Due to constant exposure to external factors and microenvironmental pressures, the adaptive process of colon cancer cells appears, which becomes particularly complicated for immunotherapy of colon cancer (2). In recent years, the role of endoplasmic reticulum stress (ERS) and its associated genes in cancer progression and therapy response has garnered significant attention. Previous study has successfully developed and validated ERS-related gene signatures for predicting overall survival (OS) in lung adenocarcinoma (3).

The reaction process of tumor cells to immunotherapy is mainly regulated by endoplasmic reticulum, which controls the synthesis and modification of various proteins. When the homeostasis of the endoplasmic reticulum is disrupted, it is referred to as ERS (4). ERS has gained significant attention in cancer research in recent years as it is present in many types of malignancies (5). ERS pathway can be triggered by various factors, including hypoxia, oxidative stress, and nutrient deficiency, ultimately leading to the activation of the unfolded protein response (UPR). ERS-related genes interact with various signaling networks implicated in colon cancer progression, such as the Wnt/β-catenin, PI3K/Akt/mTOR, MAPK/ERK, and NF-κB (6). These interactions can influence processes like cell proliferation, migration, invasion, and angiogenesis, ultimately impacting tumor behavior and response to therapy. Nonetheless, if this process persists, it may result in a maladaptive response.

ERS and its associated UPR can significantly affect the tumor microenvironment (TME), influence survival and contribute to the drug sensitivity of cancer cells, including in colon cancer. The expression of ERS-related genes like BiP/GRP78, ATF6, and IRE1α could predict survival outcomes in patients with cancer. High expression of GRP78 is often linked to poor prognoses and is indicative of a robust adaptive UPR, which can support cancer cell survival under the harsh conditions of the TME. The ERS has been found to promote various mechanisms of tumor progression, including tumor cell survival, treatment resistance, tumor invasion, and metastasis (7).

Additionally, studies have demonstrated that ERS plays a significant role in coordinating various cellular stress signals in colon cancer (8,9). Therefore, exploring the interactions between anticancer drugs and ERS could lead to potential anti-cancer strategies that may alter disease progression (5,10). ERS-related proteins may influence the immune landscape of tumors by modulating antigen presentation and the repertoire of immune cells in the TME, potentially impacting cell survival and therapy response in colon cancer (11). Furthermore, drugs associated with ERS could offer a new perspective for the treatment of colon cancer (12).

Therefore, exploring the role of ERS in colon cancer management could be beneficial. To achieve this, we conducted bioinformatics analysis to identify ERS-related genes, investigated their molecular mechanisms and evaluated their role in the sensitivity of immunotherapy for colon cancer. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2227/rc).

Methods

Data extraction and processing

The normalized RNA-sequencing count values of colon cancer and normal tissue (n=639) were downloaded from UCSC Xena (https://xenabrowser.net/datapages/), which have integrated The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression Project (GTEx) database. The count values were used to screen out the ERS-related differentially expressed genes (DEGs).

The fragments per kilobase million (FPKM) value, clinicopathological characteristics, and follow-up information of TCGA-colon adenocarcinoma (COAD) (n=443) were downloaded from TCGA website (https://portal.gdc.cancer.gov/). And the FPKM values were transformed into transcripts per kilobase million (TPM) for further analyses. The series matrix files of GSE39582 (n=579) and GSE17536 (n=232) datasets were retrospectively downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) to externally validate results from TCGA training cohort. The data from patients without complete survival information and repeated patient records were excluded.

The ERS-related genes were acquired from Molecular Signatures Database (MSigDB; http://www.broad.mit.edu/gsea/msigdb/) and GeneCards (https://www.genecards.org/). The genes from GeneCards website with relevance score >10 were selected. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Consensus unsupervised clustering

By comparing the count values of colon normal and tumor samples using Deseq2 packages of R software (version 4.1.2), ERS-related genes with adjusted P values <0.05 and the |log₂fold change| values >1 were identified as DEGs. The volcano plot of DEGs was visualized with the EnhancedVolcano R package. To further explore inherent characteristics of ERS-related genes in colon cancer, the consensus clustering was used to classify the samples according to ERS-related DEGs. The ConsensusClusterPlus R package was used for clustering, and the optimal subgroup number was accessed using cumulative distribution function (CDF). Subsequently, the correlation between the ERS clusters and clinicopathology characteristics was analyzed by chi-square test.

ERS-related prognostic model establishment and validation

The least absolute shrinkage and selection operator (LASSO) analysis was carried out to establish a prognostic model with ERS-related DEGs based on glmnet R package. The risk score was calculated based on the gene expression and LASSO coefficients. Patients were divided into high- or low-risk group according to the median value of risk score. The Eq. [1] which was used to calculate the risk score was as follow:

Riskscore=∑i=1n[expression(gene)i×coefficient(gene)i][1]

OS analyses of TCGA cohort were carried out using survival R packages, which was further validated by the GSE39582 and GSE17536 testing sets. The prognostic accuracy of model was evaluated using time-dependent receiver operating characteristic (ROC) and concordance index (C-index) with timeROC and pec R package.

To explore whether ERS-related group is an independent predictor of colon cancer, the risk score and other clinical features were brought into in the univariate and multivariate Cox analysis. Besides, a nomogram that integrated the risk group and clinical characteristic was developed to predict OS rates suing rms R package. Subsequently, calibration plots were applied to explore the predictive performance of the nomogram.

Biological function and mechanisms analysis

The Gene Ontology (GO) and the gene set enrichment analysis (GSEA) were used to investigate the potential biological function and signaling pathway between the two risk groups. The R packages of enrichplot, clusterProfiler, limma, and org.Hs.eg.db were used for the above analysis.

Single-sample GSEA (ssGSEA), ESTIMATE, immune checkpoint analysis, and drugs response

The infiltrating scores of 24 immune cell subsets of samples were calculated by the ssGSEA algorithm based on gene set variation analysis (GSVA) R package (13). In addition, ESTIMATE algorithm was used to evaluate the immune and stromal scores between two groups with ESTIMATE R package (14).

The correlation of risk score and immune checkpoints was employed to examine the immunotherapeutic responses. Ultimately, pRRophetic R package was used to evaluate chemotherapeutic drugs response by the semi-inhibitory concentration (IC₅₀) of colon cancer patients.

Statistical analysis

All statistical analyses were conducted using R (version 4.3). We considered P values less than 0.05 as statistically significant, using a two-sided approach. Descriptive statistics were employed to analyze data from colon cancer patients in TCGA, GTEx, and GEO. Categorical variables were presented using frequencies and proportions, while the correlation between the ERS clusters and clinicopathology characteristics was analyzed by chi-squared test.

Results

Consistent clustering and clinicopathology characteristics

There were 178 genes identified as ERS-related DEGs, a volcano plot revealed the top 20 of ERS-related DEGs (Figure 1A). The consistent clustering of the TCGA cohort was based on the DEGs obtained previously. According to the CDF curve, k=3 appeared to be selected optimally and TCGA-COAD patients were divided into three subtypes (Figure 1B-1D). Among the clinicopathology features, histologic type (P<0.001), tumor stage (P=0.03), metastasis (P=0.006), and colon cancer subtype (P<0.001) were associated with ERS cluster acquired previously (Figure 1E).

Figure 1 Consistent clustering based on the ERS-related DEGs and clinicopathology characteristic. (A) The volcano plot of the top 20 ERS-related DEGs. (B,C) The CDF curve, tracking plot for k=2 to 9. (D) Consensus clustering matrix of ERS patterns in TCGA-COAD. (E) The correlation between ERS cluster and clinicopathology characteristic. *, P=0.03; **, P=0.006; ***, P<0.001. CDF, cumulative distribution function; CIN, chromosomal instability; GS, genome stable; HM-SNV, hypermutated-single nucleotide variant; MSI, microsatellite instability; ERS, endoplasmic reticulum stress; DEGs, differentially expressed genes; TCGA, The Cancer Genome Atlas; COAD, colon adenocarcinoma.

Risk model establishment and validation

Using the LASSO regression, we extracted six ERS-related genes when the value of log (λ) was the minimum likelihood of deviance (Figure 2A,2B). The LASSO regression analysis was used to establish a prognostic model with ERS-related genes in COAD.

Theprognosticriskscore=PMM2×(−0.3472)+STC2×0.0603+EIF2AK1×0.0520+HSPA1A×0.0480+SLC8A1×(−0.0403)+KCNQ1×(−0.0682)[2]

Figure 2 Establishment of risk model based on the LASSO model. (A,B) Partial likelihood deviance versus log (λ) in the LASSO model and coefficients of six ERS-related DEGs. (C-E) The distribution of risk score, survival time, clinical endpoint and the heatmap of six ERS-related genes expression in TCGA-COAD dataset. (F) K-M survival plots of the two risk groups in the TCGA-COAD dataset. (G,H) K-M survival plots of the two risk groups in the GSE39582 and GSE17538 cohort. TCGA, The Cancer Genome Atlas; LASSO, least absolute shrinkage and selection operator; ERS, endoplasmic reticulum stress; DEGs, differentially expressed genes; COAD, colon adenocarcinoma; K-M, Kaplan-Meier.

To identify the ERS-related risk model in patient survival prediction, we divided the TCGA training cohort samples into two risk groups according to median risk score (Figure 2C-2E). Kaplan-Meier (K-M) survival plots of TCGA training cohort demonstrated that patients with higher risk score had significant poorer prognosis than those with lower risk score (Figure 2F). We verified it in the GSE39582 and GSE17538 cohorts from GEO database (Figure 2G,2H).

The area under the curve (AUC) values of the risk model for 1-, 3-, and 5-year were 0.721, 0.731, and 0.770 respectively, suggesting that the risk model had accurate predictive performance (Figure 3A). When compared with other clinicopathological factors, the risk score showed better predicative performance and AUC of 3-year was 0.731 (Figure 3B). In addition, the 3- and 5-year C-index of risk score was 0.701 and 0.713 respectively, which was better than other clinical characteristics (Figure 3C). Then, we further employed the univariate and multivariate Cox analyses, demonstrating that the prognostic value of ERS-related risk group (Figure 3D,3E).

Figure 3 The risk model validation, nomogram establishment and calibration curves. (A) The ROC curves of TCGA cohort. (B) The ROC curves of risk score combined with clinical characteristics in the TCGA cohort. (C) The C-index curves of risk model. (D,E) Univariate and multivariate Cox analyses of risk group and other clinical characteristic. (F) Constructed nomogram for predicting OS. (G) The calibration plot for 1-, 3-, and 5-year OS. AUC, area under the curve; CI, confidence interval; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Nomogram construction

A nomogram was established for predicting the OS probability based on risk group and other clinical features (Figure 3F). The calibration plot of the nomogram demonstrated good consistency between the actual observation and prediction (Figure 3G).

GO and GSEA

According to molecular function of GO, the DEGs between the two groups were enriched in monocarboxylic acid binding, carboxylic acid binding, phosphatidylcholine binding, and so on (Figure 4A). The results of GSEA demonstrated that the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of two risk groups were mainly concentrated in biosynthesis and metabolism, such as butanoate metabolism, fatty acid metabolism and propanoate metabolism (Figure 4B).

Figure 4 GO and GSEA. (A) GO pathway enrichment between the two groups. (B) The GSEA analysis: the KEGG pathways of two risk groups. BP, biological process; CC, cell composition; MF, molecular function; BMP, bone morphogenic protein; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; GSEA, gene set enrichment analysis.

TME, immune checkpoint analysis, and drugs sensitivity

Based on ssGSEA algorithms, the significant differences of the immune cell in TME were observed between the two risk groups (Figure 5A). The infiltration levels of cytotoxic cells, dendritic cell (DC), immature DC (iDC), macrophages, natural killer (NK) cells, plasmacytoid DC (pDC), T helper (Th)1 cells were evidently higher in the high-risk group, while B cells, Th cells, and Th17 cells had lower infiltration in the high-risk group. Moreover, the patients in high-risk group showed higher stromal scores (P=0.007) and ESTIMATE score (P=0.03) than those in low-risk group (Figure 5B).

Figure 5 TME and drugs sensitivity. (A) The ssGSEA for the association between immune cell subpopulations of two risk groups. *, P<0.05; **, P<0.01; ***, P<0.001. (B) TME score in the two risk groups. *, P=0.03; **, P=0.007. (C) Correlations between risk score and immune checkpoint. *, P<0.05. (D) The IC₅₀ values evaluate drugs sensitivity between two groups. aDC, activated dendritic cell; DC, dendritic cell; iDC, immature dendritic cell; NK, natural killer; pDC, plasmacytoid dendritic cell; Tcm, central memory T cell; Tem, effector memory T cell; TFH, follicular helper T cell; Tgd, gamma delta T cell; Th, T helper; TReg, regulatory T cell; TME, tumor microenvironment; IC₅₀, semi-inhibitory concentration; ssGSEA, single-sample gene set enrichment analysis.

The correlational analysis showed that the risk score was positively correlated with the expression of PD1 (PDCD1) (Figure 5C). By drug sensitivity comparison, the colon cancer patients in high-risk groups showed lower IC₅₀ in A.770041, ABT.888, AG.014699, and so on (Figure 5D), indicating the patients were more sensitive to these drugs.

Discussion

Colon cancer is a prevalent and deadly form of cancer. The regulation and control of ERS signaling is primarily managed by three transducers: ATF6, IRE1, and PERK (9). ERS has been linked to a range of diseases, including cardiovascular and rheumatic diseases (5,15).

The ERS has gained significant attention in cancer research due to its connection to cellular functions and its role in maintaining and restoring metabolism (16). Research in recent years has revealed that the ERS plays crucial roles in colon cancer initiation, progression, and apoptosis (17).

Activation of ERS plays a role in regulating the differentiation of epithelial stem cells in the mouse intestine (18). Additionally, ERS has been shown to induce differentiation in colon cancer stem cells and is relevant to chemotherapy sensitivity. Furthermore, increased ERS and nuclear reprogramming can lead to a pro-metastatic state in cancer (19). As such, exploring ERS further could offer a potential therapeutic approach for preventing and treating colon cancer (20,21).

In the study, we first employed consensus clustering to identify clusters with ERS-related DEGs. The differences in clinical features between these clusters confirmed the heterogeneity of ERS in colon cancer. We established a prognostic model using the LASSO method and calculated a risk score for each patient based on candidate genes to identify patterns. We determined the prediction performance of ERS risk model based on TCGA training set and verified it in the external GSE39582 and GSE17538 queues from GEO database. We also established a nomogram combining the risk model and the clinical characteristics of colon cancer, and visualized the prognostic factors by quantitative methods. The greatest advantage of the risk model developed in this study is its potential to enhance the sensitivity of immunotherapy for colon cancer. By identifying specific ERS-related genes (PMM2, STC2, EIF2AK1, HSPA1A, SLC8A1, KCNQ1) that can improve the effectiveness of immunotherapy, this model offers a personalized approach to treatment. This tailored strategy may lead to more targeted and successful interventions, improving patient outcomes and potentially reducing the need for broader, less specific treatments.

In recent years, the immune landscape has been attached great importance in clinical cancer research. Increasing evidence highlights the important role of TME in progression of colon cancer (22,23). We proved that these are closely related to tumor immunity and drug sensitivity through GO analysis and GSEA. We investigated the association between ERS risk model and tumor immune cell infiltration. The TME of the two ERS risk groups were explored by ssGSEA and ESTIMATE algorithm. These findings suggested that ERS plays an important role in immunotherapy of colon cancer.

Cumulative evidence demonstrates that ERS is relevant with chemotherapy resistance in cancers (23). We found that colon cancer patients in high-risk groups may be more sensitive to specific drugs according to drug sensitivity comparison. The roles of ERS in colon cancer may be helpful to explore potential therapeutic strategies. To achieve more precise and personalized immunotherapy responses based on the results of this study, we will establish these ERS-related genes as biomarkers for predicting immunotherapy responses in colon cancer patients.

However, there are several challenges in translating these results into clinical applications. One significant hurdle is the need for further validation and verification through in vivo and in vitro experiments to confirm the efficacy and safety of targeting these ERS-related genes in clinical settings. Additionally, the complexity of biological systems also pose challenges in predicting and ensuring the success of such targeted interventions (24).

Conclusions

We established a risk model with ERS-related genes (PMM2, STC2, EIF2AK1, HSPA1A, SLC8A1, KCNQ1), which enhance the sensitivity of immunotherapy for colon cancer. These findings may provide a new perspective for the individualized treatment of colon cancer.

Acknowledgments

Funding: This work was financially supported by the Military Logistics Research Program of PLA (No. ZZCWS23J1025).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2227/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2227/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2227/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Okugawa Y, Grady WM, Goel A. Epigenetic Alterations in Colorectal Cancer: Emerging Biomarkers. Gastroenterology 2015;149:1204-1225.e12. [Crossref] [PubMed]
2. Obacz J, Avril T, Rubio-Patiño C, et al. Regulation of tumor-stroma interactions by the unfolded protein response. FEBS J 2019;286:279-96. [Crossref] [PubMed]
3. Lin L, Zhang W. Development and validation of endoplasmic reticulum stress-related eight-gene signature for predicting the overall survival of lung adenocarcinoma. Transl Cancer Res 2022;11:1909-24. [Crossref] [PubMed]
4. Pastor-Cantizano N, Ko DK, Angelos E, et al. Functional Diversification of ER Stress Responses in Arabidopsis. Trends Biochem Sci 2020;45:123-36. [Crossref] [PubMed]
5. Liu CY, Hsu CC, Huang TT, et al. ER stress-related ATF6 upregulates CIP2A and contributes to poor prognosis of colon cancer. Mol Oncol 2018;12:1706-17. [Crossref] [PubMed]
6. Li W, Jin K, Luo J, et al. NF-κB and its crosstalk with endoplasmic reticulum stress in atherosclerosis. Front Cardiovasc Med 2022;9:988266. [Crossref] [PubMed]
7. Madden E, Logue SE, Healy SJ, et al. The role of the unfolded protein response in cancer progression: From oncogenesis to chemoresistance. Biol Cell 2019;111:1-17. [Crossref] [PubMed]
8. Guo B, Xiong X, Hasani S, et al. Downregulation of PHLPP induced by endoplasmic reticulum stress promotes eIF2α phosphorylation and chemoresistance in colon cancer. Cell Death Dis 2021;12:960. [Crossref] [PubMed]
9. Piton N, Wason J, Colasse É, et al. Endoplasmic reticulum stress, unfolded protein response and development of colon adenocarcinoma. Virchows Arch 2016;469:145-54. [Crossref] [PubMed]
10. Marciniak SJ, Chambers JE, Ron D. Pharmacological targeting of endoplasmic reticulum stress in disease. Nat Rev Drug Discov 2022;21:115-40. [Crossref] [PubMed]
11. Katoh Y, Yaguchi T, Kubo A, et al. Inhibition of stearoyl-CoA desaturase 1 (SCD1) enhances the antitumor T cell response through regulating β-catenin signaling in cancer cells and ER stress in T cells and synergizes with anti-PD-1 antibody. J Immunother Cancer 2022;10:e004616. [Crossref] [PubMed]
12. Fan F, Shen P, Ma Y, et al. Bullatacin triggers immunogenic cell death of colon cancer cells by activating endoplasmic reticulum chaperones. J Inflamm (Lond) 2021;18:23. [Crossref] [PubMed]
13. Finotello F, Trajanoski Z. Quantifying tumor-infiltrating immune cells from transcriptomics data. Cancer Immunol Immunother 2018;67:1031-40. [Crossref] [PubMed]
14. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
15. Ren J, Bi Y, Sowers JR, et al. Endoplasmic reticulum stress and unfolded protein response in cardiovascular diseases. Nat Rev Cardiol 2021;18:499-521. [Crossref] [PubMed]
16. Lemmer IL, Willemsen N, Hilal N, et al. A guide to understanding endoplasmic reticulum stress in metabolic disorders. Mol Metab 2021;47:101169. [Crossref] [PubMed]
17. Piao MJ, Han X, Kang KA, et al. The Endoplasmic Reticulum Stress Response Mediates Shikonin-Induced Apoptosis of 5-Fluorouracil-Resistant Colorectal Cancer Cells. Biomol Ther (Seoul) 2022;30:265-73. [Crossref] [PubMed]
18. Dilly AK, Honick BD, Lee YJ, et al. Synergistic apoptosis following endoplasmic reticulum stress aggravation in mucinous colon cancer. Orphanet J Rare Dis 2020;15:211. [Crossref] [PubMed]
19. Conod A, Silvano M, Ruiz I, Altaba A. On the origin of metastases: Induction of pro-metastatic states after impending cell death via ER stress, reprogramming, and a cytokine storm. Cell Rep 2022;38:110490. [Crossref] [PubMed]
20. Zhang X, Li W, Ma Y, et al. High-fat diet aggravates colitis-associated carcinogenesis by evading ferroptosis in the ER stress-mediated pathway. Free Radic Biol Med 2021;177:156-66. [Crossref] [PubMed]
21. Laudisi F, Di Grazia A, De Simone V, et al. Induction of endoplasmic reticulum stress and inhibition of colon carcinogenesis by the anti-helmintic drug rafoxanide. Cancer Lett 2019;462:1-11. [Crossref] [PubMed]
22. Angell HK, Bruni D, Barrett JC, et al. The Immunoscore: Colon Cancer and Beyond. Clin Cancer Res 2020;26:332-9. [Crossref] [PubMed]
23. Jahanafrooz Z, Mosafer J, Akbari M, et al. Colon cancer therapy by focusing on colon cancer stem cells and their tumor microenvironment. J Cell Physiol 2020;235:4153-66. [Crossref] [PubMed]
24. An Y, Wang X, Guan X, et al. Endoplasmic reticulum stress-mediated cell death in cardiovascular disease. Cell Stress Chaperones 2024;29:158-74. [Crossref] [PubMed]