A Prognostic Nomogram Based on an Immunogenic Cell Death and Endoplasmic Reticulum Stress-Related Gene Signature for Ovarian Cancer

Introduction

Ovarian cancer (OC) is one of the leading cause of cancer-related death among women worldwide,¹ with an estimated 140,000 deaths per year.² The standard treatment regimens for OC are highly relied on surgery and chemoradiotherapy. However, many patients develop chemotherapy resistance and relapse, yielding to a 5-year survival rate of less than 35%.³ The traditional FIGO staging system has been adopted as the standard assessment option for estimating the treatment and outcome of OC, but the estimation is hampered by the heterogeneity of OC that is constantly changing. Thus, it is urgent to develop a convenient and effective programmatic screening method for high-risk OC populations.

Molecular biomarkers have been well acknowledged to contribute to predict the prognosis of OC. Compared with traditional single-gene biomarkers, gene signatures comprising several genes provide strong evidence for patient prognosis and survival. In recent years, several prognostic models for OC have been proposed based on clinical characteristics or gene expression profiles;^4–6 however, few of them integrate immunogenic cell death (ICD) and endoplasmic reticulum (ER) stress, which are increasingly recognized as key processes in tumor immunity and treatment response.

As a type of regulatory cell death usually elicited by radiotherapy and some chemotherapeutic agents,⁷ ICD is capable of triggering a fully antigen-specific adaptive immune response by releasing specific molecules belonging to damage-associated molecular patterns (DAMPs) such as calreticulin, ATP, and high mobility group protein Box 1.^8,9 This immunogenic process is partly driven by the pre-apoptotic translocation of calreticulin from the ER to the cell surface, which represents a hallmark mechanism of ICD.¹⁰ Moreover, ICD based on ER stress has been widely studied as the fundamental mechanism of cancer immunotherapy.¹¹ These lead us to identify a gene signature associated with ICD and ER stress in OC patients. In the present study, we aim to screen gene signatures associated with ICD and ER stress and to develop a nomogram for predicting the prognosis of OC based on a risk score model. Such a nomogram may provide clinicians with a practical tool for individualized risk stratification and potentially guide personalized treatment strategies for OC patients.

Materials and Methods

Data Collection

