Pooled Clustering of High-Grade Serous Ovarian Cancer Gene Expression Leads to Novel Consensus Subtypes Associated with Survival and Surgical Outcomes

Ovarian cancer is the second most common and the most lethal gynecologic malignancy in the United States (1). Despite radical surgery and initial high response rates to platinum- and taxane-based chemotherapy (2), most patients experience a relapse, with an estimated 14,240 deaths in 2016 in the United States alone (3). Further understanding of the biology of ovarian cancer is required to improve clinical outcome. Microarray-based gene expression profiling has been used to classify high-grade serous ovarian carcinomas (HGSOC) into molecular subtypes but this has had limited clinical utility, in part due to lack of modifiable treatment options. In addition, there have been multiple related but different subtype definitions proposed, each derived from individual studies with potential platform-specific biases (4–6), lack of detailed clinical data, and wide variation in surgical resection rates confounding survival outcomes. We sought to investigate whether HGSOC molecular subtypes associated with surgical outcomes as well as with survival. To improve power to discern molecular subtypes, we first pooled expression analysis of common genes across 14 expression studies with adjustment for batch-effects; we then established de novo molecular subtypes and examined associations with survival and surgical outcomes. In a smaller collection of Mayo Clinic patients with more detailed baseline surgical and updated clinical follow-up data, we further assessed association with molecular subtype and clinical outcomes.

Fourteen public expression datasets were retrieved from curatedOvarianData (7), an ovarian cancer transcriptome database; we excluded studies with less than 50 tumors due to minimum sample size for batch correction. We further restricted to HGSOC cases (grade > 1, serous; refs. 8, 9) and excluded chemo-treated tumors and duplicates with the Mayo Clinic collection, resulting in 2,103 cases (Supplementary Table S1). A total of 1,017 deaths occurred, and among living cases mean follow-up was 42 months (range, 0–243). The extent of postsurgery residual disease was available as optimal debulking (≤1 cm, analogous to Mayo's RD0 and RD1 combined) and suboptimal debulking (>1 cm, analogous to Mayo's RD2).

Eligible cases (n = 381) were ascertained at the Mayo Clinic between 1992 and 2009. A gynecologic pathologist (G. Keeney) confirmed initial diagnosis of HGSOC and reviewed tissues to ensure ≥70% tumor content. Survival data were obtained from the Mayo Clinic Tumor Registry, electronic medical records, and active patient contact; 317 deaths occurred as of November 3, 2015. Among living cases mean follow-up was 101 months (range 30–202). The extent of postsurgery residual disease was available as RD0 (no macroscopic disease), RD1 (measureable but ≤1 cm), and RD2 (>1 cm). Clinical characteristics along with associations with overall survival (HRs; 95% confidence intervals, CI) are shown as Supplementary Table S2. All cases provided informed consent for use of their tissues and data in research; all protocols were approved by the Mayo Clinic Institutional Review Board.

The overall analytic flow is shown as Supplementary Fig. S1.

The probe-set with the highest mean across all datasets of a platform was utilized when multiple probe-sets mapped a gene (10), resulting in 7,361 common genes across studies. Initial study-batch effects were corrected by ComBat (11). After adjustments, the residual batch differences by study and microarray platform were assessed by association analysis of the first 10 principal component (PC) projections, revealing the third and tenth PCs still significantly associated with studies and platform. The residuals from a linear model of the expression value for each gene and both PC3 and PC10 were computed. The expression mean value for each gene was added back to the residuals to obtain an adjusted expression value for each gene. After applying the linear regression–based method to eliminate per-study and/or per-platform batch effects across 14 datasets, we checked batch-effects according to PC analysis projections versus studies, found no residual batch differences associated with studies, shown as Supplementary Fig. S2.

RNA extracted from fresh frozen tumors was assessed using Agilent Whole Human Genome 4 × 44 K Expression Arrays as previously described (4, 12). Batch-effects among two different RNA preparations and three microarray profiling dates were corrected to adjust Cy5, Cy3 labeling differences observed among experimental batches (4, 12).

