Published in Volume
123, Issue 9
(September 3, 2013)J Clin Invest.
Copyright © 2013, American Society for Clinical Investigation
Predicting time to ovarian carcinoma recurrence using protein markers
1Department of Bioinformatics and Computational Biology, The University of Texas MD Anderson Cancer Center, Houston, Texas, USA.
2Department of Obstetrics and Gynecology, Niigata University Graduate School of Medical and Dental Sciences, Niigata, Japan.
3Niigata Medical Center Hospital, Niigata, Japan.
4Department of Obstetrics and Gynecology, Kagoshima City Hospital, Kagoshima, Japan.
5Department of Obstetrics and Gynecology, School of Medicine, Nagasaki University, Nagasaki, Japan.
6Department of Obstetrics and Gynecology, Tottori University School of Medicine, Tottori, Japan.
7The Cancer Genome Atlas, National Cancer Institute, NIH, Bethesda, Maryland, USA.
8Department of Obstetrics and Gynecology, National Defense Medical College, Tokorozawa, Japan.
9Department of Obstetrics and Gynecology, Kurume University School of Medicine, Kurume, Japan.
10Division of Gynecologic Oncology, University of Pennsylvania, Philadelphia, Pennsylvania, USA.
11Department of Systems Biology, The University of Texas MD Anderson Cancer Center, Houston, Texas, USA.
Address correspondence to: Roel G.W. Verhaak, Department of Bioinformatics and Computational Biology, The University of Texas MD Anderson Cancer Center, 1400 Pressler Dr., Unit 1410, Houston, Texas 77230-1402, USA. Phone: 713.563.2293; Fax: 713.563.4242; E-mail:
First published August 15, 2013
Received for publication December 26,
2012, and accepted in revised form June 6,
Patients with ovarian cancer are at high risk of tumor recurrence. Prediction of therapy outcome may provide therapeutic avenues to improve patient outcomes. Using reverse-phase protein arrays, we generated ovarian carcinoma protein expression profiles on 412 cases from TCGA and constructed a PRotein-driven index of OVARian cancer (PROVAR). PROVAR significantly discriminated an independent cohort of 226 high-grade serous ovarian carcinomas into groups of high risk and low risk of tumor recurrence as well as short-term and long-term survivors. Comparison with gene expression–based outcome classification models showed a significantly improved capacity of the protein-based PROVAR to predict tumor progression. Identification of protein markers linked to disease recurrence may yield insights into tumor biology. When combined with features known to be associated with outcome, such as BRCA mutation, PROVAR may provide clinically useful predictions of time to tumor recurrence.
Ovarian cancer is the deadliest gynecologic cancer in the United States, with 22,240 new cases and 14,030 deaths in 2013 (1). High-grade serous cancer is the most common ovarian epithelial malignancy, accounting for approximately 70% of all cases of epithelial ovarian cancer (2). Most cases are diagnosed at an advanced stage, with this being a key contributor to an overall 5-year survival rate of less than 40% (3, 4). Although the initial response rate to standard surgery and platinum-based chemotherapy is high, 30%–40% of patients relapse within 12 months and do not respond to further platinum therapy (5, 6). Early detection of high-grade serous ovarian cancer is thus a key to reducing morbidity and mortality from ovarian cancer (7). Several clinicopathologic factors, such as age, stage, histologic grade, and tumor residuum, are considered prognostic indicators in patients with ovarian cancer, but these factors are used only in a small number of patients to guide treatment decisions, due to insufficient sensitivity and specificity (8).
With the development of microarray technologies, several studies have identified genetic markers or gene expression profiles that are associated with the prognosis of high-grade ovarian cancer (9–12). However, these signatures often contain large numbers of genes, which reduces their applicability in clinical practice. Importantly, despite the significant association of gene signatures with overall survival (OS), their predictive value of treatment response and time to tumor recurrence is limited.
The reverse-phase protein arrays (RPPA) platform allows high-throughput measurements of protein expression levels in a large number of samples. RPPA profiles have been successfully used to identify protein markers of pharmacological response and to predict prognosis in breast cancer (13, 14).
Here, we used the RPPA technology to define a PRotein-driven index of OVARian cancer (PROVAR) and show that it is able to predict time to recurrence in an independent validation cohort, outperforming several gene expression–based approaches. Our work illustrates the potential of protein-driven treatment response predictions.
Patient characteristics. Patient characteristics are described in Table 1, and detailed patient information is provided in Supplemental Table 4; supplemental material available online with this article; doi:
10.1172/JCI68509DS1. All patients included in this study had serous epithelial ovarian carcinoma. More than 95% of tumors were classified as high grade (G2 or G3). Approximately 60% of patients in TCGA and 40% of patients in the validation set underwent optimal surgical cytoreduction (<1 cm residual disease at the end of surgery). The median progression-free survival (PFS) for TCGA samples (14.9 months) was shorter than that for validation samples (19.4 months), and the difference was statistically significant (log-rank test, P < 0.001; see Supplemental Figure 1). There was no statistically significant difference in OS between TCGA and validation sets, although the P value was trending toward significance (log-rank test, P = 0.101; Supplemental Figure 1).
Clinical characteristics of the training and validation sets
Identification and validation of protein markers. The logic flow chart shown in Figure 1 summarizes the procedure used to construct and validate a protein-based index of PFS; a comparison is also shown of the protein-driven model and several gene-driven models.
Flow chart for construction and validation of PROVAR and for comparison with gene-driven models from Konstantinopoulos et al. (10), Kang et al. (11), and Verhaak et al. (12).
We used the least absolute shrinkage and selection operator (lasso) to identify protein markers most associated with PFS. After applying an L1-constrained Cox regression to the 222 TCGA samples with nonmissing annotation on PFS, with a tuning parameter chosen by 10-fold cross-validation, we identified 9 protein markers significantly associated with PFS (Table 2 and Figure 2) and termed the predictive protein set PROVAR. The expression patterns of the 9 protein markers across ovarian carcinoma are shown in Figure 2. Five proteins (AR, BID, phosphorylated TAZ [pTAZ], phosphorylated EGFR [pEGFR], and HSP70) were associated with longer PFS in the Cox regression analysis and were highly expressed in the low-risk group, and 4 proteins (STAT5α, phosphorylated PKCα [pPKCα], phosphorylated MEK1 [pMEK1], and EEF2) were associated with shorter PFS and increased expression in the high-risk group (Table 2).
Expression signatures of the 9 protein markers for the training set (TCGA) and the validation set. The first 5 proteins (AR, BID, pTAZ, pEGFR, and HSP70) were associated with longer PFS and the rest (STAT5α, pPKCα, pMEK1, and EEF2) were associated with shorter PFS (see Table 2).
Nine proteins composing PROVAR
PROVAR scores were defined as a linear combination of the expression levels of the 9 proteins weighted by the Cox regression coefficients, and those scores were calculated for each sample in the training set. Values ranged from –0.57 to 0.44 (median = –0.02) across the 222 samples and were distributed unimodally within this range. Based on PROVAR, the samples were classified into high-risk (i.e., high PROVAR score) and low-risk (i.e., low PROVAR score) groups. The median score (–0.02) was used as a cutoff.
Kaplan-Meier survival analysis showed a significant difference in PFS between the high- and low-risk groups in the training set (n = 222, P = 0.001; Figure 3A). While PROVAR was constructed based on PFS, we additionally performed Kaplan-Meier analysis of OS in order to assess the capacity of PROVAR to differentiate patients by OS. The data set used in this analysis was comprised of 387 patients with available OS data, including 165 patients not used to construct the model. A significant difference in OS was found between risk groups (n = 387, P < 0.001; Figure 3B) but also when limiting the analysis to the 165 OS-only patients not included in the training set (n = 165, P = 0.032; Supplemental Figure 2). The same cutoff was used as in the analysis of PFS. The median PFS times of the high- and low-risk groups were 12.5 months and 18.3 months, respectively, and median OS times of the 2 risk groups were 36.9 months and 52.8 months, respectively.
Kaplan-Meier curves for the training set (TCGA) and the validation set after classifying patients into 2 groups by median split using PROVAR. (A) TCGA with PFS, (B) TCGA with OS, (C) validation set with PFS, and (D) validation set with OS. Red and black lines indicate high- and low-risk groups, respectively.
Unlike gene expression profiles, in which the number of features vastly outnumbers the size of the sample cohort, the limited number of features included on the protein array reduces the risk of overfitting when correlating expression with outcome (15). To further reduce the possibility of a significant difference in outcomes as a result of an overfitted predictive model, we sought to validate the performance of PROVAR in an independent sample cohort and performed RPPA profiling on a set of 226 ovarian cancer samples that were not included in TCGA. For each patient, PROVAR was computed and patients were classified into 1 of the 2 risk groups according to the index. The same cutoff was used as that found while analyzing the TCGA set to separate the validation set into high-risk and low-risk groups. In the validation set, median PFS times were 15.1 months and 24.0 months, respectively, between high- and low-risk groups (Figure 3C), while median survival time of the high-risk group was 39.6 months. The median survival was for the low-risk group exceeded 60 months (Figure 3D).
To test whether this result was independent of known predictive variables, such as age and stage, we applied multivariate analysis using a Cox proportional hazards model with PROVAR and clinicopathologic factors (age, stage, grade, and surgery status) as covariates. PROVAR was the only factor consistently significant for both OS and PFS across different data sets (Table 3). Age was a significant factor for the TCGA set but not for the validation set. Stage was not significant for any set, and more advanced grade was associated with a decrease in PFS in TCGA samples but not with others. Surgery status was found to be predictive of PFS or OS in the validation set but not in the TCGA set. BRCA1/2 mutation was a significantly favorable factor of progression and survival in the training set, but no mutation status data were available for our validation set. A 5% of significance level was used for all tests.
Univariate and multivariate Cox proportional hazards model analyses
Comparison with gene-driven models. For the purpose of comparison with the protein-driven model, we implemented 3 gene expression–based models and tested their ability to predict outcome on the subset of validation samples for which both protein and gene expression data were available (n = 130, Konstantinopoulos et al., ref. 10; Kang et al., ref. 11; Verhaak et al., ref. 12). The difference in outcome of high-risk and low-risk groups was compared between gene-based models and PROVAR. Kaplan-Meier curves were generated for all comparisons (Supplemental Figure 3A for PFS; Supplemental Figure 3B for OS). All 3 gene-based models failed to produce a difference in PFS or OS, whereas PROVAR classification resulted in a near significant classification of PFS (P = 0.056) as well as OS (P = 0.006). To address possible batch effects when comparing models based on different data sets, cutoff values for identifying 2 groups were determined for each model separately.
Additionally, we considered a 3-group classification using the 20th and 80th percentiles as cutoffs, focusing on more extreme cases in the entire sample. Overall, trends in progression or survival patterns remained stable in the validation set (Supplemental Figure 4A). On the subset of validation samples for which both protein and gene expression data were available (n = 130), no statistically significant difference in PFS was found with any of the 4 models, but there was an improvement in prediction of OS using PROVAR or CLassification of OVARian cancer (CLOVAR) (Supplemental Figure 4, B and C). However, we note that the relatively small sample size of each group may have limited the statistical power of the 3-group classification, because the log-rank test is based on large sample approximations.
Comparative analysis of PROVAR at proteomic and genetic levels. We constructed a new index using the genes matching the 9 proteins that constructed PROVAR and using the PROVAR coefficient weights to evaluate the predictive power of the gene expression levels of PROVAR proteins in the TCGA set. The median index value, based on the matching genes, was used as a cutoff to separate the data set into 2 groups. The index at the genetic level did not preserve the predictive power found when using PROVAR at the proteomic level (Supplemental Figure 5; log-rank test, P values were 0.199 for PFS and 0.173 for OS).
Robustness of protein marker index. In an attempt to test the robustness of the 9 protein markers, the validation set was used to identify protein markers with the lasso and then those markers were compared with the original 9 markers. By using the validation set, we found 7 protein markers, FOXO3A (0.008), CDK1 (0.003), VASP (0.055), CYCLINB1 (0.042), pNFKB (–0.034), YAP (0.025), and AR (–0.037) (numbers in parentheses are lasso coefficients from Cox’s regression). The only protein from the 9 TCGA protein markers to overlap with the 7 validation set protein markers was AR. We identified 4 major protein clusters using a hierarchical clustering method on RPPA data for TCGA samples (Figure 4). Interestingly, both the 9-protein signature and the 7-protein set included representative proteins from each of the 4 clusters.
The heat map of RPPA data for the TCGA set (222 samples and 172 proteins). The row dendrogram indicates unsupervised hierarchical clustering of 172 proteins. Colored bars next to the dendrogram represent (a) 8 unique markers found by using the TCGA set; (b) the overlap, AR; and (c) 6 unique markers found by using the validation set. The clusters were cut into 4 groups as indicated by different colors. Ward’s linkage method and absolute distance were used as a group linkage method and a distance measure, respectively. The column dendrogram indicates unsupervised hierarchical clustering of 222 patients. Red bars below the dendrogram indicate the high-risk group. The sample size of each cluster and a comparison with existing molecular subtypes are provided in Supplemental Table 2.
Additionally, we wanted to assess how similar the patient classifications were using the 9- and 7-protein signatures. The patients in the training set were classified by each of 2 signatures, and, more specifically, they were split into 2 risk groups using the median as a cutoff for each case (Supplemental Table 1). The highly similar grouping found when classifying the training set using the 9-protein set and the 7-protein set suggested a significant association between the 2 signatures (χ2, P < 0.001).
Comparison with molecular subtypes of ovarian cancer. Four clusters of samples were defined by hierarchical clustering on RPPA data for the TCGA set (Figure 4). Cross-tabulation (shown in Supplemental Table 2-1) was used to correlate the 4 clusters with 4 existing molecular subtypes (12, 16). To test whether there was significant enrichment in the overlap, we performed a χ2 test, which showed a statistically significant association between RPPA clusters and gene expression subtypes (P = 0.001). Interestingly, cluster 1 mostly consisted of samples from the mesenchymal gene expression subtype, which were entirely classified by PROVAR as high risk. Consequently, a decreased PFS was observed in cluster 1 compared with that in other clusters (Supplemental Figure 7).
We have previously reported that high-grade serous ovarian cancer samples often exhibit multiple gene expression–based subtype signatures and that the classification into different mutually exclusive subtypes may therefore be less informative than in other cancers (12). Hence, we examined, for each subtype signature, whether the signature occurred more frequently in certain clusters. As shown in Supplemental Table 2-2, both mesenchymal and immunoreactive signatures were highly activated in protein cluster 1.
Sample cluster 4 included a subset of 36 samples, which appeared to form a separate group, that we classified as a low-risk group. We identified 21 significant proteins that were differentially expressed between 2 groups (36 low-risk samples vs. 43 remaining samples in cluster 4), using t tests and multiple testing correction at a significance level of 1%. Among 21 differentially expressed proteins, 4 were included in the PROVAR signature (AR, BID, HSP70, and pEGFR), and those 4 proteins were consequently highly expressed in the low-risk samples. This subcluster was potentially a major driver of the PROVAR signature.
Biological interpretation of PROVAR. To characterize the biological properties of PROVAR, we used the FatiGO tool to associate the index with gene ontology (17) as well as BioCarta (
http://www.biocarta.com/genes/index.asp) annotations. In this analysis, the annotation of proteins within the signature is compared to the annotation of other proteins on the array. Gene ontology terms and BioCarta pathways in which PROVAR proteins were most significantly overrepresented are listed in Supplemental Table 3 (P < 0.05). The top-ranked BioCarta pathway was the EGF signaling pathway, in which the EGF binds to the EGFR in the cell membrane, starting a cascade of signals. The phosphate signal activates MAP kinases (also known as ERK), which leads to mRNA transcription in the cell nucleus (18, 19). Indeed, pEGFR and 3 proteins representing different pathways downstream of EGFR are included in the PROVAR signature: STAT5α, pMEK1, and pPKCα. However, whereas protein expression of pEGFR was associated with favorable outcome, increase in protein levels of STAT5α, pMEK1, and pPKCα correlated with high risk. To gain more insights in the general patterns of protein expression, we applied hierarchical clustering to group proteins by similarity of expression in the 222 samples of the training set. We observed 4 large protein clusters, each representing a general expression trend (Figure 4). Interestingly, the 4 proteins included in the PROVAR signature as well as the EGF pathway were found to represent each of the 4-protein expression clusters. This suggests that protein expression of pEGFR, STAT5α, pMEK1, and pPKCα does not correlate and therefore does not reflect a single cellular pathway but possible reflects 4 different cellular mechanisms, with different mechanisms activating and inactivating each pathway.
Characterization of PROVAR signatures by genomic aberrations. Using Welch’s t tests and Mann-Whitney tests, we evaluated whether a difference in PROVAR scores was associated with specific genomic aberrations. Included in the analysis were TCGA samples for which PROVAR scores and genomic data were available (copy number alterations, n = 376 samples; mutation events, n = 208 samples) and genes with alterations in more than 3% of samples (1,248 focally deleted genes; 10,181 highly amplified genes; 34 mutated genes). After multiple testing correction, no copy number deletion, amplification, or gene mutation was found to result in significantly different PROVAR scores.
Additional predictive value of other prognostics factors. We investigated the effects of a combination of PROVAR and other available factors, such as age, stage, grade, surgery status, and BRCA1/2 mutation (BRCA mutation status was only available in the training set). A total of 130 TCGA samples and 208 validation samples had complete data for analysis. Three multivariate Cox models were considered with the following as covariates: PROVAR only (model A); PROVAR plus age, stage, grade, and surgery status (model B); PROVAR plus age, stage, grade, surgery status, and BRCA1/2 mutation (model C). For TCGA samples, we did not observe a significant improvement in the 2 extended models (models B and C), in terms of the differentiation of survival rates between groups, compared with the PROVAR only model (Supplemental Figure 8). Similarly, model B, with the addition of clinical factors to model A, did not add significant predictive power in the validation set.
In this study, we identified 9 protein markers for predicting time to recurrence using the protein expression data on 222 TCGA primarily high-grade serous ovarian cancers and developed a PROVAR risk classification system and successfully validated its discriminative ability to predict both PFS and OS in an independent validation set of 226 high-grade serous ovarian carcinomas. Unfortunately, no targeted therapies are currently approved for treatment of ovarian cancer, and gynecologic oncologists cannot choose a specific treatment for a subpopulation of patients with ovarian cancer like HER2 inhibitors for patients with HER2+ breast cancer. In order to improve outcome in ovarian cancer, therapies that are potentially more effective but with increased rates of treatment-related adverse events, such as dose-dense chemotherapy (20) or combination of taxane, platinum, and bevacizumab (21, 22), may be administered to those patients identified as at high risk of early recurrence. A prediction approach based on a small number of proteins that can be measured using immunohistochemistry may be more practical and accurate in clinical management than methods based on gene expression linked to OS (8, 23). Validation of PROVAR in a prospective setting is required before it can be translated to clinical use.
We found a significant association between PFS (or OS) and the classification of patients by PROVAR in an external validation set. The training set (TCGA) and the validation set were found to be heterogeneous in terms of survival, which emphasizes the relevance of this result. In contrast, several recently proposed gene-driven models for predicting survival were unable to demonstrate significance for the prediction of PFS particularly in the validation set. (10–12). OS in ovarian carcinoma is influenced by various factors, such as differences in therapy after progression. Our model presents an advantage in being constructed based on PFS. We propose that including additional factors, such as BRCA mutation status, can further improve the accuracy of our predictions.
Statistical challenges of high dimensionality arise when predicting survival based on genomics or proteomics data when extracting meaningful statistical and biological information from high-dimensional data in which each biological sample is assessed with a large number of measurements (24). In this study, the lasso was used for variable selection and coefficient estimation simultaneously. The lasso is known to improve the overall prediction accuracy by sacrificing a little bias to reduce the variance of the predicted value (25). The highly correlated nature of genomic and proteomic data has led to methods that incorporate the correlation structure among variables into a prediction model and to select a set of complementary features (24). One such approach is the group lasso, a generalization of the lasso that was designed to select the same features in all discriminant directions using a mixed-norm regularization over regression coefficients (26).
The robustness of the PROVAR signature was investigated through comparison of the 9 markers derived from the TCGA set and the 7 markers identified from the validation set. Although only the androgen receptor (AR) protein was found in both signatures, a similarity in protein expression profiles was observed between the 2 sets of markers by hierarchical cluster analysis. Presurgical androgen-deprivation therapy is first-line therapy in nonmetastatic prostate cancer (27), and the impact of AR expression in ovarian carcinoma has been the subject of several studies (28–30). AR expression was a favorable prognostic factor for survival, and malignant transformation was found to involve downregulation of AR in patients with epithelial ovarian cancer. No differences in AR expression were observed between primary and metastatic ovarian tumors, suggesting that downregulation of AR occurs during early carcinogenesis (29). Transcriptome-based predictive analysis of treatment response in breast cancer has been shown to be highly accurate by multiple gene signatures. Despite the very limited overlap between these signatures, they generally identified the same groups of responsive and nonresponsive breast tumors (31).
Interestingly, AR has been found to modulate EGFR signaling and to promote an aggressive phenotype by increasing the abundance of phosphorylated EGFR and MAPK (30). Gene ontology analysis of the protein signature further confirmed the possible importance of the pEGFR and MEK/ERK signaling cascade. STAT5α, pMEK1, and pPKCα proteins, which are representative markers of 3 intracellular signaling (STAT signaling, PKC signaling, and MEK/ERK signaling) pathways, play especially important roles associated with cell proliferation, survival, and cell invasion in cancer cells and the activation of these signaling pathways would be expected to be associated with shorter PFS/OS time in patients with ovarian cancer. Although inhibition of EGFR has not shown clear efficacy in treatment of ovarian carcinoma (32), MEK inhibitors have been reported to control low-grade disease, and combination therapy may provide a possible avenue toward improving outcome of patients with ovarian cancer (33, 34).
In conclusion, PROVAR is simple but predictive of time to progression and survival, making it potentially useful in clinical practice. Application of PROVAR could result in more direct and reliable insights into the mechanisms of recurrence after platinum-based chemotherapy in patients with high-grade serous ovarian cancer. Integrating the protein signature with genetic mutations associated with survival will likely result in the most optimal predictive model of outcome. In the era of personalized medicine, identification of patients at high risk of early recurrence may provide clinicians with opportunities for early interference and positively impact survival for a group of patients in dire need of improved prospects.
Patient samples. Fresh-frozen samples were collected from newly diagnosed patients with serous ovarian carcinoma who were undergoing surgical resection prior to chemotherapy. Ovarian tumors that were diagnosed as of the serous histologic type, with at least 70% tumor purity per pathology review, were used in the subsequent analyses (9, 16). Specifically, for TCGA project, tumor specimens were collected as follows: (a) each tumor specimen was approximately 1 cm3 in size and weighed between 100 mg and 200 mg; (b) each specimen was embedded in optimal cutting temperature (OCT) medium, and histologic sections were obtained from top and bottom portions for review. As is the custom for tissues admitted to TCGA, all sections were reviewed by a pathologist from both the tissue source site and TCGA’s biospecimen core resource, and histology and cellularity were confirmed independently. In the validation set, tissue samples were collected and snap frozen in liquid nitrogen immediately after surgical resection. In addition, adjacent tissues were formalin fixed and paraffin embedded according to standard procedures at the time of primary surgery. The histological characteristics of all specimens were assessed on hematoxylin and eosin–stained sections by pathology review.
RPPA. RPPA is a high-throughput antibody-based technique for simultaneously measuring protein expression levels in a large number of biological samples (13, 35). For this study, RPPA analyses were performed at the RPPA core facility at MD Anderson Cancer Center using standard operating procedures and highly validated antibodies (35). Lysates were prepared from the frozen tumor samples and spotted in a serial dilution onto nitrocellulose-coated slides. Each slide was then probed with a validated primary antibody followed by a secondary antibody. The signal obtained was amplified using a Dako Cytomation–catalyzed system (Dako) and visualized by DAB colorimetric reaction. The stained slides were scanned, analyzed, and quantified using MicroVigene software (VigeneTech Inc.). Finally, a 3-parameter logistic model was assumed for the dependency of the observed intensity on the unknown protein expression, and protein expressions were estimated iteratively with the logistic curve parameters using the least-squares method. This algorithm is implemented in the R package SuperCurve (36, 37) (publicly available at
http://bioinformatics.mdanderson.org/Software/supercurve). To correct for batch and sample effects, we applied Z-normalization to RPPA data derived from the training and validation sets before subsequent statistical analysis. Z-normalization was performed by sample-wise median centering, protein-wise median centering, and protein-wise scaling (division by the standard deviation).
Validation of antibodies. Antibodies were validated by Western blot, and any antibody that did not produce a single predominant band was excluded. For proteins whose expression did not show a sufficiently dynamic range to facilitate antibody validation, siRNA was used to manipulate the signal to allow evaluation of RPPA-immunoblotting correlations. Details on the validation process of antibodies and the results on the set of antibodies included in RPPA assays have been described elsewhere (38). Results of Western blot analyses for our 9 protein markers, which compose PROVAR, using the same antibodies and using lysates from several cancer cell lines, are provided in Supplemental Figure 6.
Seven markers except pTAZ and HSP70 showed monospecific bands with high correlation coefficients between RPPA and Western blot–based protein expressions. Although HSP70 had a low correlation coefficient, it showed a distinct monospecific band. The low correlation was thought to be purely due to a narrow dynamic range of expression levels. Despite the relatively poor validation, pTAZ antibody was used for RPPA, because pTAZ is one of the markers that we believe may represent damaged tumor samples, and, thus, pTAZ antibody, more than being a specifically pTAZ-targeting antibody, is a marker for damaged tumor samples. Also, as shown in Table 2, pTAZ has a fairly significant association with PFS, which has led us to retain pTAZ in our protein pool.
RPPA profiling of TCGA samples. Quantitative protein expression profiles, consisting of 172 proteins and phosphoproteins, were generated using RPPA (Supplemental Table 5) and are now available at the TCGA portal (
http://tcga-data.nci.nih.gov/tcga/). In total, 412 serous epithelial ovarian cancer samples from TCGA were included in our data cohort. After excluding stage I patients, 222 cases had no missing values for PFS and 387 cases had no missing values for OS. Patients with no missing values were included in the PFS analysis and OS analysis. Clinical information was downloaded (as of February 2, 2012) from the TCGA portal (
http://tcga-data.nci.nih.gov/tcga/). PFS was defined as the time from initial surgery to the first documented progression or recurrence or the last follow-up in the absence of progressive disease. OS was defined as the time between surgery and the last follow-up or cancer-related death. Both PFS and OS were capped at 60 months as described previously (16). The protein expression values, normalized by using Z-normalization (see RPPA), for the 387 samples are provided in Supplemental Table 6.
RPPA profiles of validation cohort. An independent data set of patients with advanced-stage serous carcinoma was obtained from Philadelphia, Pennsylvania, USA, and Japan, and the expression levels of overlapping proteins and phosphoproteins were measured by RPPA. In total, 144 proteins and phosphoproteins were in common between RPPA data sets used for training and validation. All Japanese patients provided written informed consent for the collection of samples and subsequent analysis. All Japanese patients were treated with taxane and platinum-based chemotherapy as adjuvant chemotherapy. All patients undergoing surgery, chemotherapy, or other treatment for malignant gynecological and other tumors and patients treated for benign gynecologic conditions were eligible to donate biospecimens.
Potential batch effects were adequately adjusted by using Z-normalization (see RPPA). Tumor patients with nonmissing data on PFS (209 patients) or OS (226 patients) were used for validation of predictive models of PFS or OS. In Supplemental Tables 7-1 and 7-2, the adjusted RPPA data are presented, for all validation samples.
Gene expression microarray data. Gene expression data were extracted to perform microarray analysis of gene expression profiles in the TCGA set for the purpose of comparison with protein markers. Level 3 microarray data in TCGA were obtained from the TCGA portal and were generated as described previously (16). The data were obtained from the Agilent platform and were lowess normalized. At the time of access, 270 patient samples with both gene expression data and corresponding PFS times and 500 patients with gene expression data and OS times were available for analysis. In total, 378 samples had matching RPPA and gene expression profiles.
Gene expression data and clinical annotation are publicly available for a subset of validation samples (n = 130), and the data were downloaded from the GEO database (
http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32062). The data were obtained using the Agilent platform and were normalized in GeneSpring GX by setting threshold raw signals to 1.0 and using the percentile shift algorithm to the 75th percentile.
Prediction of survival. In order to identify protein markers most associated with PFS, we applied an L1-constrained (lasso) Cox regression via cyclical coordinate descent algorithm (39) to the RPPA data in the training set (n = 222). The lasso produces sparse interpretable models by shrinking some variables to exactly zero (40). The coefficients (β) in Cox’s regression model are estimated by maximizing the partial likelihood function subject to a constraint on the L1-norm of the coefficients. The lasso estimator (β̂) maximizes the objective function given below:
l(β) – λβ1
Here l(β) is the log partial likelihood in the Cox model, and for the exact form of this function, see reference 41. The tuning parameter, λ, in Equation 1 was chosen by 10-fold cross-validation. For the implementation, we used the R package “glmnet” (39).
PROVAR was defined for each of the 222 TCGA samples as the sum of the estimated coefficients multiplied by protein expression levels, as shown below. Here i represents patients (i = 1,...,222), j represents proteins with non-zero coefficients (j = 1, ..., m), βj is the lasso coefficient of the jth protein marker, and Xij is the expression level of the jth protein for the ith patient.
PROVA = j = 1mβjXij
Specifically, PROVAR was defined as follows: –0.096 × AR –0.075 × BID + 0.029 × EEF2 –0.052 × pEGFR –0.027 × HSP70 + 0.025 × pMEK1 + 0.022 × pPKCα + 0.004 × STAT5α –0.058 × pTAZ. Each patient was classified into 1 out of 2 risk groups according to the index (PROVAR), with the median value as a cutoff.
The test set was used to validate PROVAR. For each patient, PROVAR was computed in a similar way to Equation 2, and then patients were classified into 1 out of the 2 risk groups based on the index with the same cutoff value as for the training set. The log-rank test was used to assess the patient risk classification.
In order to evaluate PROVAR as an independent predictor of progression and survival, we considered univariate and multivariate Cox proportional hazards models with PROVAR and clinicopathologic factors (age, stage, grade, and surgery status) as covariates. Among these, age, stage, and grade were treated as ordinal continuous variables with natural ordering. Further, risk scores were recalculated using regression coefficients from multivariate models, and patient classification by median split was evaluated by the log-rank test.
Gene-driven models. Restoration of DNA repair has been known to promote resistance to DNA-damaging chemotherapeutic agents (42, 43). In recent years, this mechanism has attracted attention (see Konstantinopoulos et al., ref. 10; Kang et al., ref. 11). Also it was reported that a portion of ovarian tumors are characterized by high numbers of infiltrating T-lymphocytes and stromal cells, and those tumors can be recognized by expression profiling (44, 45). A very recent work extended this result and identified subtypes and survival gene expression signatures that could provide a prognostic model (Verhaak et al., ref. 12). All these 3 models were gene driven, and in comparison with our protein-driven model, the validity of the gene-based models was assessed using the subset of validation samples with matching protein and gene expression data (n = 130).
In Konstantinopoulos et al., the index of BRCA1/2-deficient phenotype was defined as the weighted sum of the expression levels of 60 genes (10). If the index was greater than a certain threshold, the sample was classified as BRCA-like (BL), and otherwise, the sample was classified as non-BRCA-like (NBL). It was concluded that patients with the BL profile had improved PFS and OS compared with patients with a NBL profile. Using the 60 genes and the weights provided in their study, we were able to compute the index of BRCA1/2-deficient phenotype for the subset of validation samples (n = 130). Due to potential microarray platform-specific effects with different scales and distributions of expression levels, we used the median index value of BRCA1/2-deficient phenotype over the 130 patients as a cutoff when dividing patients into BL or NBL groups, instead of using the threshold proposed in Konstantinopoulos et al. (10). It was found in their study that 20 (29%) of the 70 patient validation cohort demonstrated the BL profile, and, similarly, we used the 70th percentile of index value as an additional cutoff so that the proportion of BL samples to NBL samples would be 3:7.
In Kang et al., a score was developed based on the expression levels of 23 genes in DNA repair pathways (11). Per patient, a point was given for each gene associated with better survival if the expression was higher than the median expression and, on the other hand, a point was given for each gene associated with poorer survival if the expression was lower than the median expression across samples. The score was then defined for each patient as a simple sum of points. An assessment of whether each of the 23 genes was associated with better or poorer survival was provided in their article. Patients were classified into high-risk (i.e., low score) groups and low-risk (i.e., high score) groups using the median score as a cutoff. They reported that a significantly improved survival was observed in the high-scoring group.
In Verhaak et al., a patient classification system, referred to as CLOVAR, was developed based on subtype and survival gene expression signatures (12). After selecting 100 genes whose expressions were most correlated or anticorrelated with OS, single-sample gene set enrichment analysis was performed to generate gene set activation scores for CLOVAR subtype and CLOVAR survival signatures. To adjust for batch effects, the scores were normalized by conversion into the [0, 1] interval within each data set. Each sample was then classified into 3 subtypes (immunoreactive, mesenchymal, and others) according to the normalized subtype scores with cutoffs specified (12). To combine CLOVAR subtype and CLOVAR survival signatures, a multivariate Cox proportional hazards model was considered with subtype classification and CLOVAR survival score as covariates. A risk score was computed for each sample, defined as a linear combination of covariates weighted by the Cox regression coefficients. Patients were classified into high-risk and low-risk groups using the median risk score as a cutoff. It was reported in this study that this patient classification significantly predicted OS (12). We were able to apply this classification scheme to our validation samples (n = 130) using the same Cox regression coefficients as found in Verhaak et al. (12).
PROVAR at the genetic level. In order to see whether matching genes showed similar expression patterns to the protein markers, we examined whether PROVAR would preserve its predictive power at the genetic level by using expressions of the respective genes. A new index was generated by replacing expression levels of proteins with those of matching genes and retaining the same lasso coefficients. The cutoff was accordingly set at the median index value based on the matching genes.
Robustness of protein markers. In an attempt to test the robustness of the protein markers identified in the TCGA set, we used the validation set in reverse to identify protein markers using the lasso, and then we studied how much the 2 sets of protein markers were similar. Hierarchical clustering was used for this purpose. A χ2 test was used to determine the similarity of the patient classifications using the 9- and 7-protein signatures.
Comparison with molecular subtypes of ovarian cancer. A previous study (16) identified 4 main molecular subtypes and signatures of ovarian cancer. To examine similarity between gene and protein expression–based clustering, we performed hierarchical clustering to identify sample clusters using RPPA protein expression data and correlated the protein expression–based hierarchical clusters with the previously reported gene expression–based 4 subtypes. To facilitate comparison, 4 clusters were generated by hierarchical clustering. A cross tabulation and χ2 test of independence were used to determine whether a significant association existed between them.
In our previous study (12), we reported that high-grade serous ovarian cancer samples do not consist of mutually exclusive expression subtypes but that samples often exhibit multiple signatures at different levels of activation. We therefore examined, for each subtype, whether the signature occurred more frequently in certain clusters.
Functional annotation of protein markers. To determine the functional linkages among our protein markers against the rest of the proteins in the TCGA set, we compared gene ontology of the PROVAR proteins with the remaining proteins of the RPPA platform using FatiGO (
http://bioinfo.cipf.es/babelomicswiki/tool:fatigo; refs. 46, 47). It extracts gene ontology terms or pathways that are significantly overrepresented in the gene set (or protein set) of interest compared with the genes of reference (or proteins). The reference is usually the rest of genes involved in the experiment. The statistical significance was evaluated by means of a Fisher’s exact test for 2 × 2 contingency tables.
Associations of PROVAR signatures with genomic aberrations. We investigated whether PROVAR was significantly associated with genomic alterations by comparing PROVAR scores of altered samples with PROVAR scores of nonaltered samples. Segmented copy number data for 559 TCGA ovarian carcinomas were processed using the GISTIC2.0 pipeline, and genes were categorized into 1 out of 3 groups as previously described (48): (a) deep loss, (b) neutral, (c) high-level gain. Exon sequence data on 18,500 genes and 426 TCGA samples were available for somatic mutations. Two-sample Welch’s t tests and Mann-Whitney tests were performed on genes that showed alterations in at least 3% of samples. Correction for multiple testing was carried out using Benjamini and Hochberg (49) or Benjamini and Yekutieli (50).
Combination of PROVAR with other factors. To see whether combining PROVAR with known prognostic factors would enhance the performance in prediction, we assessed the additional predictive value of available factors, such as age, stage, grade, surgery status, and BRCA1/2 mutation. We considered 3 models with the following as covariates: PROVAR only (model A); PROVAR plus age, stage, grade, and surgery status (model B); and PROVAR plus age, stage, grade, surgery status, and BRCA1/2 germline/somatic mutation (model C).
A total of 130 TCGA samples were available for testing of models A, B, and C, whereas 208 validation samples were included for evaluation of models A and B. BRCA1/2 mutation status was unavailable for the validation set. To allow comparison, Cox regression coefficients were estimated using 130 TCGA samples and all 3 models as well as 208 validation samples for models A and B. Patients were classified into 2 risk groups according to their risk scores (i.e., predicted values derived from multivariate Cox regressions), with cutoff points at the median risk score in each model.
Statistics. Statistical analysis was performed using the R statistical computing environment. The lasso was used to identify protein markers that were significantly associated with PFS. Univariate and multivariate Cox proportional hazards regression models were used with log-rank tests and Wald’s tests to assess differences in survival, and a P value of less than 0.05 was considered significant.
Study approval. TCGA ovarian cancer samples were collected and processed through TCGA Biospecimens Core Resource at the International Genomics Consortium (Phoenix, Arizona, USA) (as described in ref. 16). All specimens were collected using IRB-approved protocols and were deidentified to ensure patient confidentiality. Prior written informed consent was obtained from all Japanese patients, and the study protocol was approved by Institutional Review Boards at Niigata University (no. 239), Kagoshima City Hospital (no. H19-21), Nagasaki University (no. 080509), Tottori University (no. 1349), Kurume University (no. 78), and National Defense Medical College (no. 757).
View Supplemental data
View Supplemental Excel file 1
The results published here are in whole or part based upon data generated by the TCGA pilot project established by the National Cancer Institute (NCI) and NHGRI. Information about TCGA and the investigators and institutions who constitute TCGA research network can be found at
http://cancergenome.nih.gov/. This work is supported by NCI award no. 5 P50 CA083639-12 (to R.G.W. Verhaak), NCI grant no. CA143883 (MD Anderson Genome Data Analysis Center) and P50 CA083639 (to G.B. Mills), and MD Anderson RPPA Core Facility (supported by NIH grant CA016672).
Conflict of interest: The authors have declared that no conflict of interest exists.
Citation for this article:J Clin Invest. 2013;123(9):3740–3750. doi:10.1172/JCI68509.
Siegel R, Naishadham D, Jemal A. Cancer statistics, 2013. CA Cancer J Clin. 2013;63(1):11–30.
Kobel M, et al. Differences in tumor type in low-stage versus high-stage ovarian carcinomas. Int J Gynecol Pathol. 2010;29(3):203–211.
Jemal A, Siegel R, Xu J, Ward E. Cancer statistics, 2010. CA Cancer J Clin. 2010;60(5):277–300.
Coleman MP, et al. Cancer survival in Australia, Canada, Denmark, Norway, Sweden, and the UK, 1995–2007 (the International Cancer Benchmarking Partnership): an analysis of population-based cancer registry data. Lancet. 2011;377(9760):127–138.
Winter WE 3rd, et al. Prognostic factors for stage III epithelial ovarian cancer: A Gynecologic Oncology Group Study. J Clin Oncol. 2007;25(24):3621–3627.
du Bois A, Reuss A, Pujade-Lauraine E, Harter P, Ray-Coquard I, Pfisterer J. Role of surgical outcome as prognostic factor in advanced epithelial ovarian cancer: a combined exploratory analysis of 3 prospectively randomized phase 3 multicenter trials: by the Arbeitsgemeinschaft Gynaekologische Onkologie Studiengruppe Ovarialkarzinom (AGO-OVAR) and the Groupe d’Investigateurs Nationaux Pour les Etudes des Cancers de l’Ovaire (GINECO). Cancer. 2009;115(6):1234–1244.
Jelovac D, Armstrong DK. Recent progress in the diagnosis and treatment of ovarian cancer. CA Cancer J Clin. 2011;61(3):183–203.
Agarwal R, Kaye SB. Prognostic factors in ovarian cancer: how close are we to a complete picture? Ann Oncol. 2005;16(1):4–6.
Yoshihara K, et al. High-risk ovarian cancer based on 126-gene expression signature is uniquely characterized by downregulation of antigen presentation pathway. Clin Cancer Res. 2012;18(5):1374–1385.
Konstantinopoulos PA, et al. Gene expression profile of BRCAness that correlates with responsiveness to chemotherapy and with outcome in patients with epithelial ovarian cancer. J Clin Oncol. 2010;28(22):3555–3561.
Kang J, D’Andrea AD, Kozono D. A DNA repair pathway-focused score for prediction of outcomes in ovarian cancer treated with platinum-based chemotherapy. J Natl Cancer Inst. 2012;104(9):670–681.
Verhaak RGW, et al. Prognostically relevant gene signatures of high-grade serous ovarian carcinoma. J Clin Invest. 2013;123(1):517–525.
Tibes R, et al. Reverse phase protein array: validation of a novel proteomic technology and utility for analysis of primary leukemia specimens and hematopoietic stem cells. Mol Cancer Ther. 2006;5(10):2512–2521.
Tabchy A, Hennessy BT, Gonzalez-Angulo AM, Bernstam FM, Lu Y, Mills GB. Quantitative proteomic analysis in breast cancer. Drugs Today (Barc). 2011;47(2):169–182.
Zhang Y, Rajapakse JC. Machine Learning In Bioinformatics.Hoboken, New Jersey, USA: John Wiley & Sons, Inc.; 2008.
Cancer Genome Atlas Research Networket al. Integrated genomic analyses of ovarian carcinoma. Nature. 2011;474(7353):609–615.
Harris MA, et al. The Gene Ontology (GO) database and informatics resource. Nucleic Acids Res. 2004;32(Database issue):D258–D261.
Nakakuki T, Ide K, Yumoto N, Nagashima T, Hatakeyama M. In silico simulation of epidermal growth factor signaling in prostate cancer cells. Presented at: Genomic Signal Processing and Statistics, 2007. GENSIPS 2007. IEEE International Workshop; June 2007.
Schulze WX, Deng L, Mann M. Phosphotyrosine interactome of the ErbB-receptor kinase family. Mol Syst Biol. 2005;1:2005.008.
Katsumata N, et al. Dose-dense paclitaxel once a week in combination with carboplatin every 3 weeks for advanced ovarian cancer: a phase 3, open-label, randomised controlled trial. Lancet. 2009;374(9698):1331–1338.
Burger RA, et al. Incorporation of bevacizumab in the primary treatment of ovarian cancer. N Engl J Med. 2011;365(26):2473–2483.
Perren TJ, et al. A phase 3 trial of bevacizumab in ovarian cancer. N Engl J Med. 2011;365(26):2484–2496.
Kim DC, Wang X, Yang CR, Gao JX. A framework for personalized medicine: prediction of drug sensitivity in cancer by proteomic profiling. Proteome Sci. 2012;10(suppl 1):S13.
Donoho DL. High-dimensional data analysis: The curses and blessings of dimensionality. Presented at: Aide-Memoire of a Lecture at AMS Conference on Math Challenges of the 21st Century; August 2000; Los Angeles, California, USA.
Zou H, Hastie T, Tibshirani R. On the “degrees of freedom” of the lasso. Ann Statist. 2007;35(5):2173–2192.
Yuan M, Lin Y. Model selection and estimation in regression with grouped variables. J R Stat Soc Ser B Stat Methodol. 2006;68(1):49–67.
Sharifi N, Gulley JL, Dahut WL. Androgen deprivation therapy for prostate cancer. JAMA. 2005;294(2):238–244.
Elattar A, et al. Androgen receptor expression is a biological marker for androgen sensitivity in high grade serous epithelial ovarian cancer. Gynecol Oncol. 2012;124(1):142–147.
Nodin B, et al. Increased androgen receptor expression in serous carcinoma of the ovary is associated with an improved survival. J Ovarian Res. 2010;3:14.
Li AJ, Scoles DR, Armstrong KU, Karlan BY. Androgen receptor cytosine-adenine-guanine repeat polymorphisms modulate EGFR signaling in epithelial ovarian carcinomas. Gynecol Oncol. 2008;109(2):220–225.
Prat A, et al. Concordance among gene expression-based predictors for ER-positive breast cancer treated with adjuvant tamoxifen. Ann Oncol. 2012;23(11):2866–2873.
Farley J, et al. Selumetinib in women with recurrent low-grade serous carcinoma of the ovary or peritoneum: an open-label, single-arm, phase 2 study. Lancet Oncol. 2013;14(2):134–140.
Bartholomeusz C, et al. MEK1/2 inhibitor selumetinib (AZD6244) inhibits growth of ovarian clear cell carcinoma in a PEA-15-dependent manner in a mouse xenograft model. Mol Cancer Ther. 2012;11(2):360–369.
Bookman MA, Darcy KM, Clarke-Pearson D, Boothby RA, Horowitz IR. Evaluation of monoclonal humanized anti-HER2 antibody, trastuzumab, in patients with recurrent or refractory ovarian or primary peritoneal carcinoma with overexpression of HER2: A phase II trial of the Gynecologic Oncology Group. J Clin Oncol. 2003;21(2):283–290.
Carey MS, et al. Functional proteomic analysis of advanced serous ovarian cancer using reverse phase protein array: TGF-beta pathway signaling indicates response to primary chemotherapy. Clin Cancer Res. 2010;16(10):2852–2860.
Yang JY, He X. A multistep protein lysate array quantification method and its statistical properties. Biometrics. 2011;67(4):1197–1205.
Hu J, He X, Baggerly KA, Coombes KR, Hennessy BT, Mills GB. Non-parametric quantification of protein lysate arrays. Bioinformatics. 2007;23(15):1986–1994.
Hennessy BT, et al. A technical assessment of the utility of reverse phase protein arrays for the study of the functional proteome in non-microdissected human breast cancers. Clin Proteomics. 2010;6(4):129–151.
Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for Cox’s proportional hazards model via coordinate descent. J Stat Softw. 2011;39(5):1–13.
Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B Stat Methodol. 1994;58:267–288.
Zou H. A note on path-based variable selection in the penalized proportional hazards model. Biometrika. 2008;95(1):241–247.
Yang D, et al. Association of BRCA1 and BRCA2 mutations with survival, chemotherapy sensitivity, and gene mutator phenotype in patients with ovarian cancer. JAMA. 2011;306(14):1557–1565.
Martin LP, Hamilton TC, Schilder RJ. Platinum resistance: the role of DNA repair pathways. Clin Cancer Res. 2008;14(5):1291–1295.
Tothill RW, et al. Novel molecular subtypes of serous and endometrioid ovarian cancer linked to clinical outcome. Clin Cancer Res. 2008;14(16):5198–5208.
Marquez RT, et al. Patterns of gene expression in different histotypes of epithelial ovarian cancer correlate with those in normal fallopian tube, endometrium, and colon. Clin Cancer Res. 2005;11(17):6116–6126.
Al-Shahrour F, Minguez P, Vaquerizas JM, Conde L, Dopazo J. BABELOMICS: a suite of web tools for functional annotation and analysis of groups of genes in high-throughput experiments. Nucleic Acids Res. 2005;33(Web Server issue):W460–W464.
Al-Shahrour F, et al. BABELOMICS: a systems biology perspective in the functional annotation of genome-scale experiments. Nucleic Acids Res. 2006;34(Web Server issue):W472–W476.
Cancer Genome Atlas Research Network Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature. 2008;455(7216):1061–1068.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Stat Methodol. 1995;57(1):289–300.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Statist. 2001;29(4):1165–1188.