Common gene expression data were obtained from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. In TCGA dataset, RNA sequencing data for gene expression were downloaded from UCSC Xena. For the three datasets of GSE26712, GSE32062 and GSE18521, RNA sequencing data were downloaded from the easyGEO dataset. Among these, TCGA dataset was used as the training set, and the GEO dataset was used as the validation set. The details of the datasets are presented in Table 1, including number of patients, overall survival (OS), tumor stage, and tumor grade.

Immunity-related genes (IRGs) were obtained from the immunology and ImmPort databases. A comprehensive gene list (involving 1793 IRGs) was downloaded. Ferroptosis-related genes (FRGs) were obtained from the FerrDb database and the related literature [23]. The FerrDb database is the first database that provides updated data of regulators and markers connected with ferroptosis and disease. Finally, 267 FRGs were obtained, including 259 FRGs from the FerrDb database, and 8 genes from the literature [23] as a complement.

According to the data of the TCGA training set and the GSE26712 validation set, the IRGs and FRGs were screened using univariate Cox analysis. The intersection of the genes in the two datasets was considered as the prognosis-related genes (p $\le$ 0.05). The least absolute shrinkage and selection operator (Lasso) method was applied to calculate the coefficient ($\beta$) to screen prognostic genes in TCGA dataset for reducing the number of genes in the risk score model. Eventually, the multivariate Cox regression analysis was applied to construct the risk score model. The calculation formula was as follows: Risk score = ($\beta$A $\times$ gene A expression) + ($\beta$B $\times$ gene B expression) + $\cdot$$\cdot$$\cdot$ + ($\beta$N $\times$ gene N expression).

In this formula, $\beta$A, $\beta$B…$\beta$N represent the coefficients of genes. Genes from which the risk score was constructed were defined as signature genes. The patients were then divided into high- and low-risk groups using the median risk score as the cut-off point. Survival curves were then plotted using Kaplan-Meier from R package (version 3.5.3). Accuracy was verified by time-dependent receiver operating characteristic (ROC) curves.

The signature genes in each sample were divided into the high- and low-expression groups according to the median expression in TCGA dataset. Kaplan-Meier analysis was applied to analyze the influence of signature genes on survival. The expression of signature genes in cancer and paracancerous tissues was searched from the Gene Expression Profiling Interactive Analysis (GEPIA) database and verified in the GSE26712 and GSE18521 datasets. At the same time, the Human Protein Atlas database was applied to detect the protein expression of signature genes in cancer tissues.

To further validate the expression of signature genes at the mRNA level, RT-PCR assays were performed using the tissue samples. Tumor tissues and paracancerous tissues were collected from 8 patients with OV. The clinical data of the patients are provided in detail in Table 2, mainly including age, tumor size, stage, grade, tumor metastasis, alcohol, smoking and family history. In addition, the increase of CA125, HE4 and CA199 suggests the possibility of malignancy.

The inclusion criteria for the patients with OV were as follows: (1) Patients were diagnosed with OV for the first time; (2) patients had not undergone other treatments prior to diagnosis; (3) patients had no other malignancy; (4) patients had no other autoimmune disease; (5) the age of the patients was 18 to 70 years. The exclusion criteria were as follows: (1) Patients had other malignancies; (2) patients received other treatments prior to surgery; (3) patients had incomplete clinical data; (4) patients had a history of cancer. The Ethics Committee of The Second Affiliated Hospital of Xi’an Jiaotong University approved this study (2021241). Patient informed consent was obtained.

RNA was extracted from the tissue samples. The SYBR-GreenⅠReal Time PCR method was used to detect the changes in gene expression. Relative gene expression values were determined using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{ct}}$ method (compared to healthy controls). 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{ct}}$ $>$1 and 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{ct}}$ 1 represent up- and downregulation, respectively. ACTB was used as an internal reference gene.