BRCA1/2 germline and somatic mutations, and BRCA1 promoter methylation, collectively called as “BRCAness,” were retrieved from supplementary tables of a previous study summarizing whole-exome sequencing results from TCGA ovarian cohort (13).

Four genomic scores for ovarian TCGA cases were retrieved from a previous pan-Cancer study (14): (i) loss-of-heterozygosity (LOH): a genomic scarring score previously studied in HGSOC and defined by Abkevich and colleagues (15); (ii) large-scale transition (LST): a type of genomic scar associated with BRCA1 or BRCA2 status in triple-negative breast cancer (16); (iii) telomeric allelic Imbalance (TAI): another type of genomic scarring score which was found predictive of platinum-based chemotherapy response in both breast and ovarian cancers (17); (iv) tumor mutation count: the total number of somatic mutations per tumor sample.

To define molecular subtypes using genes with substantial expression variability, we chose 872 genes that expressed at least 1.5 log-scale fold-change differences across at least 200 samples, compared with its own median expression of 2,103 HGOSC cases in the public cohort. De novo clustering was evaluated using a consensus non-negative matrix factorization (NMF) approach (18). From cluster numbers k = 2 to k = 9, 50 independent NMF runs were performed to evaluate cluster stability. We found that k = 2 to k = 5 were comparably stable (Supplementary Fig. S3) and selected k = 5 as the stability of clustering solutions began to drop significantly starting from k = 6. To select genes most predictive of subtypes, multinomial LASSO regression was adopted (19), and an optimal parameter was selected according to the so-called “one standard error” rule after 10-fold cross-validation (20). To achieve a parsimonious set of genes for future studies, the top 100 genes with largest average absolute LASSO coefficients were chosen as the molecular subtype signature. A de novo expression centroid for each molecular subtype was computed according to average expression of samples in each molecular subtype (Supplementary Table S3). For subtype assignment in public and Mayo Clinic data, Pearson correlation coefficients were computed for each subtype centroid and given sample, and the most highly correlated molecular subtype was assigned. No assignment was made if the highest coefficient was below 0.1.

For comparison, we also assigned molecular subtypes according to the TCGA (21) and a study previously done by Tothill and colleagues (6). TCGA subtype assignment in public and Mayo data was done according to gene signatures and similar procedures as defined previously (21). Briefly, gene set activation scores for each subtype signature were computed using single-sample gene set enrichment analysis (ssGSEA; ref. 22), and raw enrichment scores were expressed as relative z-scores to account for various sizes of signatures. Given one sample, ssGSEA assigned it to one of the four TCGA molecular subtypes with the highest enrichment z-scores. For the Tothill approach, data from the samples used for subtype assignments were retrieved from Tothill and colleagues (6), and genes predictive of Tothill's molecular subtypes were extracted from a related publication (23). Expression centroid per molecular subtype (C1/C2/C4/C5) was computed according to average expression of samples in a same subtype, with detailed expression centroids comprising of 841 genes shown as Supplementary Table S4. For a given sample, Pearson correlation coefficients were computed against subtype centroids, and the most highly correlated subtype in public and Mayo Clinic data was assigned. No assignment was made if the highest coefficient was below 0.1.

For each gene, differential expression for each de novo molecular subtype compared to other molecular subtypes combined was evaluated using t tests. Pathway enrichment analysis was done using Bioconductor package “GAGE” (24). Multiple hypothesis-testing corrected q-values were computed for both analyses (25).

Cox regression was used to estimate HRs and 95% CIs, including multivariate step-wise variable selection. Univariate analysis of clinical features revealed that age, grade, stage and surgical debulking status significantly associated with survival with P value <0.05. Following step-wise variable selection for the multivariate analysis, only stage and surgical debulking status remained significantly associated with overall survival and therefore these two factors were adjusted for multivariate associations with overall survival.

For each type of genomic scores (LOH, LST, TAI, and tumor mutation count), the two-sided Wilcoxon rank sum tests were performed to test the score differences comparing s5.ANM subtype versus s4.DIF and all the other subtypes, respectively.

