The gene expression data in counts format and corresponding clinical data of TCGA was downloaded from UCSC Xena (https://xenabrowser.net/datapages/?dataset= TCGA-OV.htseq_counts.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443), which was then transformed from the raw counts to the count-per- million (cpm) value for later DRPS construction. The combined data of TCGA and GTEx in DESeq2 normalized counts format was downloaded from UCSC Xena (https://xenabrowser.net/datapages/?dataset=TCGA-GTEx-TARGET-gene-exp-counts.deseq2-normalized.log2&host=https%3A%2F%2Ftoil.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) The somatic mutation data of TCGA-OV was downloaded from GDC Data Portal (https://portal.gdc.cancer.gov/). The GSE9891 and E-MTAB-386 datasets used for validation were downloaded from https://www.ncbi.nlm.nih.gov/geo/ and https://www.ebi.ac.uk/biostudies/, separately. Totally, we included 374, 285, and 129 patients with OC from TCGA database, GSE9891, and E-MTAB-386 datasets.

The eighteen disulfidptosis-related genes were obtained from Liu et al. (2023), including four suppressor hits (SLC7A11, SLC3A2, RPN1, and NCKAP1) and fourteen core actin-related proteins (INF2, CD2AP, PDLIM1, ACTN4, MYH9, MYH10, IQGAP1, FLNA, FLNB, TLN1, MYL6, ACTB, DSTN, and CAPZB). The disulfidptosis-related prognostic model for OC was developed by two steps in the TCGA-OV cohort: First, we conducted least absolute shrinkage and selection operator (LASSO) method with the “glmnet” package (version 4.1–7) to select the robust prognosticators among 18 candidate genes using cpm format data (transformed from the raw counts) (Friedman et al., 2010). The best lambda value in the LASSO regularization was chosen by 10-fold cross-validation. Second, Cox model with forward stepwise regression was performed to build the DRPS and the corresponding prognostic model with “survival” package (version 3.5–5) (Therneau and Grambsch, 2000). Patients were stratified into high- and low-risk groups based on whether their risk scores were higher than the median of the whole cohort (high risk and worse prognosis) or less than the median of the whole cohort (low risk and better prognosis). The model was subsequently validated in two independent external datasets: GSE9891 and E-MTAB-386. The significance of the prognosis between these two groups was evaluated through the application of Kaplan-Meier survival analysis and the p-value was calculated using the log-rank test. In order to enhance the accuracy of prediction, we incorporated the FIGO stage, grade, and residual disease, along with the disulfidptosis-related risk score, to construct a nomogram using “survival” package (version 3.5–5) and “regplot” package (version 1.1) (Therneau and Grambsch, 2000). The efficacy of this nomogram was subsequently validated through the Receiver-operator characteristic (ROC) curve using “timeROC” package (version 0.4) (Blanche et al., 2013).

Mutation signature of high- and low-risk groups was analyzed by the R package “maftools” (Mayakonda et al., 2018). Six phenotype characteristics were included in the comparison, including survival status, FIGO stage, histological grade, tumor residual size, homologous recombination deficiency (HRD) status, and germline BRCA mutation status (gBRCAmut). Since SLC7A11, SLC3A2, RPN1, and NCKAP1 had been verified as suppressor hits of disulfidptosis, we treated these four genes as a geneset and calculated the gene set variation analysis (GSVA) scores per patient for it to show the disulfidptosis level (Hänzelmann et al., 2013).

Differential gene expression analysis between high- and low-risk groups stratified by DRPS followed the limma-voom pipeline using raw counts data (Ritchie et al., 2015). Significant differentially expressed genes (DEGs) were defined as genes whose FDR were less than 0.05 and absolute value of log2 fold change greater than 1. Enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathways were conducted separately on the significant DEGs between risk groups stratified by DRPS and genes included in the DRPS to identify functional pathways by the “clusterProfiler” package (version 4.8.1) (Wu et al., 2021). Gene set enrichment analysis (GSEA) in HALLMARK gene pathways was also performed with the “clusterprofiler” package. Annotated clusters were ranked according to Group Enrichment Score and p-value. Pearson correlation coefficients between genes in DRPS were calculated by “Hmisc” package (version 5.1-0). Sensitivity analysis of the DRPS was performed using copula-based methods (Yeh et al., 2023).

We calculated the DEGs between OC and normal ovary using limma method. OC samples and normal ovary samples were selected from the combined TCGA and GTEx expression matrix in DESeq2 normalized counts format for further comparison. DEGs between OC and normal ovary were identified as absolute value of log2 fold change >2 and FDR <0.001 from the result of limma.

The immune cell infiltration in TME was evaluated by TIMER 2.0 (http://timer.cistrome.org/) (Li et al., 2017), CIBERSORT (Chen et al., 2018), xCell (Aran et al., 2017), EPIC (Racle and Gfeller, 2020), quanTIseq (Finotello et al., 2019), MCPcounter (Becht et al., 2016). Since there’s no available OC cohort with immunotherapy response, the published markers as well as TIDE score were used to predict the potential immune checkpoint inhibitor (ICI) response (Jiang et al., 2018).

Single-cell RNA sequencing data was obtained from EMBL-EBI (https://www.ebi.ac.uk/biostudies/arrayexpress) (Qian et al., 2020), including 44532 cells from 5 patients with primary OC (7 OC samples and 3 normal samples). Data quality control was conducted with the primary filtering condition of nFeature_RNA >200 & nFeature_RNA <6000 & percent. mt < 25 & percent. HB < 10 and DoubletFinder (version 2.0.3). Clustering was performed using the “Seurat” package (version 4.3.0.1) at a revolution of 0.5 (Butler et al., 2018). According to published studies, we used a set of marker genes to determine cell types: EPCAM, KRT7, KRT8, KRT17, SPRR3 for epithelial/tumor cells; CD3E, CD3D, TRBC1/2, TRAC for T cells; LYZ, CD86, CD68, FCGR3A for myeloid cells; CD79A/B, JCHAIN, IGKC, IGHG3 for B cells and plasma cells; CLDN5, FLT1, CDH1, RAMP2 for endothelial cells; DCN, C1R, COL1A1, ACTA2 for fibroblasts; TAGLN, CNN1 for smooth muscle cells (Dinh et al., 2021).

The drug sensitivity was predicted by half-maximal inhibitory concentration (IC50) according to GDSC data (https://www.cancerrxgene.org/) using the “oncoPredict” package (version 0.2) (Maeser et al., 2021). The potential targeted drug for risk groups was assessed by Drug-Gene Interaction database (DGIdb) (https://www.dgidb.org/) (Freshour et al., 2021).

We downloaded the immunohistochemical (IHC) staining images of OC and normal ovaries from the Human Protein Atlas (HPA) (https://www.proteinatlas.org/). ImageJ (version 1.8.0) was used to analyze the integrated optical density (IOD) and average optical density (AOD).

R version 4.2.2 was used for statistical analyses. We utilized WilcFoxon test to evaluate the statistical significance for continuous characteristics, and Chi-square test or Fisher’s exact test for categorical characteristics as appropriate. Two-tail p-value <0.05 was set as statistically significant. The R script used in this study was uploaded in Supplementary Materials, R script.

The flowchart of this study was shown in Figure 1. The baseline characteristics of the datasets showed that epithelial (serous) was the main histologic type and grade 3 was the main pathologic grade in all the three cohorts, consistent with global real-world reports (Supplementary Table S1) (Kuroki and Guntupalli, 2020). Six genes (MYL6, PDLIM1, ACTN4, FLNB, SLC7A11, and CD2AP) were selected by the LASSO and then used as the input of multivariate Cox regression analysis to build the DRPS (Figures 2A, B). According to Figure 2C, FLNB was an independent risk factor for overall survival (OS), while SLC7A11 and CD2AP were the independent protective factors. Kaplan-Meier curves for OS illustrated that the prognosis differed markedly between the high- and low-risk groups stratified by the DRPS not only in the training set (p-value < 0.0001) but also in two validation sets (p-value = 0.0059 for GSE9891, p-value = 0.019 for E-MTAB-386) (Figures 2D–F).

Besides, we calculated the risk score of the DRPS by multivariate cox regression. The final disulfidptosis-related prognostic model was developed using the risk score, histological grade, FIGO stage, and residual disease, which was then generated as a nomogram (Figure 2G). As Figure 2H showed, the risk score was an independent risk factor for OS prediction after adjusting those clinicopathological factors. ROC curves further demonstrated that the disulfidptosis-related prognostic model achieved satisfactory performance in the training dataset and validation datasets (Figures 2I–K).

Top ten genes with the highest frequency in simple nucleotide variation (SNV) and copy number variation (CNV) were displayed in Figure 3A (high-risk group) and Figure 3B (low-risk group). None of these genes were found to have a significant difference between groups, which indicated that the difference in tumor aggressiveness between groups was more associated with the heterogeneity at the transcriptomic-level rather than the genomic-level. For clinicopathological factors, high-risk group had a higher mortality rate compared with low-risk group (Figure 3C); while FIGO stage, histological grade, residual disease, gBRCAmut, and HRD showed no significantly different distribution between groups (Figures 3D–H).

We calculated the GSVA score of disulfidptosis in each sample using the gene set which included SLC7A11, SLC3A2, RPN1, and NCKAP1 (Figure 3I). As the result, high-risk group showed a statistically lower GSVA score of disulfidptosis compared with low-risk group, suggesting a relatively low level of disulfidptosis in high-risk group. Among the four genes used in the GSVA score of disulfidptosis, the expression of SLC7A11 and SLC3A2 were significantly lower in high-risk group (Figure 3J).

Figures 4A, B showed the top DEGs between high- and low-risk groups. High expression of IGF2, IGLON5, MEGF8, HSPG2, ALDH1A2, NCCRP1 was associated with high-risk group, while high expression of SLC7A11, PLA2G12A, SST, MMP10, KRT14, PAX2 was associated with low-risk group (Figures 4A, B; Supplementary Table S2). The result of GSEA showed that multiple pathways related with malignant biological behavior of tumors were remarkably activated in high-risk group, including hypoxia, epithelial mesenchymal transition (EMT), TGF-beta signaling, WNT-beta-catenin signaling, and angiogenesis (Figures 4C, D); while pathways associated with immune response, including interferon alpha response, interferon gamma response, and oxidative phosphorylation, were enriched in low-risk group (Figure 4E).

GO enrichment analysis revealed that pathways like integrin binding, connective tissue development, and sulfur compound binding were associated with the disulfidptosis-related risk stratification, which were related to the formation and function of cytoskeleton as well as the disulfide bond (Figure 4F). KEGG enrichment analysis indicated that pathways like PI3K-Akt signaling were enriched in high-risk group, and it is worth noting that these enriched pathways have been demonstrated to be closely related to cancer development (Figure 4G).

Furthermore, correlation analysis of the DRPS showed that most of the genes in the DRPS had a significant correlation between pairwise (Supplementary Figure S1). GO/KEGG enrichment analysis showed that DRPS was mainly associated with the actin related pathways. In addition, the pathways related to sulfur transport were also significantly enriched, indicating that genes included in DRPS might participated in disulfidptosis in OC (Supplementary Figure S1). Sensitivity analysis of the DRPS revealed that good separation could be found between high- and low-risk groups in both the training and validation sets (Supplementary Table S3).

Correlation coefficients between the disulfidptosis-related risk score and various immune cells were calculated (Figure 5A). Interestingly, macrophages, neutrophils, and cancer associated fibroblasts (CAFs) were positively correlated with the disulfidptosis-related risk score of patients with OC in more than one algorithm, indicating that high-risk group exhibited more infiltration of these immune cells and these cells promoted the immunosuppressive microenvironment. Conversely, activated myeloid dendritic cells and CD4 memory T cells were negatively correlated with the disulfidptosis-related risk score, indicating that tumors in low-risk group tended to be activated myel infiltrated by these kinds of immune cells. The relative abundances of various immune cells were also estimated between high- and low-risk groups (Figure 5B), which is consistent with the results presented in Figure 5A.

Moreover, we explored biomarkers which have been widely accepted for ICI response prediction. CAF score in high-risk group was significantly higher than in low-risk group (Figure 5C); in contrast, MSI (microsatellite instability) expression signature was significantly lower in high-risk group (Figure 5D). Besides, TIDE score was significantly higher in high-risk group than in low-risk (Figure 5G), in which TIDE dysfunction scores in high-risk group was significantly higher (Figure 5E); whereas TIDE exclusion scores were higher in high-risk group than in low-risk group (Figure 5F) though without statistically significant differences. In total, these biomarkers suggested that high-risk group might have a lower response rate for ICI than low-risk group.

To further elucidate the expression of the genes in DRPS in cell level, we mapped these six genes on the UMAP plot of a public single-cell RNA dataset of OC after unsupervised clustering and annotation (Figures 6A, B). The expression bubble plots of marker genes for cell type clustering were displayed in Supplementary Figure S2. As the plot presented, the major types of cells showed distinct marker genes expression signatures (Supplementary Figure S2B), which is consistent with the previous published researches (Dinh et al., 2021; Jia et al., 2023).

MYL6 was commonly expressed in all the seven major types of cells, which had a significantly higher expression level in fibroblasts, endothelial cells, and B cells/plasma cells in tumor compared with normal tissue (Figure 6D). PDLIM1 and ACTN4 were mainly expressed in endothelial cells, fibroblasts, smooth muscle cells, and tumor/epithelial cells (Figure 6C). Besides, these two genes showed markable higher expression in the cell types mentioned about in tumor than normal tissue (Figure 6D). FLNB, SLC7A11, and CD2AP were lowly expressed in almost all major cell types (Figure 6C), but significant expression difference was also found in fibroblasts between tumor and normal tissue.

IC50 of 198 drugs in GDSC database in the TCGA-OV cohort was predicted using the “oncoPredict” package (Supplementary Table S4). IC50 of epirubicin, stauroporine, navitoclax, and tamoxifen were significantly lower in high-risk group than in low-risk group (Figures 7A–D), indicating that these drugs might be the potential treatment for high-risk patients. According to the prediction of DGIdb, riluzole (Figure 7E) as the interacted drug of SLC7A11 was the potential targeted drug for patients in high-risk group.

We took the intersection of the DEGs between high- and low-risk groups, DEGs between OC and normal ovary (Supplementary Table S5), as well as the DRPS. Finally, two genes (SLC7A11 and ACTN4) were selected as key genes for further research (Figure 7F; Supplementary Table S6).

We quantitated IHC staining images from the HPA database using ImageJ software. Since no IHC staining images of SLC7A11 could be found, only other five genes were analyzed for expression at the IHC level. The results revealed significant differences in expression of ACTN4 (p-value = 0.0394 for IOD; p-value = 0.3597 for AOD) according to Figure 7G and S3. CD2AP (p-value = 0.1142 for IOD; p-value = 0.4096 for AOD), FLNB (p-value = 0.9516 for IOD; p-value = 0.2837 for AOD), MYL6 (p-value = 0.0864 for IOD; p-value = 0.2330 for AOD), PDLIM1 (p-value = 0.9545 for IOD; p-value = 0.0522 for AOD) showed no significant differences in protein expression between tumor and normal tissues (Figure 7G; Supplementary Figure S3).