We first screen the OS status of OC patients from the gene expression omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). In total, 6 datasets (ie GSE14764, GSE19829, GSE30161, GSE26193, GSE31245 and GSE63885) were initially included in this study (Table 1). The inclusion criteria for datasets were: (1) availability of OS information, and (2) inclusion of OC patients. Datasets that did not meet these criteria or had missing OS data were excluded. In total, 448 OC patients with OS data were used as the training set, which consisted of 80 patients from GSE14764, 70 from GSE19829, 58 from GSE30161, 107 from GSE26193, 58 from GSE31245, and 75 from GSE63885. Batch effects were removed by sva package (R 4.2.1 version) according to the previous description.¹² The effectiveness of batch effect removal was visually validated using principal component analysis (PCA) plots to confirm that samples from different datasets were more closely clustered after batch effect removal. The normalization process involved log2 (FPKM+1) transformation, which was applied to the datasets for downstream analysis. The normalization process was further validated by comparing distributions across datasets both visually (via boxplots) and statistically (via Kolmogorov–Smirnov test). Besides, 357 OC patients with OS information downloaded from the Cancer Genome Atlas (TCGA) database (https://gdc-portal.nci.nih.gov/) served as the validation set. Then log2 (FPKM+1) normalization was utilized for the downstream analysis of the selected datasets.

Table 1 Information for the Six Datasets

Identification of ICD Subtypes

Thirty-four human ICD genes were retrieved from a previous study.¹³ Expression data for 31 of these genes were available. The ICD subtype analysis was performed to assess the correlation between ICD expression and OC survival using the ConsensusClusterPlus package (R 4.21.1 version), as described previously.¹⁴ ICD subtypes were classified into two clusters (subtype A and B) based on an optimal subtype K value, defined within a range of 2–6. Chi-square test was performed to assess the differences in ICD subtypes across clinicopathological features, including age, chemotherapy response, histological type, grade, and staging. Kaplan-Meier curve was performed to compare outcomes between patients with different ICD subtypes.

Identification of ICD-Based ER Stress Genes

Differentially expressed ICD genes were visualized using the Volcano plot and heatmaps generated with the Limma package.¹⁵ Multiple testing correction was performed using Benjamini & Hochberg method to obtain the adjusted P value. Differential expression was considered significant for genes with an adjusted P value of < 0.05.

ER stress-related genes were retrieved from GeneCards Suite database (https://www.genecards.org/) using the keywords of “Endoplasmic reticulum stress”. Candidate genes were selected with a relevance score of > 7. Candidate ICD-based ER stress genes were finally identified by intersecting the differential expressed ICD genes with the ER stress genes.

Functional and Pathway Analyses

The Molecular Signatures Database gene sets C5.go.mf. V7.4.symbols.gmt was utilized to evaluate the relevant pathways and molecular mechanisms of ICD genes, according to the previous description.¹⁶ Gene set enrichment analysis (GSEA) was utilized to screen the significantly related pathways.¹⁷ DAVID package¹⁸ was performed for the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses.

Calculation of Prognostic Risk Score Based on ICD-Based ER Stress Genes

Univariate Cox regression analysis was used to identify the ICD-based ER stress genes associated with the prognosis of OC using the survminer package (https://cran.rstudio.com/web/packages/survminer/index.html). Subsequently, the least absolute shrinkage and selection operator (LASSO) regression analysis was performed. The optimal lambda parameter for LASSO was chosen using 10-fold cross-validation with the minimum criteria for the mean squared error. Afterwards, Stepwise Cox regression was further utilized to construct the prognostic risk score model based on the following formula: risk score =β1X1 + β2X2 +… + βnXn, where β represented the regression coefficient and X represented the gene expression, respectively. Then the OC patients were classified into high- and low-risk groups based on the median risk score, followed by Kaplan-Meier curve to compare the survival between the two risk groups.

Comparison of Immune Microenvironment Between Low- and High-Risk Groups

We then evaluate the immune microenvironment between high-risk and low-risk groups with three algorithms. EstimateR was used to calculate the stromal and immune scores in the low-risk and high-risk groups. MCPcounter and GAVA package were used to estimate the infiltration of immune cells. Subsequently, the gene sets obtained from previous study¹⁹ were used to mark each immune cell type. As HLA plays an important role in the immune system, the tumor antigens must be presented in an HLA-dependent manner.^20,21 Therefore, we compared the expression of HLA genes between low-risk and high-risk groups.

Establishing a Nomogram for Predicting the OS of OC Patients

Univariate and multivariate Cox regression analyses were performed to select the risk factors for predicting the prognosis of OC patients. Then a nomogram was established based on the selected risk factors using the rms package (https://cran.r-project.org/web/packages/rms/). The prediction efficacy of the nomogram was evaluated using the C-index and calibration curve.

Statistical Analysis

The statistical analysis was performed using R software (4.4.2 version). Univariate logistic regression analysis was performed to compare the variables between low- and high-risk groups. Then variables with a P value of < 0.05 were adopted into the multivariate regression analysis. A P value of less than 0.05 was considered to be statistical significance.

Results

ICD Subtypes Were Effective to Distinguish the Survival of OC Patients

Consensus cumulative distribution function (CDF) was shown in Figure 1a with a K value in a range of 2–10. The obtained 31 ICD genes were classified into two subtypes (ie subtype A and subtype B) (Figure 1b), and their expression in different subtypes was shown in Figure 1c. The survival of patients with ICD genes classified into subtype A showed significantly higher OS compared with those with ICD genes classified into B subtype (Figure 1d, P = 0.01). Supplementary Figure 1 summarized the clinical data of the patients with different subtypes, which indicated no statistical differences in the clinical features between the two subtypes, including age (P=0.60), grade (P=0.315), histological type (P=0.632), stage (P=1) and response to chemotherapy (P=0.965).

Figure 1 Immunogenic cell death (ICD) subtypes (ie up-regulation and down-regulation) analysis. (a) Consensus clustering cumulative distribution function analysis. (b) Optimal consensus clustering matrix with a K value of 2. (c) Heat plot of the expression of 31 ICD genes. Red and blue color represented up-regulation and down-regulation, respectively. (d) Survival curve of different ICD subtypes.

ICD Subtypes Were Enriched to Immunity and ER Stress

GSEA analysis for the ICD subtypes showed that 74 significant pathways were significantly associated with immunity such as B cell receptor signaling pathway and T cell receptor signaling pathway (Figure 2a and b). GO enrichment analysis showed 2706 significant enriched function, among which three were associated with ER, including calcium ion homeostasis, release of sequestered calcium ion into cytosol, and intrinsic apoptotic signaling pathway in response to ER stress (Figure 2c–e).

Figure 2 Functional pathways between ICD subtypes. (a and b) GSEA indicated the genes were enriched in B cell receptor signaling pathway and T cell receptor signaling pathway. (c–e) GO analysis showed that the genes were enriched in ER calcium ion homeostasis, release of sequestered calcium ion into cytosol by ER, and intrinsic apoptotic signaling pathway in response to ER stress.

Function Analysis of ICD-Based ER Stress Genes

In total, 3485 DEGs were selected for the functional analysis of ER stress genes based on ICD genes (Figure 3a). The expression of ICD genes was presented in Figure 3b, which showed their expression was different between different subtypes. After intersecting ICD genes with ER stress genes, a total of 245 differentially expressed ICD-based ER stress genes were retrieved (Figure 3c). Functional analysis suggested that these were enriched into 148 KEGG pathways and 4002 GO terms. The most significant pathways and GO terms were visualized in Figure 3d and e.

Figure 3 Identification and analysis for ER stress genes based on ICD genes. (a) Volcano plot for identification of differential expressed ICD genes. (b) Differential expressed ICD genes in different ICD subtypes. (c) Veen plot for intersecting the differential expressed ICD genes and ER genes. (d) KEGG pathways analysis. (e) GO function analysis. Abbreviations: BP, biological process; CC, cellular component; MF, molecular function.

Prognostic Analysis Based on Risk Scores of ICD-Based ER Stress Genes

Among the 245 ICD-based ER stress genes, 25 ICD-based ER stress genes were identified to be related to the OC prognosis after univariate Cox analysis (Figure 4a). Then 18 optimal genes were identified after LASSO analysis (Figure 4b). Finally, based on Step-wise Cox regression, 7 genes (ie PTPN2, RPS27A, NR3C1, SERPINA1, CYBA, APOE, and UBXN8) were selected (Figure 4c). OC patients were divided into high- and low-risk subgroups based on the calculated the risk score of each patient. Kaplan-Meier analysis showed that the survival in the patients with high-risk subgroup was significantly lower than the counterparts of a low-risk in the training set (P < 0.01, Figure 4d) and validation set (P<0.01, Figure 4e).

Figure 4 Construction and validation of the risk score. (a) Univariate Cox regression analysis for identification of ER stress genes related to the prognosis of OC. (b) LASSO analysis. Left; distribution of LASSO coefficient, right; the likelihood deviation of the LASSO coefficient distribution, the two vertical dashed lines represented lambda.min and Lambda.1SE, respectively. (c) Multivariate Cox analysis results for the ER stress genes related to OC prognosis. (d and e) Survival curve of different risk groups in GEO dataset and TCGA dataset.

Different Risk Groups Showed Different Immune Features

In this section, we first analyzed the in infiltration between different high- and low-risk groups, which indicated that high-risk group showed higher immune and stromal scores compared with the low-risk group (Figure 5a and b). Then MCPcounter algorithm was used to compare the infiltration of the 10 selected immune cells between low- and high-risk groups including CD8 T cells, NK cells, neutrophils, T cells, B lineage, myeloid dendritic cells, cytotoxic lymphocytes, monocytic lineage, endothelial cells and fibroblasts (Figure 5c). These immune cells were significantly correlated with the genes in risk model (Figure 5d). Subsequently, ssGSEA algorithm was used to compare the enrichment score of 28 immune cells between low- and high-risk groups (Figure 5e). These immune cells were also significantly correlated with the genes in risk model (Figure 5f). Finally, we also found that the HAL-E and HAL-J were significantly correlated with the genes in risk mode (Figure 5g).

Figure 5 Immune status analysis in different risk groups. (a and b) Immune score and stromal score in different risk groups, L, low-risk group; H, high-risk group. (c) Heat plot of immune cells in MCPCounter algorithm. (d) Correlation of immune cells and genes in risk score model. The most positive and negative correlated plot were shown in the left. (e) Heat plot of immune cells in ssGSEA algorithm. (f) Correlation of immune cells and genes in risk score model. The most positive and negative correlated plot were shown in the left. (g) Correlation of human leukocyte antigen (HLA) gens and genes in risk score model. The most positive and negative correlated plot were shown in the left.

Establishing a Nomogram for Predicting the 1-, 2- and 3-Year Survival of OC Patients

Univariate and multivariate Cox regression analysis were used to screen the independent risk factors associated with the prognosis of OC. Univariate analysis results were shown in Figure 6a. Multivariate analysis further showed that risk score and grade were independent risk factors for reduced OS in OC patients, while chemotherapy response was associated with extended OS (Figure 6b). Then these factors were included into the nomogram model (Figure 6c). The prediction efficacy was acceptable with a AUC of 0.61, 0.64 and 0.67 for the 1-, 2- and 3-year OS in the training set, while that in the test set was 0.60, 0.61 and 0.60, respectively (Figure 6d). It showed a good consistency between the predictive results and the actual results as shown in calibration curve (Figure 6e).

Figure 6 Establishing a nomogram for predicting the 1-, 2- and 3-year survival of the OC patients. (a) Univariate Cox regression analysis was performed to select the factors associated with OC prognosis. (b) Multivariate Cox regression analysis for identification of independent risk factors for OC prognosis. (c) Nomogram model. (d) Evaluation of the nomogram based on AUC in the training (left) and test (right) sets. (e) Calibration curve.

Discussion

ICD driven by ER stress has been widely accepted as the fundamental mechanism of cancer immunotherapy.⁷ Cancer cells can activate specific signaling pathways to induce ER stress and expose DAMPs, thereby promoting immune responses in the anti-cancer process.²² However, the roles of ICD-related ER stress genes in OC prognosis remain poorly defined. In the present study, we established a prognostic risk score based on seven ICD- and ER stress-associated genes, including PTPN2, RPS27A, NR3C1, SERPINA1, CYBA, APOE, and UBXN8. These genes are involved in protein homeostasis and ER-associated degradation (ERAD) pathway (eg APOE, RPS27A, and UBXN8), inflammatory and immune regulation (eg PTPN2, NR3C1, and SERPINA1) and reactive oxygen species (ROS) generation and oxidative stress (CYBA). Finally, we developed a nomogram based on this seven-gene signature, which effectively predicts the 1-, 2-, and 3-year survival of OC patients.

UBXN8, RPS27A, and APOE play key roles in regulating proteostasis and modulating ER stress responses. Among these, UBXN8 acts as a critical regulator of the ERAD pathway, facilitating the retrotranslocation and degradation of misfolded proteins through its interaction with the p97/VCP complex.²³ Dysfunction of UBXN8 triggers sustained ER stress and pro-apoptotic signaling, which in turn promotes ICD. Similarly, RPS27A, which encodes a ubiquitin-ribosomal fusion protein, contributes to protein quality control by promoting proteasomal degradation.²⁴ Its deficiency may impair the clearance of unfolded proteins, activating the UPR. Although APOE is primarily known for its role in lipid transport, it can influence ER stress indirectly by modulating lipid composition and membrane integrity, both essential for ER function.²⁵ Dysregulation of APOE has been linked to increased ER stress and protein aggregation. Collectively, these genes underscore the pivotal role of protein homeostasis in connecting ER stress to cell fate decisions such as apoptosis and ICD.

PTPN2, NR3C1, and SERPINA1 genes are key regulators of inflammation and immune homeostasis, with direct implications for ER stress and ICD. PTPN2 encodes a tyrosine phosphatase that negatively regulates cytokine signaling through pathways such as JAK/STAT and NF-κB.²⁶ Loss of PTPN2 function enhances inflammatory responses and increased susceptibility to ICD.²⁷ Similarly, activation of NR3C1 reduces DAMP release and immune cell recruitment, potentially diminishing the efficacy of ICD.^28,29 SERPINA has dual immunomodulatory functions. Under normal conditions, it suppresses the inflammation, whereas the misfolded variants induce ER stress and inflammatory responses.³⁰ Taken together, these genes could modulate both immunosuppressive and immunostimulatory pathways, influencing ER stress outcomes and thereby modulating the ICD process.

Encoding a subunit of the NADPH oxidase complex, CYBA is a key contributor to ROS in immune and epithelial cells.³¹ Elevated ROS can directly induce ER stress by causing oxidative damage to proteins and disrupting the redox environment. Additionally, CYBA plays a crucial role in enhancing the immunogenicity of dying cells by promoting the release of ICD hallmarks such as calreticulin and ATP. CYBA-mediated oxidative stress has been shown to amplify these effects, linking ROS generation with pro-inflammatory and immunostimulatory outcomes. On this basis, modulating CYBA expression or activity to regulate immunogenic apoptosis and redox balance may offer therapeutic potential in ovarian cancer.

Gene signatures have been recently utilized in developing nomograms for predicting the prognosis of OC. For instance, Wang et al developed a nomogram using a seven-gene signature based on tumor stem cell-related genes to predict the OC prognosis, which showed that the nomogram combining the gene signature demonstrated strong predictive ability.³² Similarly, Pan and Ma proposed a six-gene epithelial–mesenchymal transition–related signature for prognostic prediction in OC.³³ However, few studies have integrated ICD and ER stress signatures, which are increasingly recognized as critical processes in tumor immunity and therapy response. In our study, an optimal seven-gene signature associated with ICD and ER stress was identified, and the generated risk scores effectively predicted the OC prognosis. Notably, both the signature and the risk score model were significantly correlated with the immune statues, suggesting potential utility in predicting immune response to immunotherapy. We further constructed a nomogram incorporating the risk score, tumor grade, and chemotherapy response, which demonstrated good predictive performance for 1-, 2-, and 3-year survival in OC patients. Importantly, by integrating molecular, pathological, and treatment-related factors, this model complements existing clinical tools and may facilitate individualized risk stratification. Given its association with immune features, it also has potential to guide patient selection for immunotherapies or more intensive follow-up strategies, thereby supporting precision oncology in OC.

Although this study constructed a nomogram related to ICD and ER stress for predicting the prognosis of OC, several limitations should be acknowledged. First, the model lacks external clinical validation, which is necessary to confirm its robustness and generalizability. Second, the specific roles of the identified genes in ICD and ER stress remain unclear, and further studies are needed to elucidate their functional roles. Finally, large-scale datasets and prospective clinical studies are required to validate the predictive efficiency of the model.

Conclusion

In this study, a risk score was established by seven-genes related to ICD and ER stress for predicting prognosis of OC. The risk score effectively stratified patients according to survival outcomes. A nomogram integrating the risk score, tumor grade, and chemotherapy response was then constructed, showing good predictive performance for 1-, 2-, and 3-year survival. This nomogram may serve as a practical tool in clinical practice to facilitate individualized risk assessment and support personalized treatment decisions for OC patients.