Based on the DGIdb database (https://dgidb.org/), those drugs related to signature gene were screened out, which may provide new perspectives for OV diagnosis, treatment and research.

For the purpose of further investigating the association between clinical features of OV and risk score, univariate and multivariate Cox analysis were first applied to determine the independent prognostic factors of the risk score different from other clinical factors. After testing for collinearity, the independence of risk was also verified. All independent prognostic parameters and associated clinical parameters were applied to construct prognostic nomograms through stepwise Cox regression models to predict the overall survival (OS) of patients with OV. A nomogram was used to forecast the OS. In addition, differences in risk score in the various clinical subgroups were compared.

The relative abundance of each cell infiltration in the TIME of OV was quantified using single-sample Gene Set Enrichment Analysis (ssGSEA). In the study by Charoentong, gene sets that marked each immune cell type infiltrating the TIME, which is enriched in multiple human immune cell subtypes, were obtained [24, 25]. The stromal score, tumor purity and ESTIMATE score of each patient with OV were calculated using ESTIMATE from R package (version 3.5.3). The Wilcoxon test was applied to compare the differences in immune cell infiltration, stromal score, tumor purity and the ESTIMATE score.

Differences in biological processes were investigated using GSVA between the high- and low-risk groups. “C2.cp.kegg.v7.2.Symbols” was downloaded from the Molecular Signature Database. The R package “limma” was applied to calculate significantly enriched pathways. A false discovery rate (FDR) 0.05 was considered as statistically significant.

The R package “limma” was applied to identify differentially expressed genes (DEGs) between the high- and low-risk groups under the screening thresholds of FDR 0.05 and $|$log${}_2$(FC)$|$$>$0.5. The Database for Annotation, Visualization and Integrated Discovery (DAVID) was applied to perform enrichment analysis, and the screening threshold was FDR 0.05.

First, the “maftools” package in R was utilized to calculate the TMB score. Tumor samples were divided into the high- and low-TMB groups based on the median of TMB. The survival curve between risk and TMB was analyzed using the Kaplan-Meier method. In addition, a correlation analysis between the risk score and TMB was performed.

For the purpose of exploring the predictive power of tumor risk score in patients treated with immunotherapy, the expression of multiple immune checkpoints was assessed in the high- and low-risk groups and their correlation with the risk score was determined. Based on immunotherapy cohorts, the responses of patients were analyzed in the high- and low-risk groups. The difference in the risk score in the response and non-response groups was analyzed.

Based on the clinical data of OV cases in TCGA and the GEO datasets, univariate Cox proportional hazards regression analysis was used to screen genes associated with immunity and ferroptosis for survival (p 0.05). The Lasso method was used to screen variables to prevent over-fitting and build a prognostic model. Survival curves for risk score were analyzed using Kaplan-Meier analysis, and the 1-, 3-, and 5-year time ROC curves of the model were plotted using the R package “timeROC”. Univariate and multivariate Cox analyses were used to determine whether risk was an independent prognostic factor distinct from other clinical factors. Statistically significant differences were assessed using the two-tailed t-test in the RT-PCR experiments. The ssGSEA algorithm was used to quantify the relative abundance of each cell infiltration in the OV TIME, and the stromal score, tumor purity and the ESTIMATE score were calculated for each patient with OV using the ESTIMATE algorithm through the R package “ESTIMATE”. In addition, based on the Wilcoxon test, the differences in various immune cells, stromal scores, tumor purity and the ESTIMATE score between the high- and low-risk groups were compared. In addition, the correlation between key genes and immune cells was explored based on Pearson’s correlation analysis. GSVA enrichment analysis was performed using the R package “GSVA” and based on the R package “limma” was used to calculate differentially expressed pathways (FDR 0.05). DEG analysis was performed using the R package “limma” between the high- and low-risk groups (FDR 0.05 and $|$log${}_2$(FC)$|$ $>$0.5). Multiple testing correction was performed in the GSVA and GO enrichment analysis in the high- and low-risk groups. An FDR 0.05 was considered statistically significant. DAVID was used for the enrichment analysis of DEGs between the high- and low-risk groups (FDR 0.05). TMB was calculated using the “maftools” package in R.

To investigate the impact of immunity and ferroptosis genes on the prognosis of patients with OV, univariate Cox regression analysis was applied in TCGA and the GSE26712 datasets with a screening criterion of p $\le$ 0.05; a total of 220 and 118 genes were identified. After taking the intersection of the genes in the two datasets, 12 prognostic genes were finally identified (Fig. 1), including C5AR1, COLEC12, EGFR, GZMB, IGF2R, ISG20, PPP3CA, RB1, STAT1, TAP1, TRIM27 and TSHR.

Fig. 1.

Screening of prognostic genes in OV. (A) Forest diagram of 12 prognosis-related genes in TCGA training set. (B) Forest diagram of 12 prognosis-related genes in the GEO validation set.

Given the possible interaction between immunity and ferroptosis, 12 prognostic genes were used for the construction of the prediction model of OV. The optimal $\lambda$ (lambda) was obtained by the minimum partial likelihood deviation. Finally, 10 signature genes were obtained, including 8 immunity-related genes (C5AR1, GZMB, IGF2R, ISG20, PPP3CA, STAT1, TRIM27 and TSHR), 1 ferroptosis-related gene (RB1) and 1 immunity and ferroptosis-related gene (EGFR). The immunity and ferroptosis-related risk score model was constructed, as shown in Fig. 2A,B. The formula is as following:

Risk score = (C5AR1 $\times$ 0.167018542) + (EGFR $\times$ 0.024944649) + (GZMB $\times$ –0.114825674) + (IGF2R $\times$ 0.223153693) + (ISG20 $\times$ –0.145695011) + (PPP3CA $\times$ 0.137679588) + (RB1 $\times$ 0.176152429) + (STAT1 $\times$ –0.157629949) + (TRIM27 $\times$ –0.446534567) + (TSHR $\times$ 0.14287582)

Fig. 2.

Construction of risk model in TCGA training set. (A) Change trajectory of each independent variable. The horizontal axis and vertical axis represents the log value of the independent variable lambda and the coefficient of the independent variable, respectively. (B) LASSO Cox regression with 10-time cross-validation for optimal parameter ($\lambda$) selection, in which the vertical dashed lines showed minimum $\lambda$ value and 1 $\times$ standard error $\lambda$ value, respectively. (C) Risk score distribution map, the dotted line indicates the optimal cut-off value of the risk score. All OV patients are divided into high and low risk groups. In survival status, the blue dots represent the surviving OV patients and the red dots represent death, and the corresponding risk gene cluster expression profiles were visualized as a heatmap. (D) Survival Curve, blue represents low risk score, red represents high risk score. (E) The yellow, blue and red lines represent 1, 3 and 5 years of the ROC curve, respectively.

The influence of the risk score obtained from the construction of 10 immunity and ferroptosis-related genes on OS was further investigated. The risk score was calculated using a formula (detailed results are presented in Supplementary Table 1), the samples were divided into the high- and low-risk groups based on the median risk score. The low-risk group had a higher survival proportion, as shown in Fig. 2C. The OS was significantly lower in patients in the high-risk group (Fig. 2D). The AUC values at 1, 3 and 5 years were 0.639, 0.645 and 0.716, respectively (Fig. 2E). In TCGA training set, the risk score had good predictive power for the OS of patients.

To further evaluate the robustness of the risk score in predicting the OS of patients with OV, two datasets (GSE26712 and GSE32062) were applied to test the prediction capacity of the risk score for OS. First, the risk scores of the tumor samples in the two datasets were calculated separately using the same formula. The samples were divided into the high and low-risk group based on the median risk score (Supplementary Fig. 1A,D). The low-risk group had a higher OS rate. The OS of the patients in the high-risk group was significantly lower (Supplementary Fig. 1B,E). In the GSE26712 dataset, the AUC values were 0.713, 0.73 and 0.731 at 1, 3 and 5 years, respectively (Supplementary Fig. 1C). In the GSE32062 dataset, the AUC values were 0.643, 0.620 and 0.629 at 1, 3 and 5 years, respectively (Supplementary Fig. 1F). In the validation set, it was further demonstrated that the risk score model had a better predictive ability for patients with OV.

The results of survival analysis with the risk score revealed that the OS of the high-risk group was significantly lower than that of the low-risk group. Based on the expression values of signature genes in each sample of the TCGA training set, the samples were divided into the high- and low-expression groups according to the median expression. The influence of signature genes on survival was analyzed using the Kaplan-Meier method. The C5AR1, EGFR, RB1 and TSHR genes had no significant effect on survival (p $>$ 0.05), while patients with a high expression of GZMB, ISG20, STAT1 and TRIM27 had a longer OS. In addition, patients with a low expression of IGF2R and PPP3CA had a longer OS (Supplementary Fig. 2).

Since there were no data available for the control group in TCGA dataset, the expression of signature genes in cancer and paracancerous tissues was searched from the Gene Expression Profiling Interactive Analysis (GEPIA) database. STAT1 and GZMB were highly expressed in patients with OV, while EGFR was expressed at low levels in patients with OV. No significant differences were found for the remaining genes examined (Fig. 3). Subsequently, the expression of 10 signature genes was verified in the GSE26712 and GSE18521 datasets. GZMB, ISG20 and TSHR were significantly upregulated in patients with OV, while EGFR, PPP3CA and RB1 in GSE26712 were significantly downregulated in patients with OV. In the GSE18521 dataset, C5AR1 was significantly expressed at low levels in patients with OV, while GZMB, ISG20, STAT1, and TSHR were significantly highly expressed in patients with OV (Supplementary Fig. 3). Of note, GZMB was significantly highly expressed in patients with OV in all data analyzed.

Fig. 3.

Signature genes expression of C5AR1, GZMB, IGF2R, ISG20, PPP3CA, STAT1, TRIM27, TSHR, RB1, and EGFR in the GEPIA database. Red and gray represent the patients with tumor group and normal control group (T: tumor; N: normal control), respectively. ${}^{*}$p 0.05.

At the same time, the Human Protein Atlas database was applied to analyze the protein expression levels of STAT1, GZMB and EGFR (Supplementary Fig. 4A). Based on the analysis using the UALCAN database, no significant expression difference was observed for STAT1; however, EGFR and PPP3CA proteins were expressed at low levels in patients with OV (Supplementary Fig. 4B–D).

In addition, RT-PCR detection indicated that GZMB, STAT1 and ISG20 were markedly upregulated in patients with OV (Supplementary Fig. 4E). The OS analysis with individual signature genes revealed that a high expression of the GZMB, STAT1, ISG20 and TRIM27 genes, and a low expression of the IGF2R and PPP3CA genes corresponded to a longer OS rate. However, the analysis of signature genes in cancer and paracancerous tissues revealed that the GZMB, STAT1 and ISG20 genes still had higher expression levels, and that the EGFR and PPP3CAgenes had lower expression levels in cancer tissues.

Based on the DGIdb database, drugs related to anti-tumor and immunotherapy for 10 signature genes were screened out (Supplementary Fig. 5), including 32 drugs and 4 genes (EGFR, PPP3CA, RB1 and STAT1). Among these, cisplatin is widely used in the treatment of various solid tumors [26]. Sunitinib is a novel, multi-targeted oral treatment for cancer. Etoposide, docetaxel and decitabine are cell cycle-specific antitumor agents. Topotecan is indicated for metastatic advanced OV that has failed chemotherapy. Paclitaxel, vorinostat and other drugs are associated with OV treatment [27, 28].

Univariate and multivariate Cox analyses determined the association of the risk score with clinical signatures. Age and risk could be used as independent prognosis factors (Fig. 4A). The independence of risk was verified in the validation set (Fig. 4B). A nomogram was constructed to predict OS probabilities by combining risk and age with significance for independent prognostic indicators in a multivariate Cox analysis (Fig. 4C,D). Each factor was allocated proportionally to its risk contribution to survival. The calibration curve indicated that the combined model (nomogram) revealed a high accuracy in OS at 1 year. In addition, differences in risk scores across clinical subgroups were compared. The risk scores significantly differed by age, status and stage (Fig. 4E).

Fig. 4.

Association of risk score with clinical features. (A) Multivariate Cox analysis in TCGA training set. (B) Multivariate Cox analysis in the GSE32062 validation set. (C) Nomogram of clinical features and risk score. (D) Calibration plot of 1-, 3-, and 5-years. (E) Risk score differences in age, status, stage and grade subgroups. ${}^{*}$p 0.05; ${}^{**}$p 0.01 and ${}^{****}$p 0.0001; ns, not significant.

The infiltration status of 23 immune cells in TCGA dataset was evaluated using ssGSEA. The majority of the immune cells in the high-risk group had more infiltration (Fig. 5A). However, there was no corresponding survival advantage for patients in the high-risk group, and the OS of patients in the high-risk group was significantly lower (Fig. 2D). Epithelial-mesenchymal transition (EMT) 1, EMT2 and EMT3 was significantly higher in the high-risk group subtypes (Fig. 5B). Moreover, the ESTIMATE score, stromal score and tumor purity for each patient with OV was calculated using ESTIMATE. The ESTIMATE score and stromal score were significantly higher in the high-risk group, although tumor purity was low (Fig. 5C–E).

Fig. 5.

Comparison of the immune microenvironment in the high- and low-risk groups. (A) Difference in the infiltration degree of 23 types of immune cells. Red and blue represent high-risk and low-risk groups, respectively. (B) Difference in infiltration degree of EMT, green and pink represent high-risk and low-risk groups, respectively. (C) Difference in ESTIMATE score, orange and blue represent high-risk and low-risk groups, respectively. (D) Difference in stromal score, brown and blue represent high-risk and low-risk groups, respectively. (E) Difference in tumor purity, purple and green represent high-risk and low-risk groups, respectively. ${}^{*}$p 0.05; ${}^{**}$p 0.01; ${}^{***}$p 0.001 and ${}^{****}$p 0.0001; ns, not significant.

Furthermore, the ferroptosis status was also predicted based on suppressor genes (SOFs) and activators (DOFs) of ferroptosis in the literature [29]. Among the 16 repressors of ferroptosis, CD44 and HMOX1 were significantly highly expressed in the high-risk group, while HELLS, HSF1, HSPA5, NQO1, OTUB1 and SLC7A11 were significantly expressed at low levels. In the 22 ferroptosis activators, ANO6, EGFR and ZEB1 were significantly highly expressed in the high-risk group, and CHAC1, ELAVL1, PEBP1 and SAT1 was significantly expressed at low levels. Taken together, the findings indicated that ferroptosis may be more pronounced in the high-risk group (Supplementary Fig. 6).

Significantly enriched pathways were calculated using the R package “limma” and screened on the condition of FDR 0.05. A total of 103 metabolic pathways were screened out. The top 10 pathways were selected from the high- and the low-risk groups, respectively. Some metabolic pathways, such as ECM receptor interaction, focal adhesion and axon guidance were more active in the high-risk group; while in the low-risk group, proteasome, DNA replication, spliceosome, protein export, aminoacyl TRNA biosynthesis and other metabolic pathways were more active (Supplementary Fig. 7A).

The cancer tissue samples were grouped according to the risk score in TCGA dataset and DEG analysis was performed. Finally, 158 DEGs were screened, of which 102 upregulated genes and 56 downregulated genes were identified (Supplementary Fig. 7B). DAVID was used to analyze the function of the DEGs under the screening threshold of FDR 0.05. GO enrichment analysis revealed that the DEGs were mainly enriched in the extracellular region under cellular component (CC). As regards molecular function (MF), the DEGs were mainly enriched in collagen binding. In terms of biological process (BP), the DEGs were significantly enriched in the type I interferon signaling pathway, as shown in Supplementary Fig. 7C.

It was considered meaningful to explore the inner link between the risk score and TMB in order to elucidate the genetic characteristic of each subgroup. Kaplan-Meier analysis of the survival curve revealed that the high-TMB group had a longer survival rate. The increase in the TMB level led to an improvement in survival, while the low-risk and high-TMB group survived the longest (Fig. 6A,B). There was no significant difference in TMB between the high- and low risk groups (Fig. 6C). Furthermore, there is no correlation between risk score and TMB (Fig. 6E). A high TMB induces antitumor immune cell activation and improves the prognosis of patients with OV.

Fig. 6.

Differences in TMB between high and low-risk groups. (A,B) Survival curves of the different groups. (C) Violin plot of TMB. (D) Bar chart of proportional distribution. (E) Correlation between TMB and risk score.

The expression of immune checkpoints, including CD274, CTLA4, KLRC1, LAG3, PDCD1 and TIGIT was assessed in the high- and low-risk groups, followed by the evaluation of the relevance between these immune checkpoints, and risk score and signature genes. In the low-risk group, CD274, CTLA4, KLRC1, LAG3 and PDCD1 were significantly highly expressed. GZMB exhibited a positive association with each immune checkpoint, while the risk score was negatively associated with each immune checkpoint (Fig. 7A,B).

Fig. 7.

Immune checkpoints and anti-PD-1/L1 immunotherapy. (A) Differences in immune checkpoints between high- and low- risk groups. Red and blue represent high-risk and low-risk groups, respectively. (B) Correlation analysis of immune checkpoints. (C) Survival curves of high- and low- risk groups. (D) Differences in risk scores between different anti-PD-1 clinical response groups. (E) Proportion of patients responding to anti-PD-1/L1 immunotherapy in high- and low- risk groups. (F) Accuracy of risk score in differentiating the PD/SD and CR/PR groups. (G) Correlation of risk score with clinical response to anti-PD-1/L1 immunotherapy. ${}^{*}$p 0.05; ${}^{**}$p 0.01; ${}^{***}$p 0.001 and ${}^{****}$p 0.0001; ns, not significant. SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response.

Immunotherapy represented by the blocking of programmed cell death ligand 1 (PD-L1) and programmed cell death 1 (PD-1) has become a key breakthrough in the treatment of cancer. The present study then investigated whether the risk score can forecast the response of patients to immune checkpoint blockade treatment based on an immunotherapy cohort. In the anti-PD-L1 cohort, patients demonstrated significant clinical benefits and a notable prolonged survival was observed in the low-risk group (Fig. 7D). The therapeutic advantage was confirmed in the low-risk patients (Fig. 7D–F). Finally, the ability of the risk score to distinguish between the disease progression group (PD/SD) and disease remission group (CR/PR) treated with anti-PD-L1 immunotherapy was examined using the ROC curve (Fig. 7G). Overall, the findings of the present study strongly suggest that the risk score is significantly associated with the tumor response to immunotherapy. The risk score established herein may contribute to the prediction of the response to immunotherapy.