Five stable de novo molecular subtypes (s1, s2, s3, s4, s5) of HGSOC were defined (Supplementary Table S3 and Supplementary Figs. S1–S3) using NMF-based consensus clustering of pooled data from 2,103 publicly available HGSOC expression profiles corrected for batch-effects. We developed a parsimonious 100-gene signature using LASSO-based analyses that discriminated these molecular subtypes, and assigned Mayo Clinic HGSOC tumors to the most likely subtype according to expression similarities as described in Materials and Methods. Independently in the public and Mayo Clinic datasets, we compared the assignments of de novo molecular subtypes with those of previous molecular subtyping systems defined by TCGA and Tothill studies (refs. 5, 6; Fig. 1). Overall, the de novo subtyping system was significantly associated with prior subtyping systems (χ² test P < 2.2e−16; Supplementary Table S5). Three de novo molecular subtypes were largely consistent with previous molecular subtype definitions (s1.MES: “mesenchymal”; s2.IMM: “immunoreactive”; s3.PRO: “proliferative”). However, we report two de novo molecular subtypes emerging from “TCGA-differentiated/Tothill-C4” subtype. We examined differential expression patterns between each subtype versus others, and found that the top differentially expressed genes of the novel s4 subtype was more consistent with the previously defined signatures of Tothill's C3/TCGA-differentiated subtype (Supplementary Tables S6–S10). Therefore, we refer to this as the s4.DIF (differentiated) subtype. The fifth novel subtype was named the “anti-mesenchymal” (s5.ANM) subtype because of its downregulation of genes that were typically upregulated in the mesenchymal subtype.

Classifications based on this five-level molecular subtyping system was significantly associated with overall survival of advanced stage patients in both public cohort (stage III or IV; P < 0.001; Supplementary Fig. S4), and Mayo Clinic cohort (stage IIIC or IV; P = 0.001; Fig. 2A). In univariate analysis of Mayo Clinic data, all three subtyping systems (five-level, Tothill, TCGA) showed significant associations with survival. For five-level subtyping system, the new s5.ANM subtype was observed with noticeably improved survival compared with s1.MES (Fig. 2A). In multivariate analysis adjusting for age, stage, and surgical outcome, the five-level molecular subtypes had the most significant associations with overall survival (P = 0.0004; Table 1). s2.IMM and s5.ANM tumors were significantly associated with longer survival time and lower death-risk compared with s1.MES (s2.IMM: HR 0.69; 95% CI, 0.49–0.97; P = 0.03; s5.ANM: HR 0.43; 95% CI, 0.28–0.64; P = 4.6E−5). In contrast, the Tothill subtyping system did not have statistical significance in multivariate analysis after adjusting for other clinical factors; for TCGA subtyping system, only the immunoreactive subtype was significantly associated with improved outcome compared with the mesenchymal subtype.

Among late-stage HGSOC patients with an optimal surgical debulking outcome (residual disease ≤ 1 cm), de novo molecular subtypes were significantly associated with overall survival in both the public (n = 936, P < 0.001) and Mayo Clinic cohorts (n = 250, P = 0.014; Supplementary Fig. S5). Among suboptimally debulked late-stage patients, these molecular subtypes were still significantly associated with survival in the public cohort (P = 0.005) but not in the Mayo Clinic cohort (P = 0.167); however, the de novo s5.ANM subtype showed a nonsignificant trend toward increased survival compared with other molecular subtypes especially in Mayo Clinic cohort (median survival 25.8 months vs. 19.6 months; Supplementary Fig. S6). A potential advantage of Mayo Clinic cohort is the availability of complete cytoreduction data to investigate the relationship between molecular subtypes and residual disease on survival. Noticeably, molecular subtype was not associated with overall survival for patients who had no residual macroscopic disease (RD0; P = 0.23; Fig. 2B) or high volume (RD2) surgical outcome (P = 0.167; Fig. 2D). We also evaluated debulking-survival associations within each molecular subtype, to determine if specific subtypes benefit more or less from surgical resection. We found unequivocal survival benefits of having RD0 surgical outcome regardless of de novo subtype membership (Supplementary Fig. S7). Interestingly, s2.IMM tumors with RD1 and RD2 surgical outcomes resulted in similar survival (Supplementary Fig. S7C), as did s4.DIF and s5.ANM tumors with RD0 and RD1 surgical outcomes (Supplementary Fig. S7E and S7F). This implies surgical goals may not be the same for all subtypes and that molecular subtype may be tested in combination with neoadjuvant treatment.

