End-to-End Integrative Segmentation and Radiomics Prognostic Models Improve Risk Stratification of High-Grade Serous Ovarian Cancer: A Retrospective Multi-Cohort Study

Study design and participating cohorts

In this multi-national retrospective study, we used data from 3 cohorts of primary HGSOC (aged ≥18 years). In the Hammersmith Hospital (HH) cohort (n=211) clinical data (including data related to fresh frozen tissue, imaging, and clinical annotations) were obtained from patients treated at HH between June 2004 and November 2015 (Table 1). The HH dataset was split 7:3 into training (n=147), and HH validation (n=64) and used in all machine learning and deep learning models to predict OS, as well as the nn-U-Net segmentation model.

For assessment of our models, we used data from two published benchmark datasets. The first is from The Cancer Imaging Archive (TCIA) cohort (n=73). The second test dataset consisted of 323 cases from Kliniken Essen-Mitte (KEM) and was previously used in (9) to validate the radiomics prognostic vector identified by Lu and colleagues (8).

We used the Radiomics Quality Score (RQS) (12) and Transparent Reporting of a Multi-variable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD) guidelines (https://www.tripod-statement.org/) (Appendix A Tables 1-2) for reporting the development and validation of the prediction models.

Data Collection and Procedures

A retrospective cohort study was conducted with ethical approval for the analysis of human data, which was obtained from the Hammersmith and Queen Charlotte’s & Chelsea Research Ethics Committee (approval 05/QO406/178) and the Kliniken Essen-Mitte Research Ethics Committee (informed consent was waived). The study was carried out in accordance with the ethical standards of the relevant committees and with the principles of the 1964 Declaration of Helsinki and its subsequent amendments.

The HH training and validation data consisted of patients with HGSOC who received treatment at the Hammersmith Hospital (HH) between June 2004 and November 2015. These patients were selected based on the availability of fresh frozen tumour tissue samples and preoperative CT images. The HH and TCIA cohorts’ patient demographics, surgical and tumour related data were collected retrospectively from medical records and the multidisciplinary team (MDT) notes. Demographic, surgical, and tumour-related data were collected retrospectively from medical records and the multidisciplinary team. With the HH cohort, management of all patients and the indications for surgery were discussed within a multidisciplinary team as per the UK National Health Service (NHS) guidelines. All operations were performed through a midline laparotomy by a specialised dedicated multidisciplinary team within a maximal effort approach aiming to achieve total macroscopic tumour clearance. Standard surgical procedures included peritoneal cytology, extrafascial hysterectomy, bilateral salpingoophorectomy and infra-gastric omentectomy. When indicated, additional procedures, such as dissection of macroscopically suspicious pelvic and paraaortic lymph nodes, bowel resection, splenectomy, diaphragmatic stripping/resection and/or partial resection of other affected organs (e.g., urinary bladder, liver/liver capsule, pancreas, lesser sack) were performed to achieve optimal tumour debulking. No systematic pelvic and paraaortic lymph node dissection was performed routinely in the absence of suspicious bulky lymph nodes (<1 cm). Ninety-seven percent of patients were treated with a platinum-based chemotherapy mainly in a combination regimen with paclitaxel or as monotherapy in isolated cases. Patients were regularly evaluated at the end of their treatment for evidence of disease recurrence. Clinical examination and CA-125 assessment (if the preoperative value was elevated) were performed every 3 months for the first 2 years and then 6-monthly. A CT/MRI-scan was ordered if the above examinations revealed any pathology. An isolated CA-125 increase was not regarded as a recurrence.

In the KEM dataset all consecutive HGSOC patients who underwent primary debulking surgery in between January 2011 and July 2018 (9). Patients were eligible for inclusion if they were confirmed to have HGSOC histology and an evaluable portal venous CT through the primary ovarian tumour mass prior to undergoing upfront cytoreductive surgery, i.e., before any systemic chemotherapy. Any suspicious bulky visible lesions (≥1 cm) from the CT scan were included in the segmentation process. Patients were not eligible if no primary tumour mass was visible due to the initial absence of adnexal mass or due to surgical excision or if the patient has undergone preoperative chemotherapy. All patients in the KEM dataset were to subsequently receive platinum-based combination chemotherapy unless contraindications applied such as poor performance status. Maintenance regimens and availability of clinical trials differed between the KEM and the primary HH study setting, mainly due to funding and licensing differences between the European Medical Agency (EMA) and the UK National Institute for Health and Care Excellence (NICE UK). Oncologic follow-up of the KEM data was performed according to the German follow-up recommendations, mainly symptom-guided and based on clinical and ultrasonographic examination in combination with CA-125 measurement in most patients; initially 3-monthly for the first 3 years and then 6-monthly. A CT or MRI scan was performed if the above examinations revealed any pathology. Isolated CA-125 increase was not regarded as a recurrence. For the HH cohort, follow-up patterns for patient care were similar. Patients were routinely evaluated at the end of their treatment for evidence of disease recurrence. Clinical examination and CA-125 assessment (if the preoperative value was elevated) were performed every 3 months for the first 2 years and then 6-monthly. Even though a CT/MRI scan was ordered if the above examinations revealed any pathology, no routine ultrasonographic examinations were performed at follow-up in asymptomatic patients. Isolated CA-125 increase was not regarded as a recurrence.

The primary outcome measure was OS, defined as the time from the date of surgery until the date of death or last observation. Staging was determined using the FIGO criteria for ovarian epithelial carcinoma (13). In the HH and TCIA dataset optimal debulking was defined as postoperative residual disease less than 10 mm, as this criterion was applied to most of the retrospective patients. Primary chemotherapy resistance was defined as stable disease, partial response, or progressive disease during first-line chemotherapy. Residual disease was dichotomised as either total macroscopic tumour clearance (tumour-free) or the presence of any macroscopic postoperative residual disease (non-tumour-free). A subset of HGSOC patients from TCGA was used as a validation cohort. Preoperative CT images and clinical and histological data for these cases were obtained from the Cancer Imaging Archive and the UCSC Cancer Browser, respectively.

CE-CT scans, Segmentations and Radiomics Collection

As patients were referred to the cancer centre from a network of cancer units, CE-CT scans were acquired at multiple institutions using different manufacturers and different imaging protocols. The primary tumour masses, of the HH, TCIA, and KEM datasets, were segmented separately by experienced radiologists using ITK snap (Version 3.2, 2015). All segmentations were then checked in consensus with two experienced radiologists specializing in ovarian cancer imaging (GA and AR, with 6- and 16-years’ experience at the time of HH, and TCIA segmentation). The KEM dataset was manually segmented 3 years after HH and TCIA segmentations.

As primary lesions alone harbour prognostic information (8), we did not include secondary lesions in our analysis. As described previously, the entire primary adnexal mass volume (both cystic and solid components) was included in the analysis (8,9). If both adnexa were involved, then both were included in the analysis either as two separate segmentations or as a single segmentation if the mass was confluent. We segmented the entire primary mass including cystic and solid components but excluded ascites. The segmentations only included tissue that was considered highly likely to be cancer by the expert reader. Areas of doubt on CT were not included in any segmentations. In aggregate, inclusion criteria related to the CT images were as follows: primary adnexal mass visible, portal venous phase CT through the adnexal mass, no previous surgical or medical treatment for ovarian cancer. Exclusion criteria related to images were non-contrast or arterial phase CT with no portal venous phase, adnexal mass not included on CT, previous surgery for resection of an adnexal mass, neoadjuvant chemotherapy, the significant artefact for example from metal prostheses that precluded meaningful segmentation of adnexal mass. In accordance with the Image Biomarkers Standardization Initiative (14) CE-CT scans were resampled to isotropic voxel dimensions of 1×1×1 mm.

We developed and compared several different CE-CT, convolutional neural network (CNN) and radiomics based models, independently and in combination with traditional clinical biomarkers, namely age, stage and residual disease. All continuous variables were scaled and mean-centred using training dataset statistics. All categorical variables were ‘one-hot-encoded’.

Copy Number, Transcriptomics and Proteomics Data Collection

Copy Number, transcriptomics and proteomics data used in this study were previously reported in Lu and colleagues radiomics prognostic vector study (8). For each tumour in the study, one frozen tumour piece was placed into a tube containing 500 μl RLT buffer from RNeasy kit (QIAGEN) and one Retsch 6 mm steel core bead. Tubes were placed into well adapters of a Tissuelyser II (QIAGEN), and tissues were lysed at 15 Hz for 2 min. Tubes were centrifuged briefly and 320 μl was removed for subsequent RNA extraction using the RNeasy kit (QIAGEN) according to the manufacturer’s instructions. RNA concentrations were quantified using the Bioanalyzer system (Agilent).

For DNA extraction, 450 μl of Buffer ATL from the QIAAMP DNA kit (QIAGEN) was added to the centrifuge tube, and DNA was extracted following the manufacturer’s instructions and quantified using QuBit (Thermo Fisher Scientific). Reverse Phase Protein Array (RPPA) arrays were carried out and analysed by MD Anderson Cancer Centre. Here Protein lysates were diluted and loaded onto nitrocellulose-coated slides that had been pre-conjugated with primary antibodies. Each protein was then visualised via a colorimetric reaction and quantified by Array-Pro Analyzer. The raw expression values were then normalised to protein loading and quantified by means of standard curves. Log₂ transformed and median-centred data were used for the downstream analyses. To perform molecular subtyping, total RNA from each individual case was reverse transcribed into cDNA, followed by amplification with a pool of indexed primers that target a predefined gene list (42 genes). The primers were selected from the Illumina Design Studios. The cleaned PCR product underwent QC by Tapestation (Agilent) to confirm the amplicon sizes. Forty-eight samples were multiplexed in one single MiSeq run. SR 50 bp were used to generate approximately 20 million reads per run.

Deep learning model development

We developed a U-Net-based segmentation model to facilitate end-to-end application of our models that require segmentation. We base our method on nn-UNet, a recently proposed fully automatic framework for the configuration of U-Net based segmentation methods, that has achieved state-of-the-art performance on 23 public datasets used in international biomedical segmentation competitions (15). We trained the nn-UNet stack of 2D, 3D, and 3D low resolution, and evaluated an ensemble of the four architectures (15).

As there are currently no CNN prognostic models to directly predict OS from CE-CT scans, we built several custom CNNs. ResNet-based architectures were applied without a predefined segmentation region-of-interest (16). ResNet-18, ResNet-32, and a simple CNN were adapted to predict OS by backpropagating the negative loglikelihood loss or the Weibull log loss, similar to that of the Cox proportional hazard and Weibull accelerated failure time model, in the final dense layer. We trained the models for 1,000 epochs with the AdaDelta optimizer and reduce-on-plateau call-back parameters. In the deep learning models to predict OS, we created a cropped around the region of interest, as standard, to prevent exploding or vanishing gradients often caused by non-heterogenous voxel intensity regions.

We developed our radiomics and CNN-based models on radiologist-checked adnexal lesion ground truths as a baseline.

Radiomics model development

With the 666 textural based radiomics features, computed using TexLab 2.0 (8), we evaluated a combination of 13 different feature reduction techniques with 12 machine learning algorithms, including linear and tree-based techniques, together with boosted and regularised variations. Linearly dependent features were removed, and corresponding lower ordered features were retained.

Several feature selection methods were employed: (I) correlation-based techniques: Pearson, Spearman, Kendall; (II) feature selection algorithms based on mutual information optimisation: joint information maximisation (JIM), mutual information feature selection (MIFS), minimum redundancy maximum relevance (MRMR); and (III) model-based approaches: a univariate Cox-regression model, a random forest variable hunting (RFVH), a random forest variable hunting with maximal depth (RF-MD), a random forest with variable hunting and variable importance (RFVH-VIMP), a random forest with variable hunting and Gini impurity corrected variable importance (IMPRF), and a random forest based on permutation variable importance (PVIRF) (17). Additionally, we also selected features at random for comparative purposes.

For each technique, radiomics feature selection was repeated 100 times by bootstrapping samples of the training cohort. Feature ranks were derived from each of the feature selection technique-specific by averaging the respective technique specific statistics. The number of radiomics features to include was treated as a hyperparameter based on this ranking.

The comparison of different machine learning algorithms included the following non-parametric models: (I) the Cox model, the Cox-Net, and Cox Lasso methods with lasso and elastic-net regularisations; (II) models based on boosting trees (BT): BT_Cox, BT_Cindex; (III) boosting gradient linear models (BGLM): BGLM_Cox, BGLM_Cindex; and (IV) random forest based methods: random survival forest (RSF), random forest using maximally selected rank statistics (MSR_RF), and a random forest with extra trees (ET_RF). Furthermore, we investigated the following full-parametric models (V): survival regression (Survival-Regression) and models based on the Weibull distribution: BT- and BGLM-Weibull.

Hyper-parameter optimisation was performed via grid-search with ten-fold cross-validation, optimising for Concordance index (C-index). Hyper-parameter ranges of the models are listed in Supplementary Table 1. A description of algorithms and feature reduction methods are detailed in the Supplementary Material. A bootstrapping (b=100) strategy was performed at the feature selection and model training stages (Figure 1), and then ensemble models were created. The best model feature importances were computed using the permutation variable importance technique.

Statistical analysis

To assess the model’s performance for the prediction of OS, we used C-index, hazard ratios (HR), the average time-dependent area-under-the-receiver operating characteristic (tAUC) curve over 1, 2, 3, 4 and 5 years, and Fisher’s exact test P-values. Kaplan-Meier analysis and log-rank tests were applied to assess survival outcomes among sub-groups with pathological complete response or incomplete response.

Quantitative statistics were presented as mean (SD) or median (IQR). Continuous variable distributions between datasets were compared using Kolmogorov-Smirnov (KS) tests, and categorical variables were compared using the χ² test. Spearman Correlation was performed to assess associations between the best model’s predicted probabilities and proteomics, copy number, and transcriptomics data.

The best performing medical image deep learning and radiomics models were selected based on C-index, as it is identical to the area under the curve (AUC) measure. K-means clustering was used on the training dataset derived predicted probabilities from our best performing CE-CT based radiomics model to define risk groups.

Gene Set Enrichment Analysis (GSEA)

In addition to predicting prognosis in HGSOC patients and performing Kaplan-Meier analysis, we performed GSEA analysis of genes associated with the high-risk group using the clusterProfiler R package. Using univariate logistic regression, we identified genes associated with the high-risk group; in addition we evaluated with the homo sapiens KEGG pathways (hsa), and Genome wide annotation for Human, primarily based on mapping using homo sapiens Entrez Gene identifier-based pathways (ord.HS.ef.db), that are provided by Bioconductor.

Clinical Characteristics of data

A total of 618 patients were included in this analysis Patient demographics and clinical parameters are a summarised in Table 2.

nn-UNet Segmentation models

The 2D, 3D-full resolution, and 3D-low resolution nn-UNet model took 81, 130, and 96 hours to train on a 24GB NVIDIA TITAN RTX. We assessed our models using dice score across the HH training, HH validation, TCIA test, and KEM test datasets. Boxplots of the dice scores are shown in Figure 3A. The median scores were 0.96, 0.90, 0.88, and 0.80 for the HH training, HH validation, KEM test, and TCIA test sets, respectively. In Figure 3B we show an example, whereby (in red) the nn-UNet algorithm segmented a bi-lateral lesion, together with the original radiologist segmentation coloured blue. Performance of the 3D-low resolution and 2D are shown in Appendix Figure 1; examples of a high and low performing 3D full-resolution nn-UNet segmentations (dice= 0.89, 0.55) are shown in Appendix Figure 2.

• Download figure
• Open in new tab

Figure 3. Illustration of 3D high -resolution nn-UNet performance.

(A) Boxplots of dice scores for all samples in Hammersmith Hospital (HH) (training (Train) and validation (Val)), Kliniken Essen-Mitte (KEM) and The Cancer Imaging Atlas (TCIA) test datasets. The median scores were 0.96, 0.90, 0.88, and 0.81 for the HH training, HH validation, KEM test, and TCIA test sets, respectively. (B) We show an example: in red is a predicted segmentation and in purple is the original segmentation.

Preliminary Analysis of Radiomics Data Structure

Data segmented as described above were further used in radiomics analyses. To further understand the radiomic characteristics of the HGSOC subtype, we performed unsupervised Spearman correlation analysis using the radiomic profiles in the HH training cohort (Figure 4A). We show that several radiomics features were highly correlated, and that the CT scanner manufacturer is not a main source of variance (Figure 4B). In the radiomics cohort from HH training cohort we extracted 666 radiomics features from 147 patients performed hierarchal clustering and defined 4 clusters (Figure 4C). Univariate Cox proportional hazard regression shows several radiomics features were significantly associated with OS (Figure 4D), indicating supervised machine learning models would yield value.

• Download figure
• Open in new tab

Figure 4 Radiomics Data Structure.

(A) Spearman correlation heatmap between radiomics profiles of the training dataset. (B) Principal components 1 and 2 derived from HH training data radiomics highlight the different scanner vendors. (C) Sample-wise radiomics profiles with hierarchal clustering based dendrograms annotated patient-wise (left) and radiomics feature wise (top), radiomics feature type, FIGO stage, CA-125, primary chemotherapy outcome and residual disease annotated respectively. (D) Volcano plot of univariate Cox proportional hazard model log₂ transformed P-values against the corresponding hazard ratios of the training dataset.

Supervised Models

In Figures 5A-D, we display the radiomics results for HH training, HH validation, TCIA test, and KEM test sets; for all models used. For OS, the best performing radiomics bootstrap-ensemble model based in terms of external-validation was the Random Survival Forest, built on features selected from the Random Forest-based Feature Selection, using permutation technique described by Janitza and colleagues (18) to calculate the variable importance, i.e. two random forests were implemented, one for feature selection, and a second for prediction on the selected features. We termed this combination – Permutation-Variable Importance Random Forest – Random Survival Forest (PVIRF-RSF). This radiomics model consistently demonstrated high C-index (HH train: 0.85 ± 0.01, HH validation: 0.66 ± 0.06, TCIA: 0.72 72 ± 0.05, KEM: 0.60 ± 0.01). In addition, by dichotomizing the predicted risk (probabilities) from the PVIRF-RSF model using k-means clustering, (threshold, 0.425), we show that this model is capable of stratifying individuals into significantly distinct high and low risk sub-groups, with log-rank P-values: <0.000001, 0.0057, 0.0044, 0.0055, for HH training, HH validation, TCIA and KEM datasets, respectively.

Whilst k-means clustering showed two distinct high and low risk groups, the number of patients in each group did not capture the most at risk patients. Aiming to derive individuals with highly aggressive phenotypes we also stratified the predictions of PVIRF-RSF model in to three distinct groups termed high, intermediate and low using k-means (k=3 clusters) derived thresholds of 0.33 and 0.44. These risk groups were significantly distinct (P-values: <0.000001, 0.039, 0.014, 0.035) for HH train, HH validation, TCIA, and KEM, respectively) and could capture the 5-30% of patients with increased risk.

