For ovarian cancer prognostic analysis, we collected 106 patients’ H&E diagnostic WSIs with corresponding clinic data from TCGA-OV (https://portal.gdc.cancer.gov/projects/TCGA-OV) via the National Cancer Institute GDC Data Portal. The inclusion criteria herein consist of three aspects: 1) The quality of WSIs was assessed by an experienced clinic doctor; 2) retaining only one WSI for each case; 3) the case contained both clinical information and WSIs. As a result, 90 cases (60 uncensored patients and 30 censored ones) were incorporated to obtain a prediction model. Additionally, we exploited the five-cross validation method to train the model.

Herein, we utilized an open-source tool, CLAM WSI-analysis toolbox (Lu et al., 2021), to implement the segmentation and feature extraction tasks for each WSI. In this scenario, each original WSI consists of four levels of different resolutions. First, we segmented the tissue region of interest from the level 0 with the highest resolution in each WSI. Second, the whole slide image for each patient was split into M patches with 256 × 256 pixels without overlap. Third, a pre-trained deep network ResNet50 (trained on the ImageNet dataset) was exploited to extract feature tensors by feeding patches into it. Finally, the third block in the ResNet50 model was selected to output the feature tensor with 1,024 dimensions. Additionally, normalized gene expression was measured as Transcripts Per Kilobase Millions, and we processed the genetic mutation data of the TCGA dataset using the R package “maftools”.

For each gigapixel WSI, it is a challenging task to provide pixel-level labels via human labor. The goal of survival data analysis is to train a predictor for hazard probability in a time interval based on plenty of image patches. Thus, we designed a weakly-supervised deep learning architecture as shown in Figure 1A, and the representative images were shown in Figure 1B.

We split the entire pipeline into two parts. In the first part, we extract M feature tensors X from M patches with 256 × 256 pixels via the pre-trained model. As we knew, the tissue region of interest varies with each WSI. In the second part, an attention-based network is designed to assign weights for each feature tensor as the input of the subsequent fully connected layer (FC). A set of learnable weight parameters are obtained for all patches feature for one case.

Due to the size of the dataset being slightly small, a two-layer network is designed to build the prediction layer. The overall network is optimized to predict the four hazard probabilities in four intervals via the loss function. The patient-level risk score is calculated by summing the four hazard scores. For subsequent analysis, the high-group and low-group are divided according to the median of all risk scores.

To realize the survival analysis from patient-level data, we divide evenly the overall patient survival time into four intervals [t _i ,t _i+1), i = 0,1,2,3 according to the uncensored cases. Y j ∈ 0,1,2,3 denotes the ground truth for the jth case. The subscript j denotes the jth patient. The prediction layer outputs the corresponding hazard probability as shown in Eq. 1. f h a z a r d r | x j = P T j = r | T j ≥ r , x j (1) where x j denotes the feature tensor for the jth case. Given T j ≥ r and x j , the conditional probability f h a z a r d r | x j denotes the probability that the time period T j (in months) ends at time r. The survival function f s u r v is described in Eq. 2. f s u r v r | x j = P T j > r | x j = ∏ u = 1 r 1 − f h a z a r d u | x j (2)

The loss function is exploited as shown in Eq. 3 (Zadeh and Schmid, 2021). L = L c e n s o r e d + L u n c e n s o r e d = − c j ⁡ log f s u r v Y j | x j − 1 − c j log f s u r v Y j − 1 | x j ∙ f h a z a r d Y j | x j (3)

To balance the differences between the uncensored cases and censored counterparts, a hyper-parameter α is utilized in the final loss function   L s u r v . L s u r v = 1 − α L + α L u n c e n s o r e d (4)

The HRD score is calculated as the sum of telomeric allelic imbalance (TAI), large-scale state transitions (LST), and loss of heterozygosity (LOH) scores (Shi Z. et al., 2021). HRD scores are derived from research by Thorsson et al. (2018). HRD+ was defined as a high HRD score (threshold >42 score).

The relative infiltration level of immune cell types was quantified via single sample gene set enrichment analysis (ssGSEA) by the “GSVA” R package (Hanzelmann et al., 2013). GSEA was performed to assess related pathways.

Chemotherapeutic response prediction for OC samples was conducted in R by using the “oncoPredict” package (Maeser et al., 2021) from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013). The ridge regression model was applied to evaluate the half maximal inhibitory concentration (IC50).