Further analysis examined whether HGSOC molecular subtypes correlated with surgical debulking outcomes scored on a trinary scale in the Mayo Clinic cohort. We found that s1.MES subtype was associated with the lowest RD0 rates across different subtyping systems, ranging from 12% to 16%, whereas the immunoreactive subtype consistently had the highest RD0 debulking rate, ranging from 40% to 46% (Fig. 3A; Supplementary Table S11). For 339 stage III-C or IV HGSOC patients in TCGA cohort with trinary scale of residual disease data, we again confirmed lowest RD0 rate in s1.MES subtype (6%) versus other subtypes ranging from 19% to 25% (Fig. 3B). In the public cohort using binary surgical outcomes (optimal vs. suboptimal debulking), surgical outcomes were also found to be significantly associated with HGSOC molecular subtypes (P < 2.2E−16): the s1.MES/C1/mesenchymal subtype was again associated with lowest optimal-debulking rates ranging from 49.1% to 51.6% across different subtyping systems.

Given our observation that tumors of the s1.MES molecular subtype showed less favorable surgical outcome, we hypothesized that the s1.MES subtype would also be associated with more complex surgical procedures than other molecular subtypes. To test this hypothesis, we enumerated surgical procedures in the Mayo Clinic cohort and computed surgical complexity scores (SCS; ref. 26) for advanced-staged HGSOC patients with optimal debulking outcomes (i.e., RD0 or RD1). We excluded suboptimal/RD2 cases, which often involve unresectable cases and aborted surgical attempts, resulting in paradoxical lower SCSs despite more advanced disease. We found that the s1.MES subtype tended to have higher complexity scores than other molecular subtypes in RD0 and RD1 cases in the Mayo cohort, supporting the hypothesis (Fig. 3B). Several recent reports have suggested specific gene signatures may be predictive of the ability to surgically resect HGSOC, potentially serving as biomarkers for surgical decisions (27–29). When expression levels of these genes were compared between the s1.MES subtype and others, 9 of 10 genes in these previously defined surgical resection signatures were significantly upregulated in the s1.MES subtype (Supplementary Table S12). This observation further supports the relationship between the s1.MES subtype and surgical debulking outcome.

To shed light on potential biological mechanisms underlying the relationship between molecular subtypes and outcomes, we performed differential expression analysis and pathway enrichment (Supplementary Figs. S8–S17; Supplementary Tables S7–S11, and S13–S17). In these analyses, the s1.MES subtype showed strong upregulation of cancer-associated fibroblast signatures (e.g., FAP) and EMT/stem-cell drivers (e.g., ZEB1 and SNAIL2), as well as significantly increased expression levels of fibronectin (FN1). The s2.IMM subtype showed significant upregulation of antigen presentation machinery (e.g., TAP1, PSMB8, and PSMB9); and its top subtype marker TAP1 was highly correlated with expression of PD-L1 (R_Spearman = 0.72; P < 2.2E−16; Supplementary Fig. S18), a marker has recently been associated with increasing tumor infiltrating T cells and favorable prognosis in HGSOC (30). Finally, our analysis involving unsupervised clustering showed that the TGFβ pathway is upregulated in the mesenchymal subtype (q-value = 1.5 × 10⁻⁴; Supplementary Table S13), which carries a worse prognosis (Fig. 2A).

