The processed RNA sequencing profiles (UCSC TOIL RNA-seq) with clinical information were extracted from TCGA(Goldman et al., 2020). There were 416 OC samples and 88 normal samples in our study. After converting Ensemble Gene to Gene Symbol and filtering through genes with low expression, the lncRNAs were identified according to the annotation information in GENCODE database (Harrow et al., 2012). We obtained 288 ferroptosis-related genes (FRGs), of which 108 FRGs were drivers, 69 FRGs were suppressors and 111 FRGs were markers, from FerrDb database (Zhou et al., 2020). In the aggregate, 215 mRNAs associated with ferroptosis were screened via matching it to the filtered mRNAs above. Differential analysis of ferroptosis-related mRNAs and annotated lncRNAs between tumor and normal samples was according to linear regression and Empirical Bayes (p < 0.05, |logFC|>1) (Ritchie et al., 2015). In addition, we used the Benjamini and Hochberg for multiple testing correction to get adjusted p value and assessed the differences at multiple and significance levels (p-adjusted<0.05, |logFC|>1) (Ferreira, 2007; Ghosh, 2019). We calculated the Pearson correlation coefficient between FRGs and lncRNAs to filter the lncRNAs that are significantly related to mRNAs (|r|>0.5, p < 0.001).

According to R’s survival package (Wang et al., 2016), univariate Cox regression was performed on candidate FRLs to screened prognostic lncRNAs. The R-language random function was used to randomly divide the total OC samples (total set, TOS) into two sets, training set (TS) and validation set (VS). Whereafter, using the LASSO Cox regression (Tibshirani, 1997)in the glmnet software package of R (Engebretsen and Bohlin, 2019)and 20 times cross-validation analyses, the optimal combination of FRL markers was screened based on the TS set. We constructed a risk model for OC patients based on formula: Risk score ( RS ) = ∑ β 1 n c R N A × EXP 1 n c R N A

In the formula, ß _lncRNA represented the Lasso prognostic coefficient for optimized FRLs and Exp_lncRNA meant the expression level. In the TS set, the OC samples were classified into high-risk (with a RS higher than or equal to the median value of RSs, n = 104) and low-risk (with a RS lower than the median value of RSs, n = 104) according to the mean value of RSs. In addition, in order to verify the accuracy of the model, RS values of each sample in TOS (high-risk: n = 208; low-risk: n = 208) and VS (high-risk: n = 104; low-risk: n = 104) sets were calculated using the same regression coefficient. The association between the risk groups and the prognosis statuses was evaluated by using the Kaplan-Meier curve method of survival package in R.

A visual nomogram was constructed after univariate and multivariate Cox regression analysis of clinical features and risk groups, and verified by correction curves to determine the accuracy of the risk model.

Based on ssGSEA algorithm, we calculated the enrichment fraction of FRGs in different samples to obtain the ferroptosis score (FS). Based on the median of FSs, the samples were divided into high-score (n = 208) and low-score (n = 208) groups. We compared the differences in survival between patients with high-score and low-score and, in conjunction with the risk scores, we also analyzed survival among patients in different subgroups (low-risk + low-score: n = 104; low-risk + high-score: n = 104; high-risk + low-score: n = 104; high-risk + high-score: n = 104). DAVID was used to analyze the biological functions of potential FRGs (Huang et al., 2009a; Huang et al., 2009b). Differentially expressed genes (DEGs, |logFC| and FDR<0.05) were filtered between the two risk groups on the basis of the limma R package (Clark et al., 2014). Candidate small molecular drugs and mechanisms of action were predicted via uploading DEGs to CMAP database (Lamb et al., 2006; Lamb, 2007).

We used five algorithms (Estimate (Yoshihara et al., 2013), ssGSEA (Jia et al., 2018); (Hänzelmann et al., 2013), CIBERSORT (Newman et al., 2015; Charoentong et al., 2017), MCPcounter (Shi et al., 2020), and xCell (Aran et al., 2017)) to estimate differences in immune microenvironment between risk groups. The expression levels of immune checkpoint genes were extracted, and the expression differences among groups were compared by inter-group T test.

