Gene expression profiles from individual HGS-OvCa tumor samples exhibit features of multiple subtype signatures.
It was previously reported that a portion of ovarian tumor samples are characterized by high numbers of infiltrating T lymphocytes or stromal cells and that such tumor samples can be recognized by expression profiling (23
). Using immunohistochemistry, we identified an increased number of samples with infiltrating T lymphocytes in the Immunoreactive group (42% [n
= 19] vs. 16% [n
= 80]; P
< 0.01, Fisher exact test), whereas desmoplasia, associated with infiltrating stromal cells, was found more often in the Mesenchymal cluster, with borderline significance (27% [n
= 26] vs. 12% [n
= 63]; P
= 0.07, Fisher exact test). To investigate the possibility that tumor samples expressing the Immunoreactive and Mesenchymal signatures were additionally expressing the non–infiltrating cell–associated Differentiated or Proliferative signature, we used single-sample gene set enrichment analysis (ssGSEA) to assess gene set activation scores for all 489 tumor samples. ssGSEA is a rank-based comparison of the expression levels of genes in the gene set with all other genes in an expression profile. In this analysis, expression profiles clustering with, for instance, the Differentiated group will have a high Differentiated signature score. Signatures of the 4 subtypes were described previously (20
). To define ssGSEA score cutoffs for inclusion in the set of samples expressing a subtype signature, we first assessed the silhouette width for each clustered profile. The silhouette width is defined as the ratio between the smallest correlation between a profile and all other profiles within its cluster and the largest correlation between that profile and all other profiles outside its cluster. The minimum ssGSEA score among samples that defined the original clusters and that clustered with a positive silhouette width to establish the ssGSEA score was used as cutoff value for inclusion into a subtype set. In this analysis, 82% of the 489 expression profiles were assigned to at least 2 subtypes (Figure and Supplemental Table 1; supplemental material available online with this article; doi:
). For comparison, a similar analysis on our previously reported subtypes of glioblastoma (25
) assigned 24% of expression profiles to more than 1 class (Supplemental Figure 1). Every possible signature combination was observed in at least 1 tumor sample, and negative correlations were observed between the Immunoreactive and Proliferative signatures (r2
= –0.81) and between the Differentiated and Mesenchymal signatures (r2
= –0.57). The overlap in gene signature scores suggests that HGS-OvCa does not consist of mutually exclusive expression subtypes, but that each tumor sample is represented by multiple signatures at different levels of activation. This pattern may reflect a higher level of homogeneity than is seen in other tumor types, such as glioblastoma and breast cancer.
Figure 1 HGS-OvCa signature ssGSEA scores. ssGSEA scores for 489 HGS-OvCa tumor samples were generated using 4 previously described gene expression signatures (20): differentiated (D), immunoreactive (I), mesenchymal (M), and proliferative (P).
Characterization of gene expression signatures.
To assess whether the 4 HGS-OvCa subtype expression signatures were associated with specific genomic aberrations, we investigated whether gene set activation scores were significantly increased in altered versus nonaltered samples. Genes with alterations in more than 3% of samples were included in the analysis (copy number alterations, n
= 489; mutation events, n
= 316). After Benjamini-Hochberg correction of gene set score t
values for multiple testing, no gene mutation or copy number deletion resulted in a significant change of gene set score. Copy number amplifications of 367 genes, mapping to 13 different commonly amplified regions, were found to be associated with a significant increase of activation levels of at least 1 of the 4 ovarian cancer signatures (Supplemental Table 2). To gain further insight into the possible target of amplification for the 13 amplification regions, we calculated Pearson correlation coefficients among 279 amplified genes with matching expression data and gene set activation scores. Of 11 genes with matching expression data at 1q32.1, LAD1
, which regulates stability of epithelial layer stability with mesenchyme, showed a correlation with the Differentiated gene set that was distinctly higher than other genes. When linking 105 genes at the 19q13 locus with expression levels to the Proliferative signature, we found NOTCH3
to be among the most highly correlated. NOTCH3 inactivation in samples expressing the Proliferative signature may provide a therapeutic opportunity (26
). The 12q14 gene whose expression correlated most highly with the Proliferative gene set was HMGA2
, which suggests a role for this chromatin-modifying gene that has previously been related to progression in lung adenocarcinoma (27
). Amplification of the transcription factor MYC
was associated with a significant difference in Differentiated and Immunoreactive gene set scores, but correlation was much higher with Differentiated than with Immunoreactive signature scores.
The gene signatures were further characterized by correlating signature activation scores to the activation scores of 7,565 gene sets from the Molecular Signatures Database (
; ref. 28
) and gene expression levels (see Supplemental Methods).
Predicting outcome using an integrated signature classifier.
We previously described a transcriptional signature of 193 genes predictive of overall survival, which was developed using univariate Cox regression analysis on the expression profiles of 215 HGS-OvCa samples (20
). In order to maximize the discovery power of the TCGA dataset, we used all 481 TCGA HGS-OvCa–derived expression profiles with matching clinical data, rather than the 215 profiles used in our previous work (20
). To generate the new signature, we selected the 100 genes whose expression was most correlated (good prognosis genes, n
= 47), or anticorrelated (poor prognosis genes, n
= 53), with survival (Supplemental Table 3). We chose a signature of 100 genes to prevent sensitivity to overfitting of the data and to permit analysis of new samples with a medium-throughput expression assay. We named this signature the “Cl
assification of Ovar
ian Cancer” (CLOVAR) survival signature. Gene set enrichment analysis of the 100-gene CLOVAR survival signature, which includes genes such as RB1
, and RXRB
, showed an overlap between the signature and Molecular Signatures Database gene sets such as BENPORATH_MYC_MAX_TARGETS and LY_AGING_OLD_UP (Supplemental Table 4). Of the 100 genes in the CLOVAR survival signature, 36 were included from our recently reported 193-gene survival signature (20
To verify the prognostic power of the CLOVAR survival signature, we used ssGSEA to analyze 64 TCGA expression profiles not included in the training dataset and 815 expression profiles with matching survival annotation from 6 published studies (9
). CLOVAR survival signature ssGSEA scores were calculated for each of the 7 validation datasets. The validation cohort included data generated on 5 different expression platforms, and the scale of the ssGSEA scores was dependent on the total number of genes included in the expression platform used. To make ssGSEA scores comparable between different expression platforms, scores were converted to a [0,1] scale within each dataset (see Methods). The validation set was stratified into CLOVAR poor prognosis or CLOVAR good prognosis groups using a cutoff score of 0.5. The difference in survival between the 2 prognostic groups was highly statistically significant (P
< 0.0001; Figure A).
Survival per CLOVAR subtype, CLOVAR survival class, and double classifiers.
A significant difference in recurrence-free survival was observed between patients in the TCGA cohort whose tumor samples expressed either the Immunoreactive or the Mesenchymal signature (Supplemental Figure 2). Next, we set out to verify this survival characteristic in the validation cohort. As reducing the size of the signature allows for the use of medium-high throughput gene expression quantification systems, we first reduced the 4 subtype signatures from 800 genes to gene signatures that, combined, contained 100 genes. We refer to these herein as CLOVAR subtype signatures. The CLOVAR subtype signatures classified the TCGA cohort into the 4 groups used to define the signatures with a cross-validation error rate of 8%, which indicates that the original clustering was captured by the reduced signatures. In order to determine whether validation cohort samples expressed the subtype signatures, we assessed ssGSEA scores of CLOVAR subtype signatures in the TCGA expression dataset. The minimum ssGSEA score among core samples was used to establish the ssGSEA score used as a cutoff value for inclusion into a subtype set. Using these ssGSEA cutoff values, tumor samples in the validation cohort were identified as expressing the CLOVAR Immunoreactive or CLOVAR Mesenchymal signatures. Cases expressing both signatures were assigned to a class according to the higher of the 2 scores. Those samples that did not express one of the signatures related to infiltrating cells were assigned to the CLOVAR Proliferative, CLOVAR Differentiated, CLOVAR Immunoreactive, or CLOVAR Mesenchymal class depending on the highest normalized gene set score for the 4 CLOVAR gene set signatures. Kaplan-Meier analysis revealed a highly significant difference in survival among the 4 validation set expression groups established by this approach (P < 0.0001; Figure B).
To combine the CLOVAR survival and CLOVAR subtype signatures for outcome prediction, we established a multiple covariate model using the TCGA dataset and including 3 signature scores as variables: CLOVAR survival, CLOVAR Immunoreactive, and CLOVAR Mesenchymal. In this model, CLOVAR survival was the dominant variable (Supplemental Table 5). Based on this model with CLOVAR survival and CLOVAR subtype as covariates, we computed risk scores of patients in training and validation sets. Patients in each set were then classified into good, intermediate, and poor outcome groups according the risk scores, using tertiles of risk scores for the TCGA set as cutoffs. We used this model to compute risk scores and classified the validation set (n = 879) into good, intermediate, and poor outcome groups, which showed 5-year survival of 55%, 38%, and 29%, respectively. Next, we tested whether including additional prognostic factors could improve the predictions. CLOVAR survival, CLOVAR Immunoreactive, and CLOVAR Mesenchymal scores, stage, grade, residual disease, and age annotation were available on 428 TCGA patients, which were used to calibrate a model (Supplemental Table 5). Using this extended model, we classified 328 patients from the validation set for which these parameters were available. The 5-year survival of the good, intermediate, and poor outcome groups using the extended model was 60%, 39%, and 27%, respectively. Lastly, we tested whether addition of “presence of BRCA1/BRCA2 germline mutations” to the model could further inform prognostic classification. BRCA1/BRCA2 germline data were available for 273 patients in the TCGA set and 83 patients in the validation set. The 5-year survival of the good, intermediate, and poor outcome groups in the validation set using the extended model including BRCA1/BRCA2 germline mutations was 61%, 25%, and 0%, respectively. To compare the performance of the 3 models, we computed likelihoods for the TCGA and validation sets using high/intermediate/low classification as a covariate, only including samples with nonmissing data on all variables (TCGA, n = 273; validation set, n = 83; Figure ). Both in the TCGA set, which was used to train the models, and in the validation set, the extended model combined with BRCA1/BRCA2 germline mutation status showed the highest likelihood, suggesting most accurate performance in outcome prediction. Although the comparison may be hampered by the specific set of cases for which the correct annotation was available, this result suggests that a model combining CLOVAR survival, CLOVAR Immunoreactive, and CLOVAR Mesenchymal scores and BRCA1/BRCA2 germline mutation status provides optimal outcome predictions.
Survival analysis of multiple covariate outcome predictions.
germline mutation status was only available for 9.4% of samples in the validation set, we sought to further optimize outcome prediction through combination of the CLOVAR survival and CLOVAR subtype signatures. We performed Kaplan-Meier survival analysis between samples with CLOVAR poor prognosis and CLOVAR good prognosis survival signatures within each of the CLOVAR Immunoreactive, CLOVAR Mesenchymal, CLOVAR Proliferative, and CLOVAR Differentiated subtypes. Importantly, the CLOVAR survival classifier predicted outcome with high significance within each of the 4 CLOVAR subtype groups (Figure C). The worst outcome was found in patients with tumors classified as both CLOVAR Mesenchymal and CLOVAR poor prognosis, which represented 23% of the validation set and showed a median overall survival of 23 months. The hazard ratio for death within 5 years of diagnosis for this group compared with the remainder of the sample cohort was 1.95 (95% confidence interval, 1.59 to 2.40), and the hazard ratio for recurrence within 12 months was 1.66 (95% confidence interval, 1.26 to 2.18) (Supplemental Table 6). Of the CLOVAR Mesenchymal/CLOVAR poor prognosis cases, 63% were determined to be resistant to platinum therapy, defined as having tumor recurrence within 6 months after the end of platinum therapy (Supplemental Table 1 and Supplemental Figure 3). All other samples — those not classified as CLOVAR Mesenchymal/CLOVAR poor prognosis — showed a platinum resistance rate of 23%. To refine the classification, we sought to understand those CLOVAR Mesenchymal/CLOVAR poor prognosis patients that had a relatively favorable outcome, with greater than 12 months progression-free survival (PFS). Germline BRCA1
mutations are associated with increased response rates to chemotherapy and overall survival (32
). However, germline BRCA1
mutations were not found more frequently among tumor samples predicted as CLOVAR Mesenchymal/CLOVAR poor prognosis with greater than 12 months PFS compared with those CLOVAR Mesenchymal/CLOVAR poor prognosis samples with less than 12 months PFS. Similarly, such patients were not more likely to have been optimally debulked, were not younger, and did not receive a different treatment regimen compared with patients for which there was concordance between molecular prediction and poor clinical outcome. Interestingly, the amount of residual tumor after surgical debulking was independently related to prognosis in the CLOVAR Mesenchymal/CLOVAR poor prognosis, CLOVAR Immunoreactive/CLOVAR poor prognosis, and CLOVAR Differentiated/CLOVAR poor prognosis groups (Figure D and Supplemental Figure 4). Understanding the misclassification of these patients requires further investigation and will assist in more accurate prediction of PFS.