We performed a univariate analysis based on clinic parameters and risk scores. Afterward, multivariate Cox regression was conducted using the significant prognostic variables (p < 0.05). The nomogram was generated using the R package “rms”. We used calibration curves to test the consistency between predicted and actual survival rates. A time-dependent Receiver operating characteristic (ROC) curve was also used to assess the predictive accuracy of the nomogram. In addition, the Decision Curve Analysis (DCA) was used to demonstrate the advantage of the prediction curve using the R package “ggDCA”.

The statistical significance for variables with non-normal distribution was analyzed using the Wilcoxon rank sum test. The comparison between two groups of variables with normal distribution was estimated using an unpaired Student’s t-test. Non-parametric correlation analyses were conducted based on Spearman’s rank correlation coefficient. The prognostic analysis was performed by the Kaplan-Meier analysis method, and the log-rank test was used to evaluate significant differences. All statistical analyses were conducted using Python software (version 3.7) and R software (version 4.1.3). p < 0.05 was considered statistically significant.

Overall 90 cases picked from the original TCGA-OV dataset were divided randomly into training (72 cases) and test (18 cases) datasets, respectively. We initialized the network parameters with “nn.init ()” and adopted the Adam solver with a momentum of 0.9 for the training process. Batch size, epoch number, and α are initialized as 1.40, 0.31, respectively. Additionally, the initial learning rate is set to 0.0002 and decayed with a Cosine Annealing schedule (Loshchilov and Hutter, 2017). To verify the effectiveness of the proposed network model proposed, we conducted a 5-fold cross-validation experiment on one NVIDIA GeForce RTX 3090 GPU. We utilize Pytorch based on Python 3.7 to implement all training and test tasks.

C-index was originally proposed to evaluate predictions for binary responses (Harrell et al., 1996). As an evaluation metric utilized widely, it is herein utilized to evaluate the network performance. Figure 2A demonstrated that the C-index overall increased in the process of training. The C-index in the cross-validation experiment varies from 0.5096 to 0.6053 and the corresponding mean value was 0.5789. Figure 2B showed that the loss function reduces overall in 5-fold cross-validation results.

To verify further the effectiveness of the deep survival network, the KM curve was employed to analyze the prediction results. The high- and low-risk cohort was classified via the median value of the sorted hazard value associated with each patient. The KM cure of overall survival and recurrence-free survival were plotted as shown in Figures 3A, B. In addition, the p-value was calculated to analyze these two distributions based on the log-rank test. p-value of 0.00845 demonstrated that there was a significant difference between the high- and the low-risk score group The clinical characteristics between the high and low risk score subgroups, including age, stage, grade, and residuals were presented in Table1 to provide a clearer understanding of the sample distribution. Detailed risk score and clinical characteristics were shown in Supplementary Table S1.

Platinum-based chemotherapy is the essential treatment for OC. Patients with HRD+ (HRD score> 42) and BRCA1/2 mutation are more beneficial from chemotherapy, thus, we evaluated the prognostic applications of the risk score in HRD+ and HRD-subgroups. In the TCGA-OV cohort, the HRD + subgroup displayed a significantly better OS than HRD-group (p < 0.0001, Figure 4A). 61 overlapped patients with WSIs and HRD scores were detected (Figure 4B). The HRD + subgroup harbored mainly high-risk score population, while the HRD-group revealed an opposite result (Figure 4C), which suggested that the risk score may work in a different way from the method based on HRD in predicting survival outcomes. Moreover, survival analysis demonstrated a significantly better OS in patients with high-risk score than that with low-risk score in HRD + subgroup (p = 0.013; Figure 4D). Interestingly, there was no significant difference in HRD-group comparing the OS of the two risk populations (p = 0.92, Figure 4E). ROC curve analysis is used to assess the sensitivity and specificity of prediction models and validate the results of risk prediction values. The time-dependent ROC curve proved the reliable performance of risk score in HRD + subgroup (2-years AUC = 0.743, 3-years AUC = 0.718, 5-years AUC = 0.775, Figure 4F). To present the HRD score, FIGO staging, survival status, and risk score calculated by WSI as a unified system, Sankey diagram was constructed to describe the relationship between these features (Figure 4G).