FRG-FRL co-expression (|PCC|>5, p < 0.001) network was constructed and visualized by Cytoscape3.6.1 (Shannon et al., 2003). The miRWalk speculated that the ferroptosis-related mRNAs targeted by the corresponding miRNAs (Dweep and Gretz, 2015). Furthermore, we combined the results of six commonly used databases (Microt4, miRWalk, miRDB, miRanda, Targetscan and RNA22) to obtain the miRNA-FRG relationship pairs. The miRNAs targeted by corresponding FRLs of the risk model were speculated by miranda (Enright et al., 2003). FRLs and FRGs regulated by the same miRNA were defined as ceRNAs mutually. The UCSC data set (Yuan et al., 2019) searched for potential cis-interactions between FRLs and mRNAs transcribed from the same chromosome within 200 kb. Then, further construct a network based on overlapped TFs to detect the trans-regulatory function of FRLs and FRGs (Shu et al., 2019).

Moreover, we assessed the Tumor Immune Dysfunction and Exclusion (TIDE) score through the TIDE database to predict the response to Immune Checkpoint Blockade (Fu et al., 2020). The chemotherapy drugs were extracted from Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013) and the pRRophetic R package is used to calculate the 50% maximum inhibitory concentration (IC50) (Geeleher et al., 2014).

With the approval of the Ethics Committee, specimens were taken from 30 patients. Fresh specimens (≤0.5 cm) were rapidly immersed in RNA Stabilization Solution (GenePharma, China) and incubated the samples overnight at 4°C. The tissue was completely permeated and then transferred to a -80 °C refrigerator for subsequent experiments.

Using ultrasonic crusher, the tissues were cracked by lysate. Total RNA was extracted from OC cells and tissues by using RNA Extraction Kit (TAKARA, Japan). Moreover, cDNA was gained after reverse transcription by using the PrimeScriptTMRT (TAKARA, Japan). PCR was performed by ABI 7500 instrument to evaluate the expression level of the nine FRLs by using TB Green® Premix Ex TaqTM (TAKARA, Japan). Total 1 μL cDNA was used as the template in a final 10 μL PCR reaction volume containing 5 μL TB Green Premix EX Taq, 0.4 μL corresponding primers (GenePharma, China), 0.2ul ROX Reference Dye II and 3.4 μL DEPC water. PCR was run as follows: 95.0 °C for 5s, 40 cycles of 95.0 °C for 5s and 60.0 °C for 34s. The melting curve was obtained at 95.0 °C for 15s, 60.0 °C for 1min and 95.0 °C for 15s. The primer sequences are shown in Supplementary Table S1.

We used R software package, TBtools, and GraphPad Prism for statistical analysis and visualization. Log-rank test was used to evaluate the statistical significance of the difference between the survival curves of two subgroups (Bland and Altman, 2004). Wilcoxon test was used to compare the differences between the two groups. The experiment was made in triplicate, and each experiment was repeated three times. Student’s T test or ANOVA were used for the experimental data. p < 0.05 was considered statistically significant.

A total of 98 differential FRGs (Figure 2A) and 1,145 exegetical lncRNAs were identified (Figure 2B). Enrichment analysis further exhibited that 131 GO and 24 KEGG pathways were enriched based on aforementioned 98 FRGs. The result showed that the identified FRGs were closely related to several important biological processes or pathways, such as cellular response to tumor necrosis factor (Figure 2C), response to hypoxia (Figure 2C), PI3K-Akt/MAPK/HIF-1/p53 signaling pathway (Figure 2D), etc.

548 FRLs were obtained by Pearson correlation analysis (|PCC|>0.5, p < 0.001) and FRL with significant prognosis was mined (p < 0.05). The FRL-signature was constructed by the LASSO Cox analysis of 70 prognostic FRLs. The λ selection diagram is shown in Figures 3A,B.

For 70 prognostic FRLs obtained, the LASSO Cox Regression model from glmnet package and 20-fold cross-validation analysis were further used to screen the optimal biomarker combinations of FRLs. Nine optimal FRLs were filtrated and FRL CTC-246B18.8 was a risk factor with HR > 1, while FRLs AC133644.2, RP1-223E5.4, RP1-313I6.12, RP11-872J21.3, RP3-512B11.3, RP5-1120P11.1, SNHG10, and USP30-AS1 were protective factors with HR < 1 (Figure 3C, Supplementary Figure S1). To further validate the expression levels of the nine FRLs in OC, the qPCR analysis was detected in our collected specimens and cells. The results showed that FRL CTC-246B18.8 was significantly upregulated in OC cells and samples (Supplementary Figures S2,3), while FRLs SNHG10 and USP30-AS1 were downregulated compared with normal control (Supplementary Figures S2,3).