To shed some lights on novel s5.ANM subtype, we performed additional analysis to compare this novel subtype versus other subtypes in TCGA cases for which genomic data are available. We first examined distribution of BRCAness, including BRCA1/2 mutation and BRCA1 promoter methylation status, among these five de novo subtypes, shown as Supplementary Table S18. We found a trend that s3.PRO and s5.ANM subtypes had lower mutation rate in BRCA1 and BRCA2 (15.3% and 16.7%) comparing with all the other three subtypes (25%–29.6%). After accounting for BRCA1 methylation status, s3.PRO and s5.ANM had higher wild-type BRCA/2 rates (76.3% and 75%, respectively), comparing with s1.MES (62.9%), s4.DIF (68.2%), and s2.IMM (57.4%). In TCGA cohort, we further evaluated associations between de novo subtypes and four tumor genomic scores summarized by a previous study (14). Among four genomic scores considered, 5.ANM subtype was associated with statistically higher telomeric allelic imbalance scores comparing to s4.DIF subtype (P = 2.83e−2). In addition, s5.ANM subtype was found associated with statistically lower tumor mutation count comparing with s4.DIF subtype (P = 2.97e−3), and comparing with all the other subtypes (P = 4.89e−5), shown as Supplementary Fig. S19.

This study utilized transcriptional profiling from over 2,100 cases and unsupervised clustering to derive molecular subtypes of HGSOC, then applied this revised classification scheme to the original cases as well as a carefully annotated Mayo Clinic cohort to assess the relationship between subtype, surgical resectability, and survival. To our knowledge, this study is the most comprehensive HGSOC molecular subtype analysis using a study- and platform-independent approach. This large-scale analysis confirmed the previously described mesenchymal (s1.MES), immunoreactive (s2.IMM), and proliferative (s3.PRO) subtypes as distinct entities. While our work was in progress, a cross-population analysis has also confirmed at least mesenchymal and proliferative subtypes were stable through clustering comparisons across individual expression studies (31). In this study, a new subtype “anti-mesenchymal” (s5.ANM) was derived from TCGA-differentiated/Tothill-C4 subtype, with noticeable downregulation of mesenchymal signature genes and better survival outcomes in both the public and Mayo Clinic HGSOC cohorts. To summarize our current analysis and linkage of the different subtyping systems, we made a schematic summary plot shown as Fig. 4.

In view of the prognostic importance of surgical debulking for HGSOC and the modifiable use of neoadjuvant treatment, we examined the relationship between subtype and extent of residual disease. In the Mayo Clinic cohort, we found that the mesenchymal subtype had the lowest RD0 rates across different HGSOC subtyping systems. Among optimally debulked cases, we also observed that the mesenchymal subtype tumors had higher surgical complexity scores than other molecular subtypes. It does appear that the mesenchymal subtype is associated with higher rates of unsuccessful and/or more complex surgical resections. This is also supported by the fact that 9 of 10 genes previously reported to be correlated with debulking outcome are significantly upregulated in mesenchymal subtype (27–29). Collectively these observations suggest a potential role for subtyping in surgical planning or referral to expert centers. Not surprisingly, the de novo molecular subtypes were significantly associated with survival in both the public and Mayo Clinic cohorts when surgical debulking status was classified in a binary manner (≤1 cm vs. >1 cm remaining). Theoretically, if tumors were completely removed, molecular subtypes should have no impact on survival; in fact, we observed that molecular subtypes were not associated with survival (P = 0.23) in Mayo Clinic cases with no macroscopic diseases (RD0 using the trinary system, n = 97). Because this observation was made in a small single-site population, further validation in additional cohorts is required to understand possible interactions between surgical efforts and tumor molecular subtypes in a clinical trial setting.