When we implement the pipeline (including resampling) for the previous modelling technique that derived the radiomics prognostic vector (Table 5) described in Lu and colleagues (2019) within the cells corresponding to “Cox” on the x-axis and “Cox_Lasso” on the y-axis (Figure 5A-D) we obtain low C-indexes (0.56, 0.52, 0.54, 0.51 for HH train, HH validation, TCIA, and KEM, respectively).

We show univariate Cox Proportional Hazzard performance (C-index, tAUC, and HR) for conventional indicators of HGSOC prognosis, namely age, CA-125, FIGO stage, residual disease CCNE1, BRCA1 and BRCA2 (copy number variation) for the HH train, HH validation, TCIA and KEM test datasets (Appendix Table 3). The deep learning-based (Simple-CNN, ResNet18, and ResNet34) model performance (C-index, tAUC, and HR) are presented in Appendix Table 4. In Table 4 we compare the PVIRF-RSF model with several traditional clinical models and with the custom-built deep learning-based models (Simple-CNN, ResNet18, and ResNet34).

Biological interpretation of the radiomics Model

To understand tumour biological characteristics linked to the PVIRF-RSF predictions (Figure 6A), we performed a gene-set enrichment analysis on the corresponding TCIA and HH mRNA data. A volcano plot of the TCIA data with significantly enriched genes highlighted, is displayed in Figure 6B. Here, we found 10 genes to be significantly associated with the high-risk group. We found that Amoebiasis, alpha-Linolenic acid metabolism, and MAPK signalling pathway were the three pathways most significantly enriched for associations with high PVIRF-RF risk groups (P-value= 0.003, 0.003, 0.006). Since several amoebic proteins including lectins, cysteine proteineases, and amoebapores are associated with the invasion process, the results suggest emphasis of an invasive phenotype in the high-risk group. Emphasis of SNARE interactions in vesicular transport, basement membrane, bicellular tight junction and response to fibroblast growth factor (Figure 6 C, D) further support this notion. We show there is a significant increase in Tumour cellularity in the high-risk group the low-risk group (Figure 6 E). In addition, there is a small but insignificant decrease in Rad51 protein expression in the high-risk group compared to the low-risk group within the HH cohort (Figure 6 F), although there is a small significant decrease in Rad51 protein expression in the high-risk group compared to the low-risk group in the TCIA data set (Appendix Figure In the TCIA proteomics data, the features that were negatively correlated JNK2, YB1, RAD51, MTORPS448, and HER3PY1298 were the most correlated with the predictions. The most significantly associated protein, STAT5ALPHA was positively correlated with the predicted probabilities (Figure 7 and supplementary Table 3).

• Download figure
• Open in new tab

Figure 6 Volcano Plot and GSEA of Permutation-Variable Importance Random Forest – Random Survival Forest (PVIRF-RSF) Model.

(A) Radiomics features and there premutation computed importance in the PVIRF-RSF model (B) Volcano plot of genes expressed that are associated with the predicted probabilities of PVIRF-RSF radiomics model (TCIA data). (C) Top 10 KEGG based activated and suppressed pathways probabilities of PVIRF-RSF (TCIA data). (D) Top 10 Genome wide annotation for Human, primarily based on mapping using Entrez Gene identifiers (ord.HS.ef.db) activated and suppressed pathways probabilities of PVIRF-RSF (TCIA data). (E) HH data Tumour cellularity box plot across risk groups. (F) HH data Rad51protein expression across PVI-RSF risk groups.

• Download figure
• Open in new tab

Figure 7 Univariate Analysis Testing for key variables associated with predictions of the Permutation-Variable Importance Random Forest – Random Survival Forest (PVIRF-RSF) Model Derived from TCGA data

(A) The predicted probabilities for all 71 subjects (x-axis) in the TCIA cohort based on inference from the PVIRF-RSF Model and Clinical Variables. (B) Top ten RPPA variables most correlated with predictions of PVIRF-RSF (x-axis corresponds to each sample). (C) Top ten RNA variables most correlated with predictions of PVIRF-RSF (x-axis corresponds to each sample). (D) Top ten Copy Number variables most correlated with predictions of PVIRF-RSF (x-axis corresponds to each sample).

In the TCIA copy number variation data, the genes LRRC7, RN7SL392P, GNG12.x, DIRAS3, WLS, MIR1262, RPE65, LEPROT, snoU13|ENSG00000238931.1, RALGPS1 were the most correlated (negative correlation) (Figure 7 and Appendix Table 4). In the TCIA RNA expression data, the genes TF, FNDC4, ZNF671, PC, FAAH, MIZF, SOX30, PTGES2, B3GNT4, FPGS were the most correlated (negative correlation) (Figure 7 and Appendix Table 4).

In https://data.mendeley.com/xxx we provide the corresponding coefficients, Spearman Rho coefficients, and adjusted P-values. In the HH cohort we display a significant increase in tumour cellularity in the high-risk group compared to the low-risk group (Supplementary Figure 4).

The study by Lu et al. (2019) suggested that RPV emphasised an activated stromal phenotype. From the foregoing, demonstrate that the new model, PVIRF-RSF, emphasised different biology, we compared examined its relationship to histology data. PVIRF-RS is strongly positively correlated with tumour cellularity and not fibronectin indicating marked differences between the two models, with the PVIRF-RS emphasising the tumour (invasive) rather than the stromal compartment.