According to the mean value of RSs, the patients with OC from TCGA database were divided into two risk subgroups. The K-M survival curves indicated a worse prognosis in the high-risk group (Figures 3D–F). The univariate and multivariate Cox regression analysis of clinical features and the risk model showed that age, “primary therapy outcome success”, and risk model were independent prognostic factors (Figures 4A,B). A nomogram was further constructed based on age, “primary therapy outcome success”, and risk model (Figure 4C). The calibration curve was drawn to prove the accuracy of the model (Figure 4D).

We calculated the ferroptosis enrichment score of each sample. Then, we found that patients in the high-risk group had a higher ferroptosis score (Figure 5A). Patients with OC were divided into two subgroups based on the median of ferroptosis scores (FSs). The Kaplan-Meier survival curve showed that the OS of patients with the high ferroptosis score was significantly lower (Figure 5B), indicating that both high ferroptosis score and the high-risk score are risk factors (Figure 5C). Heatmap of the associations among the expression levels of nine FRLs, risk scores, ferroptosis scores and clinical features (age, stage, grade, therapy outcome, and tumor residual disease) were displayed in Figure 5D. We found that low expression of CTC-246B18.8 was associated with lower risk score, while the opposite was true for the remaining eight FRLs (Supplementary Figure S4). Low expression of CTC-246B18.8 showed lower ferroptosis score, while RP1-313I6.12, RP1-313I6.12, RP3-512B11.3, and RP5-1120P11.1 did the opposite (Supplementary Figure S5).

The enrichment analysis of functional pathways showed that 91 pathways showed significant differences between the two risk subgroups. The KEGG pathways were sorted according to the p value, and the top 10 were selected for display (Figure 6A). GSEA enrichment analysis showed that the risk groups were associated with tumor-related pathways, such as Hedgehog Signaling Pathway, Homologous Recombination, base Excision Repair, and immune-associated pathway, Antigen Processing and Presentation (Figure 6B). The differentially expressed genes (DEGs) were uploaded to CMAP database and compared with the potential drugs. The small molecule drugs were identified according to the p value, and the top ten were selected for display (Figure 6C).

The matrix score, immune score, estimated score, and relative infiltrate abundance of immune and stromal cells per sample were estimated based on five algorithms which were combined and shown in the heat map (Supplementary Figure S6). Wilcoxon was further used to compare the significance of each cell between the two groups. The results showed that there were statistically significant differences in interstitial fraction among 40 microenvironment cells (Figure 7A).

We compared the expression of seven checkpoint genes. The block diagram of the expression distribution of seven immune checkpoint genes (CD47, CD274, LAG3, HAVCR2, TNFRSF4, SIRPA, and VTCN1) between the two risk groups was shown in Figure 7B. The results showed that the other six genes emerged remarkable differences (p < 0.05), except HAVCR2 (p = 0.66). SIRPA had a higher expression, while the other five genes had a lower expression in the high-risk group. In Figure 7C, TIDE score of OC patients in the low-risk group was lower than that in the high-risk group, suggesting that OC patients with low RSs were more sensitive to ICB treatment.

We matched 20 m6A regulators and saw that the expression levels of FMR1, HNRNPC, METTL16, METTL3, and METL5 were significantly lower in the high-risk group. However, the differences of the remaining m6A regulators between the two groups did not reach a significant level (Figure 7D).

Figure 8 showed the results of nine commonly used OC chemotherapy agents. Our data showed that the IC50 levels of Rucaparib in the low-risk group were significantly higher than that in the high-risk group. Inversely, the IC50 levels of Cisplatin, Paclitaxel, Veliparib, and Vinblastine in low-risk group were significantly lower than that in high-risk group, indicating that the OC patients in the low-risk group were more sensitive to these drugs. However, there were no significant differences in sensitivity between the two risk groups to Bleomycin, Docetaxel, Gemcitabine, and Vinorelbine.

We selected 21 FRLs with significant prognosis to fabricate the FRG-FRL co-expression network (Figure 9A). A total of 105 FRL-miRNA-FRG relationship pairs were obtained, including nine FRLs, 80 miRNAs, and 12 FRGs (Figure 9B). The 12 FRGs included: CEBPG, CISD2, CYBB, GPT2, HMOX1, NCF2, NRAS, OTUB1, SLC2A1, SRC, TAZ, and ULK1. The networks according to the cis and trans interaction were displayed in Figures 9C,D. And, we found that p53 may regulate FRGs and FRLs (CTC-246B18.8-SRC/ULK1, RP11-872J21.3-OTUB1, RP5-1120P11.1-SLC2A1/CEBPG) through trans interactions. These findings may provide a basis for understanding the regulatory mechanisms of FRLs correlated with FRGs.