The relationship between risk score and mutation landscape has been assessed in OC patients. In both high- and low-risk score groups, the top 20 mutated genes were presented in Figures 5A,B. TP53 (93%), USH2A (20%), and TTN (17%) exhibited the most frequent mutations in the high-risk score group, while the highly mutant genes in the low-risk score group were TP53 (79%), TTN (18%) and AHNAK (11%). Moreover, the high-risk score group genes were enriched in HOMOLOGOUS RECOMBINATION, OXIDATIVE_PHOSPHORYLATION, and RIBOSOME pathway (Figure 5C), as well as the low-risk score group genes were enriched in FOCAL ADHESION, JAK STAT SIGNALING PATHWAY, and ECM RECEPTOR INTERACTION pathway (Figure 5D). The mutation (Supplementary Tables S2, S3) and GSEA data (Supplementary Table S4) were shown in Supplementary Material. We also tried to analyze the correlation between the risk score and TMB, but got a negative result (Figure 5E). Since the high antigenicity induced by tumor mutations can recruit a large number of immune cells, we studied the relationship between risk and TIDE score, to evaluate the predictive value of risk score in immunotherapy outcomes. Unfortunately, although the risk score was correlated with TIDE (p = 0.0091, Figure 5F), the spearman correlation coefficient was relatively low. These results were consistent with the poor response of ovarian cancer to immunotherapy. Next, we analyzed the relationship between risk score and immune cell infiltration by ssGSEA. Several types of immune cells, such as the central memory CD4 T-cell, central memory CD8 T-cell, and effector memory CD8 T-cell, were found significantly less recruited in high-risk group tumor environment (Figure 5G). The risk score may be a supplemented as an indicative tool to further investigate immune cell infiltration in ovarian cancer.

We attempted to identify whether the risk score could be applied to predict the sensitivity of response to chemotherapies using the GDSC database. The results revealed that the low-risk score group had a lower half maximal inhibitory concentration of BMS-754807, doramapimod, JAK1_8709, JQ1, NU7441, RO3306, SB216763, and WZ4003, while the high-risk score group had a lower half maximal inhibitory concentration of cisplatin, leflunomide (Figure 6). The p38MAPK inhibitors ralimetinib have been shown to play an anti-cancer role in ovarian cancer (Campbell et al., 2014). While Doramapimod also showed potent anti-inflammatory effects as p38MAPK inhibitors (Schreiber et al., 2006), perhaps our study will explore a new alternative for the application of Doramapimod in ovarian cancer. Similarly, RO3306 (Yang et al., 2016), SB216763 (Kaltofen et al., 2020) were also shown to play an anti-cancer role in ovarian cancer. Moreover, the combination of JQ1 and cisplatin helped ovarian cancer-bearing mice survive (Yokoyama et al., 2016). However, previous study showed that NU7441 could induce resistance to PARP inhibitor in BRCA1-defective cells (McCormick et al., 2017), and BMS-754807 combined with carboplatin/paclitaxel was observed resistance in ovarian carcinosarcoma of patient-derived xenograft (Glaser et al., 2015). Therefore, we should actively explore new strategies for more drug combinations to avoid the resistance in ovarian cancer. Our model provided possibility for novel pathways of drugs. Also, we hope that more pre-clinical models will prove our predicted results.

Based on the available clinic features, Cox regression analyses were conducted to identify the possibility that risk score was an independent prognostic factor for OS. The univariate Cox regression analysis revealed significant associations between risk score and OS (HR = 0.48, 95% CI = 0.266–0.86, p = 0.0138, Figure 7A). When other confounding factors were corrected, multivariate Cox regression analysis proved the risk score of WSIs was an independent predictor of prognosis (HR = 0.505, 95% CI = 0.275–0.928, p = 0.028, Figure 7B).

Following the results of the Cox regression analysis, we have further developed a nomogram incorporating two independent prognostic factors (risk score and age) to offer a quantitative method for estimating 2-, 3-, and 5-years survival rates of OC patients ( Figure 7C). In addition, calibration plots showed that the nomogram was comparable to an ideal model ( Figure 7D ). At 2, 3, and 5 years, the AUC of the risk score was 0.605, 0.578, and 0.680. Nomogram’s accuracy of prediction over 2, 3, and 5 years was significantly higher, at 0.680, 0.613, and 0.745, respectively ( Figure 7E ). DCA result also indicated our nomogram had a promising potential for clinical application ( Figure 7F).