Identification of immunogenic cell death gene-related subtypes and risk model predicts prognosis and response to immunotherapy in ovarian cancer

Datasets proceeding

EOC expression data were extracted from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. For the TCGA cohort, read count files from each sample were obtained by the GDC portal tool, and genes with a minimum count of zero in all samples were filtered out, followed by conducting the trimmed mean of M values (TMM) method and calculating the normalized factor supplied by edgeR software. Considering the same histopathology, a large cohort GSE9891 was selected for validation which samples were diagnosed as serous or papillary serous carcinoma, and duplicated genes with the distinct probes were merged by their median values following logarithmic transformation. A total of 362 and 255 samples from TCGA (https://portal.gdc.cancer.gov/) and GSE9891 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE9891) were remained, after excluding samples with unrecorded or less than 30 survival days separately. Moreover, thirty-four ICD genes were derived from a previous study by Garg, Ruysscher & Agostinis (2016).

Consensus clustering

Specimens from TCGA set as the training cohort were categorized by k-means clustering method based on Euclidean distance among genes performed by the ConsensusClusterPlus package. This process was iterated 1,000 times with 0.8 proportion of resampling to guarantee consistent outcomes. The optimal cluster number (k = 2) ranging from k = 2 to 10 was determined in the elbow and heatmap plots, defined as ICD-high and ICD-low subtypes.

Screening for differentially expressed genes (DEGs)

The expressed differences of all protein-coding genes between clusters were implemented by the limma package. The filtering thresholds were determined as a corrected P-value by Benjamini Hochberg (BH) less than 0.05 and an absolute log2FC value exceeding 1. In addition, the Wilcoxon rank-sum test was also conducted for specific genes or cells from distinct groups.

Enrichment analysis

Interactions among ICD-related proteins were drawn from the STRING database (Szklarczyk et al., 2021) (https://string-db.org/) and customized by Cytoscape. After conversion from DEG symbols between two ICD subtypes to Entrez IDs applied by org.Hs.eg.db package, clusterProfiler (Wu et al., 2021) software was carried out to gain the functional domains of Gene Ontology (GO) from biological process, cellular component, molecular function, and curated pathways of the Kyoto Encyclopedia of Genes and Genomes (KEGG).

Gene set enrichment analysis

Gene Set Enrichment Analysis (GSEA) offers statistical insights into the biological significance of gene sets (Liberzon et al., 2015). To gain a further understanding of the relevant signaling pathways and molecular mechanisms activated in the ICD-high subtype, GSEA was used to analyze the enrichment differences and the enrichment scores (ES) or normalized ES were calculated with curated c2.cp.kegg.v7.4.symbols.gmt datasets chosen as the reference set, which was available from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/downloads.jsp).

Differences in immune characteristics of the two ICD subtypes and correlation analysis

To elucidate the tumor microenvironment (TME) between two ICD subtypes, two estimation methods, as CIBERSORT (https://cibersort.stanford.edu/), ESTIMATE (Yoshihara et al., 2013), were applied to deconvolute the bulk transcriptome for estimation of the relative or absolute component of 22 immune cell types, and calculate the stromal or immune scores representing the infiltrated extent of stroma or immune cells in tumor tissues, as well as estimate scores, from which tumor purity could be inferred. Subsequently, comparisons were conducted to compare the association of 22 immune cell types and expression levels of immune checkpoints (ICP) (Hu et al., 2021), human leukocyte antigen (HLA) genes (Ju et al., 2021), as well as immunomodulators related to tumor-infiltrating lymphocytes, immune cytolytic activity, and interferon response. Finally, the tumor immune dysfunction and exclusion (TIDE, http://tide.dfci.harvard.edu/) algorithm based on the online platform was employed to measure the response of immunotherapy, thus evaluating the interaction among tumor cells, immune cells, and immune escape mechanisms.

Construction and validation of the prognostic ICD risk signature

First, the prognostic values of ICD cluster-related DEGs were extracted from the univariate Cox regression method. Then least absolute shrinkage and selection operator (LASSO) Cox regression analysis provided in glmnet package was applied performing continuous shrinkage and feature selection to formulate an ICD prognostic risk signature. All differentially expressed ICD genes were considered here, extracting and ranking the optimal lambda values by introducing a penalty term throughout 1,000 iterations, during each of which ten-fold cross-validation was set to derive the best-fit lambda value while minimizing the mean cross-validated error. Finally, the LASSO Cox model was constructed through the formula:

$$RiskScore=\sum _{i=1}^n{k}_i\ast A_i.$$

Here, A_i represents the expression value of the ICD genes, k_i is designated as the regression coefficient, and n denotes the number of optimal solutions derived from the LASSO results. The frequency of genes in each iteration was additionally counted representing their importance. Meanwhile, this ICD risk signature was also validated in the GSE9891 cohort.

Survival analysis

Kaplan-Meier (KM) analysis and log-rank test were performed using the survminer and survival packages to compare overall survival (OS) between the high- and low-risk groups based on the ICD risk signature in TCGA training and GSE9891 validation cohorts. Univariate Cox regression analysis was adopted to determine potential prognostic capacity among clinical pathological variables and the ICD risk signature, following an examination of prognostic independence by multivariate Cox analysis for OS.

Human EOC and normal ovarian epithelial tissues

After obtaining approval from the Medical Ethics Committee of the Second Affiliated Hospital of Harbin Medical University (approval number, YJSKY2022-074) and written informed consent from the patients, tissues from 10 women diagnosed as EOC at the Second Affiliated Hospital of Harbin Medical University were retrieved, of whom six cases also manifested omental metastasis and other four cases showed normal omentum. Additionally, six normal ovarian epithelial tissues were obtained. Inclusion criteria were EOC patients without chemotherapy, radiotherapy, or other neoadjuvant therapy; without other malignant tumors; and showing a clear pathological diagnosis. The collected tissues were excised and immobilized in formalin, cut into thin slices, and embedded in paraffin.

Immunohistochemistry (IHC)

Paraffin sections of ovarian and greater omentum tissues were first dewaxed and hydrated, then subjected to heat-induced antigen retrieval with EDTA (pH = 8.0), and washed five times with distilled water. Sections were blocked with goat serous and incubated with primary antibody (rabbit anti-human NLRP3 antibody (1:500, CY5651; Abways, San Diego, CA, USA)) and murine anti-human IFNG monoclonal antibody (1:1,000, MH45241; Abmart, Shanghai, China) at 4 °C overnight. We used 100 µL of goat anti-rabbit/mouse IgG (Kit-0038; Nakasugi Jinqiao) as a secondary antibody and sections were incubated at room temperature for 30 min. Freshly configured DAB (1:1) was added to stain the tissues, which were then counterstained with hematoxylin. Sections were dehydrated and sealed, and then images were obtained using an Olympus microscope. Three randomly 200× fields-of-view were photographed and preserved. Immunohistochemical semi-quantitative expression analysis was performed using a double-blinded method. A total of 10 ovarian cancer tissues, six normal ovarian epithelial tissues, six metastatic omental, and four normal greater omental cases were analyzed. ImageJ (NIH, Bethesda, MD, USA) was used to measure integrated optical density (IOD) and area (Area). The average absorbance (AOD) was calculated as AOD = IOD/Area.

Statistical analysis

Analyses and visualization of high-throughput experimental data were all processed in R software (version 4.3.0) and Linux. Experimental data by IHC were analyzed and plotted by GraphPad Prism 9.0, and the results were expressed as the mean ± standard deviation. All the results were repeated at least three or more times. The t-test or Wilcoxon rank-sum test was used for comparison between the two groups, and the one-way ANOVA or Kruskal–Wallis test was used for comparisons of more than two variables. All statistical studies used P < 0.05 or the corrected P value < 0.05 as the threshold for statistical significance.

Expression patterns and consensus clustering of ICD-related genes in EOC

ICD-related genes in this study were compiled from a large-scale meta-analysis and their expression profiles in EOC samples compared to normal ovarian tissues were first investigated. Elevated expression of ICD genes, including IL10, CALR, IL1B, BAX, FOXP3, IFNG, ATG5, CXCR3, CD8A, CD8B, NLRP3, CASP1, CD4, PRF1, HSP90AA1, IFNB1, PDIA3, CASP8, MYD88, TNF, and IL6, was observed in EOC (Fig. 1A). Then protein-protein interaction (PPI) network was examined and visualized via STRING database and CytoScape software, elucidating tight connections among ICD genes (Fig. 1B). A total of 362 eligible OC patients from the TCGA cohort were used as the training cohort to identify potential clusters based on expression patterns of ICD-associated genes. Through the consensus clustering method, two OC-ICD clusters were identified based on Euclidean distance of which cluster C2 exhibited higher expression of ICD-associated genes compared to C1 (Figs. 1C–1E). Hence C2 was labeled as ICD-high subtype and C1 as ICD-low subtype (Fig. 1F). After excluding samples less than 30 survival days, the KM plot revealed different prognoses between the two ICD subtypes (log-rank P = 0.026, Fig. 1G). The ICD-high subtype had a better prognosis than the low subtype.

Exploration of DEGs and functional enrichment analyses between different ICD subtypes

The identification of DEGs and signaling pathways specific to each subtype aided in understanding the molecular mechanisms affecting prognosis. Based on linear fitting and empirical Bayes estimation provided by limma, 1440 DEGs were identified, with 1,220 genes exhibiting higher expression in the ICD-high subtype and 220 genes exhibiting lower (Fig. 2A). GO and KEGG enrichment analyses revealed that overexpressed genes in the ICD-high subtype primarily participated in immunological activities, including lymphocyte-mediated immunity, activation of immune response via cell surface receptor signaling pathways, interaction between cytokines and cytokine receptors, and humoral immune responses (Figs. 2B, 2C). Furthermore, 17 DEGs were identified as ICD-associated genes with elevated expression in the ICD-high subtype (Fig. 2D). GSEA was utilized to compare the contribution of different gene sets to DEGs between the two subtypes. The results revealed significant activation of several immune-associated signaling pathways in the ICD-high subtype, including the toll-like receptor pathway, Fc gamma R-mediated phagocytosis, the T cell receptor signaling pathway, and leukocyte transendothelial migration (Fig. 2E).

Characterization of the tumor immune microenvironment landscape in ICD-high and ICD-low subtypes

Growing evidence has indicated that ICD could elicit specific antitumor immune responses. It was further demonstrated substantial activation of immune signaling pathways, particularly in the ICD-high subtype. Analyses of the TME composition here showed the ICD-high subtype had increased stromal, immune, and ESTIMATE scores, and reduced tumor purity compared to the ICD-low subtype (P < 0.0001, Fig. 3A). Using the CIBERSORT deconvolution method combined with the LM22 feature matrix, the relative proportions of 22 immune cells were further assessed between the two subtypes in the TCGA training cohort (Fig. 3B). This revealed the ICD-high subtype had notably elevated percentages of M1 and M2 macrophages, resting NK cells, plasma cells, and CD8⁺ T cells (Fig. 3C). Correlation analysis revealed strong associations among infiltrating immune cells, demonstrating complexity in OC tissues. For instance, naive B cells, plasma cells, and CD8⁺ T cells positively correlated with M1 macrophages, while memory B cells negatively correlated with M0 and M1 macrophages, naive B cells, and plasma cells (Fig. 3D). These infiltration of 22 immune cell types and correlations were further validated using the GSE9891 dataset (Fig. 3E). Furthermore, comparing the expression pattern of the 25 ICP and 17 HLA genes showed virtually upregulated patterns in the ICD-high subtype, which might mean strong immunogenicity in the high subtype (P < 0.0001, Figs. 3F, 3G).

Construction and validation of the prognostic ICD risk signature

The potential capacities of ICD-associated DEGs in predicting prognosis were further explored. Demonstrated by univariate Cox analysis two distinct genes NLRP3 and IFNG out of 17 target ICD genes from all DEGs were identified and showed significant association with OC prognosis (Fig. 4A). Based on machine learning algorithms, the 17 ICD-related DEGs were considered into LASSO models by introducing 1,000 random iterations, each performed by 10-fold cross-validation to identify optimal lambda values and genes, counting the frequency of models in each iteration. The results indicated a model comprised of four distinct genes appearing with the highest frequency (Fig. 4B). Components of these four-gene models were further examined conforming to a consistent composition (Figs. 4C, 4D). Specifically, the ranking frequency of IFNG, NLRP3, FOXP3, and IL1B was settled after 1,000 iterations, suggesting sequential significance for prognosis among these ICD genes (Fig. 4E). A risk score model was thus developed based on these genes, with coefficients derived from the LASSO algorithm, fitting as risk score = −0.02139749 * FOXP3 − 0.30089077 * IFNG + 0.05682286 * NLRP3 + 0.02748733 * IL1B. Correlations between risk score and survival status were investigated in the TCGA training cohorts, of which two risk groups were categorized based on the median risk value, and results showed patients with higher risk scores preferred more deaths and shorter survival time (Fig. 4F). Furthermore, a significant trend also revealed poorer OS in the high-risk group compared to the low-risk group (log-rank P = 0.0012, Fig. 4G), as well as in the GSE9891 validation cohort (log-rank P = 0.014, Fig. 4H).

Univariate and multivariate Cox regression analyses were utilized to assess the predictive performance of the ICD risk model compared with other clinical signatures. Both tumor remnant and the risk model were identified as significant prognostic indicators associated with survival, as a high ICD risk score was indicated as a risk factor (P = 0.002, Fig. 4I). Multivariate Cox assessment showed that the ICD risk score was an independent predictor among other clinical parameters (P = 0.005, Fig. 4J).

Exploration of the association between ICD risk signature and tumor immunotherapy

Since the significance of ICD in anti-tumor immune response has been acknowledged, the association between the ICD risk score and the tumor immunotherapy was examined. Results indicated a link between the two indices, with higher risk scores negatively correlated with CD8 cells, follicular helper T cells, and M1 macrophages in the TCGA cohort (P < 0.05, Fig. 5A), as evidenced in the GSE9891 cohort (P < 0.05, Fig. 5B).

The TIDE score was utilized to evaluate the clinical efficacy of immunotherapy across different ICD risk-score groups. As the TIDE score increased, susceptibility to immune evasion was also escalated, suggesting a potential reduced benefit from immune-checkpoint inhibitor (ICI) treatment. Results revealed a higher TIDE score in the high-risk group compared to the low-risk group (P < 0.001, Fig. 5C), whereas the microsatellite instability (MSI) score and the T-cell dysfunction score were higher in the low-risk group (P < 0.0001, Figs. 5D, 5E), against the higher T-cell exclusion score in the high-risk group (P < 0.0001, Fig. 5F). Moreover, excluding samples with less than one survival month, there was no significant difference in risk scores between the immunotherapeutic non-responder and responder groups (Fig. 5G); however, samples in the immunotherapeutic non-response group proved higher ICD risk scores after excluding samples with less than three survival months (P < 0.05, Fig. 5H), suggesting that immunotherapy might be more beneficial for the low-risk group.

Validation for expression of ICD-related genes in clinical tissues

Distinct expression patterns of specific ICD genes were verified by immunohistochemical analyses in clinical tissues. NLRP3 and IFNG, the top two genes in the risk model, were selected and statistically elevated in OC tissues compared to normal ovarian tissues (P < 0.05). Both genes also showed an increasing tendency in metastatic vs. non-metastatic omental tissues, with a statistically significant difference in NLRP3 expression (Figs. 6A–6D). The experiment results confirmed the upregulation of NLRP3 in OC and metastatic omentum, enhancing understanding of the related characterization of ICD genes involved in OC progression.