Although the molecular subtyping of HGSOC by transcriptional profile is appealing, how this subtyping can facilitate basic and clinical research into etiology, prognosis, or therapeutic potential is still unclear (32–34). For example, existing prognostic signatures or molecular subtypes have not been able to identify platinum-refractory HGSOC with sufficient accuracy to allow patients with a pertinent profile to be triaged to a clinical trial in place of standard chemotherapy (29). Nevertheless, HGSOC molecular subtyping could provide important biological insight and facilitate therapeutic studies. For example, this study identified several subtype-specific pathways that were potentially targetable, and therefore worthy of future study. In our analysis, TGFβ pathway was upregulated in tumors assigned to mesenchymal subtype, which carries a worse prognosis. Several other supervised expression studies have also reported that TGFβ pathway activities are associated with worse clinical outcomes and ovarian cancer metastasis (29, 35–37). Therefore, tumors with the mesenchymal gene expression pattern might be considered for future trials containing TGFβ inhibitors. With other developing therapeutic options such as PARP inhibitors and immune therapies for ovarian cancer (38, 39), it might also be important to consider treatment responses in the context of HGSOC molecular subtypes. Noticeably, laparoscopy assessment has been used in MD Anderson study to subjectively evaluate disease burden and resectability before cytoreduction surgery (40), making presurgery biopsy and molecular subtyping possible. In the future, knowledge of HGSOC molecular subtypes might help to identify patients for the most appropriate treatment options.

To utilize molecular subtyping clinically, a robust assay with reproducible subtype assignment is needed. There are examples of clinically implementing reliable and reproducible subtyping assays in other cancers, including breast cancer and lymphoma (41, 42). Noticeably, research evaluation of clinically feasible panels of ovarian subtypes is underway, pinpointing a 48-gene based NanoString panel as a promising platform also working for formalin-fixed paraffin-embedded samples (43). Once successfully implemented, molecular subtypes could be captured in prospective clinical trials to evaluate subtype-specific treatment benefit and serve as biomarkers for facilitating patient stratification. For example, we recently found PRO and MES molecular subtypes may derive a comparably greater benefit than DIF and IMM subtypes, in a retrospective study in a randomized controlled phase III trial of primary ovarian cancer with bevacizumab treatment (44).

In summary, unsupervised clustering of gene expression profiles from a large public HGSOC cohort (n = 2,103) identified five molecular subtypes of HGSOC: three of these molecular subtypes, proliferative, mesenchymal, and immunoreactive, were largely consistent across different subtyping systems (Tothill and TCGA), but two additional molecular subtypes were split from the previously described Tothill-C4/TCGA-differentiated subtype. Importantly, the novel anti-mesenchymal subtype, which was previously included in the differentiated subtype, is associated with improved survival in both the public and Mayo Clinic HGSOC cohort (n = 381). Further analysis demonstrated the potential prognostic significance of this subtyping analysis. In particular, the mesenchymal subtype was associated with a much lower RD0 rate and increased surgical complexity in Mayo Clinic cohort and a lower optimal debulking rate in the public cohort. Future work in other clinically annotated HGSOC cases is needed to evaluate whether molecular subtyping can be used to guide treatment decisions in HGSOC (e.g., neoadjuvant vs. adjuvant options).

No potential conflicts of interest were disclosed.

Conception and design: C. Wang, E.P. Heinzen, G.L. Keeney, W.A. Cliby, A.L. Oberg, E.L. Goode

Development of methodology: C. Wang, E.P. Heinzen, G.L. Keeney, W.A. Cliby

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K.R. Kalli, M.J. Maurer, E.P. Heinzen, G.L. Keeney, E.L. Goode

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C. Wang, S.M. Armasu, E.P. Heinzen, G.L. Keeney, A.L. Oberg, S.H. Kaufmann

Writing, review, and/or revision of the manuscript: C. Wang, S.M. Armasu, K.R. Kalli, M.J. Maurer, E.P. Heinzen, G.L. Keeney, W.A. Cliby, A.L. Oberg, S.H. Kaufmann, E.L. Goode

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): K.R. Kalli, E.P. Heinzen, G.L. Keeney

Study supervision: C. Wang, E.P. Heinzen, G.L. Keeney, S.H. Kaufmann

This work was supported by a career enhancement award from the Mayo Clinic Ovarian Cancer SPORE (P50 CA136393; to to C. Wang), National Cancer Institute funding (R01 CA122443; to E.L. Goode), and funds from the Mayo Clinic Cancer Center (P30 CA015083) and the Fred C. and Katherine B. Andersen Foundation to all authors.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.