The RNA-sequencing expression profiles (level 3) and corresponding clinical information of 376 OV tumor samples were downloaded from the TCGA database (https://portal.gdc.com) up to April 2022. Meanwhile, the transcriptome data from 180 normal tissues were downloaded as controls from the GTEx database (https://gtexportal.org/home/datasets), as depicted in Figure 1A. We normalized the gene expression profiles via the limma package in the R software. Based on the set cut-off criteria of |Log2 (Fold Change)| >1 and adjusted p-value < .05, the differentially-expressed genes (DEGs) between OV and normal tissues were identified. Additionally, the RNA-sequencing expression profiles and corresponding clinical information of 111 OV patients were downloaded from the ICGC database (https://dcc.icgc.org/releases/current/Projects) as the validation set. The workflow of the study is shown in Figure 1A.

To identify ferroptosis-related mRNAs of OV patients, we downloaded a total of 182 FRGs (Relevance score ≥ 1) from the Genecards database (https://www.genecards.org/). Through the Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/), we selected FRGs, which are differentially expressed between OV tumor tissues and normal tissues. The multi-gene correlation heatmap was displayed through the “ggstatsplot” package in the R software package. For Kaplan-Meier (K-M) curves of the genes, the p-value and hazard ratio (HR) with 95% confidence interval (95% CI) were generated via Log-rank tests.

The least absolute shrinkage and selection operator (LASSO) regression algorithm with 10-fold cross-validation was performed for the feature selection of FRGs with prognostic value. Multivariate Cox regression analyses were used to build the FRGs-based prognostic model. The “glmnet” R package was used to identify prognostic gene signatures and calculate the risk score of each patient in the datasets based on the signature. For survival analysis, the samples were divided into low-risk and high-risk groups based on the medium cut-off value, and the K-M analysis was used to explore the prognostic significance of the signature. Time-dependent receiver operating characteristic curve (ROC) analysis of 1-year, 3-year, and 5-year survival was analyzed using the “timeROC” R package.

Next, both univariate and multivariate Cox regression analyses were performed to assess the independent prognostic factors to build the nomogram. The forestplot was used to show the p-value and HR (95% CI) of each parameter, using the “forestplot” R package. Based on the selected prognostic parameters, a prognostic nomogram was developed to predict the 1-year, 3-year, and 5-year OS of OV patients in the TCGA dataset, using the “rms” R package. We measured the discrimination of the nomogram model through calibration curves, which could overlay the observed probabilities and nomogram-predicted probabilities with 95% CI. Moreover, we validated the model via the Harrell’s concordance index (C-index) with a 10-fold cross-validation.

OV tissues were obtained from 36 OV patients, who have received cytoreductive surgery, followed by platinum-based chemotherapy, at the Department of Gynecology in Renji Hospital from 5 Marc 2019 to 16 December 2019. OS (overall survival) was measured from the date of initial surgery to the last follow-up or death. Progression-free survival (PFS) was identified from initial treatment to the last follow-up or OV progression, which was assessed by clinical and radiographic evidence. This project was approved by the Ethics Committee of Renji Hospital Affiliated to Shanghai Jiaotong University School of Medicine. All patients provided informed consents, in compliance with the declaration of Helsinki, for the usage of their samples for research purposes.

Total RNA of the tissues was extracted by the Trizol Reagent (T9424, Merk) and then reverse transcribed into cDNA by the RevertAid First Strand Cdna Synthesis Kit (K1622, Thermo Fisher Scientific) following the protocols. Then, we conducted the real-time quantitative reverse transcription-polymerase chain reaction (RT-PCR) analysis through the SYBR Green Master Mix (A25742, Thermo Fisher Scientific) using the QuantStudio™ 7 Flex Real-Time PCR System (Life technologies, United States) following manufacturer’s instructions. Primer sequences were designed as follows: SLC7A11, Forward: 5′—TCA​TTG​GAG​CAG​GAA​TCT​TCA—3′ and Reverse: 5′—TTC​AGC​ATA​AGA​CAA​AGC​TCC​A—3′; ZFP36, Forward: 5′—CCC​AAA​TAC​AAG​ACG​GAA​CT—3′ and Reverse: 5′—GCT​CTG​GCG​AAG​CAC​A—3′; TTBK2, Forward: 5′—ATG​CTC​ACC​AGG​GAG​AAT​GT—3′ and Reverse: 5′—TGC​ATG​ACC​ACG​TAG​TTG​AAA—3′; and GAPDH, Forward: 5′—GGC​AAA​TTC​CAT​GGC​ACC​G—3′ and Reverse: 5′—TCG​CCC​CAC​TTG​ATT​TTG​GA—3′. The comparative expression level was calculated by the 2^−ΔΔCT method. GAPDH was used as an internal control.

To confirm underlying functions and associated high-level genome information of potential genes, we analyzed the data via functional enrichment. We conducted Gene Ontology (GO, including molecular function, biological pathways, and cellular components) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses by Metascape (https://metascape.org), a public-available online tool to analyze DEGs from multiple data sets (Zhou et al., 2019). Stepwise, we processed the DEGs through the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org), a website that could provide screens for human–protein interactions (von Mering et al., 2003), after which the Cytoscape software was used to generate a visualized protein-protein interaction (PPI) network.

To characterize the tumor immune microenvironment, we estimated the abundance of tumor-infiltrating immune cells (TIICs) of each sample through the CIBERSORT algorithm (Newman et al., 2015), a immunological computing method based on a gene expression signature matrix with various marker genes. Original gene expression data from TCGA was normalized before CIBERSORT analysis. Then, at the CIBERSORTx website (https://cibersortx.stanford.edu/), we downloaded the gene signature matrix with interpretation, which outlined 22 subtypes of TIICs. In order to enhance deconvolution algorithm’s accuracy, we accounted for the p-value and root mean squared error of the CIBERSORT. The correlation between the risk score of the signature and immune cells was calculated via the Spearman’s correlation test.

The Pearson’s test was performed to assess the association between the signature and expression of immune checkpoint genes, such as cytotoxic T-lymphocyte-associated protein 4 (CTLA4), programmed cell death protein 1 (PDCD1), and programmed cell death protein 1 ligand 2 (PDCD1LG2), etc. Additionally, based on the RNA-sequencing expression profiles and corresponding clinical data from the TCGA database, we evaluated the potential Immune Checkpoint Blockade (ICB) response. Stepwise, we predicted the potential ICB response of OV patients through the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu), a public-available computational framework which was developed by Jiang et al. (2018) to model tumor immune escape and predict ICB response.

Furthermore, we predict the chemotherapeutic response of each patient based on the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/), the largest public-available pharmacogenomics database. The prediction process was conducted using the “pRRophetic” R package. To identify effective drugs for OV treatment, the samples’ half-maximal inhibitory concentration values (IC50) were obtained from the GDSC database and estimated by ridge regression.

3 All the bioinformatic statistical analyses were implemented by the R software (foundation for statistical computing 2020, version 4.0.3). All the p-values had been passed a two-tailed test, and p-value < .05 was considered statistically significant. Differences between the low-risk and high-risk groups were compared using the Wilcoxon test, and p-values were adjusted through the BH method. Spearman correlation analysis was applied to estimate the correlation between quantitative variables without normal distribution.

The transcriptome data and corresponding clinical characteristics of 376 OV patients were obtained from the TCGA database. Meanwhile, the transcriptome data from 180 normal tissues were also downloaded as controls from the GTEx database. A total of 6,406 DEGs were identified, among which 2,333 genes were upregulated, and 4,073 genes were downregulated in OV, compared to the normal controls (Figures 1B,C). Then, 182 FRGs (Relevance score ≥ 1) were downloaded from the Genecards database for our further analyses, among which 63 FRGs were differentially expressed between tumor tissues and normal tissues, as shown in the Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) (Figure 1D). The enrichment analysis of the 63 differentially expressed ferroptosis-related genes (DE-FRGs) was processed by Metascape (https://metascape.org) (Zhou et al., 2019), which exported the top 20 most significant KEGG and GO pathways in the TCGA-OV cohort (Figure 1E). The pathways were mainly enriched in cellular response to chemical stress, response to inorganic substances, ferroptosis, etc. In order to obtain a PPI network, the DE-FRGs were processed through the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org), a website that could provide screens for human–protein interactions (von Mering et al., 2003) (Figure 1F).

For multiple regression analysis, the LASSO-penalized Cox analysis is a common method to enhance model explicability and forecast accuracy. From the above-mentioned 63 DE-FRGs, six genes (ZFP36, ITGB8, SLC7A11, MYCN, SREBF1, and TTBK2) were filtered to be related to OV patient prognosis through the LASSO regression model (Figures 2A, B). The overview of the function of the six potential DE-FRGs was listed in Table 1 (Helland et al., 2011; Nie et al., 2013; Liao et al., 2015; Cui et al., 2018; Moore et al., 2018; Zhu et al., 2020; Hong et al., 2021). The expression distribution of the six optimal prognostic FRGs in OV tumor tissues and normal tissues was shown in Figure 2C, and the correlation among the FRGs was also presented (Figure 2D). Next, we used the multivariate Cox hazard regression analysis to distinguish prognostic genes, namely SLC7A11, ZFP36, and TTBK2 (Figure 2E). Therefore, the three-gene prognostic signature model was ultimately constructed as follows: Risk score = (−.2084) * SLC7A11+ (.0945) * ZFP36 + (.3619) * TTBK2. The K-M survival curves showed that OV patients with high expression of ZFP36 and TTBK2 suffered worse OS, while those with high expression of SLC7A11 had better OS (Figure 2F).

We focused on the TCGA database as the training set (n = 374) and the ICGC database as the validation set (n = 111). Then, the risk score of every OV patient was calculated via the formula mentioned above. Based on the median cut-off point, we divided OV patients into two groups: high-risk and low-risk in both training and validation sets (Figures 3A, B, top). We also showed the survival status of all OV patients and the expression profiles of the three prognostic genes in high-risk and low-risk groups (Figures 3A, B, middle and bottom). Most of the death cases were mainly distributed in the high-risk group. Additionally, TTBK2 was highly expressed in the high-risk group, while SLC7A11 was highly expressed in the low-risk group. The K-M survival curves showed that the high-risk group suffered worse 1-year, 3-year, and 5-year OS, compared with the low-risk group in the training set (p-value = .00362) and validation set (p-value = .00579) (Figures 3C, D). Additionally, the ferroptosis-associated three-gene prognostic signature showed promising AUC values in the time-dependent ROC analysis for 1-year, 3-year, and 5-year OS ((Figures 3E, F). Taken together, the ferroptosis-associated three genes were prognostic signature for OV patients.

We analyzed the relationship between the ferroptosis-associated 3-gene signature and clinical characteristics (Figures 4A–D). The results indicated that features including age, race, grade, and FIGO stage had no significance with the signature (p-value ≥ .05). Figure 4E showed the distribution of each OV patient, referring to different clinical variables and the risk groups stratified by the 3-gene signature. In addition, we conducted both univariate and multivariate Cox regression analyses to find out whether risk score was an independent prognostic factor for OV patients (Figures 5A, B). Through the multivariate Cox regression analyses, we found that in addition to risk score (p-value < .001), FIGO clinical stage (p-value = .044) and age (p-value = .008) were also confirmed as prognostic factors. Based on the integration of risk score, FIGO clinical stage, and age, we constructed a nomogram of 1-year, 3-year, and 5-year OS probability among OV patients, with the concordance index (C-index) of .6334 (Figure 5C). Calibration curves of the nomogram implied excellent consistency with standard curves between observed and predicted 1-year, 3-year, and 5-year outcomes (Figure 5D). Furthermore, we calculated the nomogram score of every OV patient and divided them into two groups based on the median cut-off point. The K-M survival curves showed that OV patients with high nomogram scores suffered worse OS in both training and validation sets (Figure 5E). Thus, our results indicated that the nomogram model based on the 3 FRGs could predict OV patient prognosis efficiently.

However, considering that the findings of bioinformatics analysis based on public database have many uncertainties, we tried to verify the correctness through experiments on human tissue. We enrolled in 36 OV patients at our institution, with the medium follow-up time of 37.88 (31.48–42.72) months. The clinical features of the involved patients were listed in Table 2. We conducted the qRT-PCR analysis to measure the expression of SLC7A11, ZFP36, and TTBK2 in the OV tissues. The results revealed that higher expression of ZFP36 and TTBK2 was found in OV patients who suffered poor prognosis, while lower SLC7A11 expression were found among them (p-value < .05, Figure 6A).

Stepwise, to better evaluate the significance of ferroptosis-related genes among OV patients, we performed both univariate and multivariate Cox regression analyses for clinical features in relation to OV patient prognosis (Figures 6B,C). The results indicated that SLC7A11, ZFP36, and TTBK2 (p-value = .029, .004, and .015, respectively) were prognostic factors, in addition to FIGO clinical stage (p-value = .028). The K-M survival curves showed that patients with higher expression of ZFP36 and TTBK2 suffered worse OS, while those with high expression of SLC7A11 had better prognosis (p-value < .05, Figure 6D), which was consist with the results of bioinformatics analysis.

Stepwise, we defined the FRGs-related pathways through enrichment analysis among the groups stratified by the 3-gene signature. The KEGG pathway enrichment analysis identified several critical pathways, such as the microRNAs in cancer, PI3K-Akt signaling pathway, MAPK signaling pathway, and others (Figure 7A). As shown in Figure 7B, interestingly, the GO biological process (BP) pathways were mainly enriched in those related to gene silencing, G-protein coupled receptor signaling pathway, response to growth factors, etc. The GO cellular component (CC) pathway analysis identified significantly enriched pathways, including chromatin, terminal bouton, intrinsic component of the plasma membrane, and others (Figure 7C). In Figure 7D, the GO molecular function (MF) pathways were mainly enriched in RNA binding involved in post-transcriptional gene silencing, mRNA binding, DNA binding transcription factor activity, etc. Collectively, pathways related to gene silence and RNA modification were significantly enriched.

Recently, researchers reported that the m6A, a common type of RNA modification, played a critical role in cancer development and progression (Shi et al., 2019). Accordingly, we derived 19 typical m6A-related genes (including METTL3, METTL14, RBM15B, RBM15, MTAP, YTHDC1, YTHDC2, ZC3H13, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPC, HNRNPA2B1, RBMX, ALKBH5, and FTO) from a study on the clinical significance and molecular characterization of m6A modulators of 33 cancer types in the TCGA pan-cancer project (Li et al., 2019). Interestingly, all these m6A-related genes, expected for YTHDC2, HNRNPC, and RBMX, were significant-highly expressed in the high-risk group, compared to the low-risk group (p-value < .05, Figure 7E).

Additionally, we evaluated the landscape of immune infiltration among OV patients through the CIBERSORT deconvolution algorithm in order to determine whether this 3-gene signature was related to the tumor immune microenvironment. Figure 8A summarized the composition of the 22 immune cells infiltrating in OV patients from both low-risk and high-risk groups. Based on the CIBERSORT analysis, 3 out of the 22 immune cell proportions, including activated Myeloid dendritic cells (DCs), plasma cells, and M0 macrophages, were significantly different between the two risk groups. As shown in Figure 8B, activated DCs and plasma cells were downregulated in the high-risk group compared to the low-risk group, while M0 macrophages were upregulated. Correlation analysis between the 22 immune cells implied that Macrophage M1 and CD8+ T cells had the highest positive relationship, with a correlation coefficient of .44 (p-value < .0001) (Figure 8C). Except for the intense negative correlation between all kinds of activated cells and resting cells, Macrophage M2 and Follicular helper T cells had the strongest negative relationship, with a correlation coefficient of .46 (p-value < .0001).

Nowadays, mountains of evidence support the clinical implications of immune checkpoint molecules in the realm of immunotherapy for OV patients. Accordingly, we evaluated the relationships between the risk score and expressions of immune checkpoint molecules. The results implied that CD274, CTLA4, HAVCR2, and PDCD1LG2 were significantly higher in the high-risk group than in the low-risk group (Figure 9A, p-value < .05), suggesting that OV patients in the high-risk group could be more likely to benefit from immunotherapies related to these critical immune checkpoints. Stepwise, we applied the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to predict clinical response toward immune checkpoint blockade (ICB) in two risk groups. Interestingly, we found that OV patients in the high-risk group had greater TIDE score, which represents poorer efficacy of ICB and shorter survival after the ICB treatment (Figure 9B, p-value = .0001).

To assess chemotherapy sensitivity between the two risk groups, we estimated the half-maximal inhibitory concentration (IC50) of eight commonly used OV chemotherapy agents by ridge regression based on the Genomics of Drug Sensitivity in Cancer (GDSC) database (Figure 9C). Our data showed that the estimated IC50 levels of Paclitaxel, Vinblastine, Docetaxel, Sorafenib, and Veliparib in the low-risk group were significantly higher than that in the high-risk group, indicating that the OC patients in the high-risk group were more sensitive to these drugs. However, there was no significant difference in sensitivity to Cisplatin, Bleomycin, and Gemcitabine between the two risk groups.