Abstract
Purpose: Gleason score strongly predicts prostate cancer mortality; however, scoring varies among pathologists, and many men are diagnosed with intermediate-risk Gleason score 7. We previously developed a 157-gene signature for Gleason score using a limited gene panel. Using a new whole-transcriptome expression dataset, we verified the previous signature's performance and developed a new Gleason signature to improve lethal outcome prediction among men with Gleason score 7.
Experimental Design: We generated mRNA expression data from prostate tumor tissue from men in the Physicians' Health Study and Health Professionals Follow-Up Study (N = 404) using the Affymetrix Human Gene 1.0 ST microarray. The Prediction Analysis for Microarrays method was used to develop a signature to distinguish high (≥8) versus low (≤6) Gleason score. We evaluated the signature's ability to improve prediction of lethality among men with Gleason score 7, adjusting for 3 + 4/4 + 3 status, by quantifying the area under the receiver operating characteristic (ROC) curve (AUC).
Results: We identified a 30-gene signature that best distinguished Gleason score ≤6 from ≥8. The AUC to predict lethal disease among Gleason score 7 men was 0.76 [95% confidence interval (CI), 0.67–0.84] compared with 0.68 (95% CI, 0.59–0.76) using 3 + 4/4 + 3 status alone (P = 0.0001). This signature was a nonsignificant (P = 0.09) improvement over our previous signature (AUC = 0.72).
Conclusions: Our new 30-gene signature improved prediction of lethality among men with Gleason score 7. This signature can potentially become a useful prognostic tool for physicians to improve treatment decision making. Clin Cancer Res; 23(1); 81–87. ©2016 AACR.
See related commentary by Yin et al., p. 6
Translational Relevance
Prostate cancer is a very heterogeneous disease; many men experience an indolent course of disease, whereas others have aggressive disease that progresses to metastases and death. Gleason score is one of the best tools to predict the risk of lethal prostate cancer. However, its predictive ability is limited both by high interobserver variability among pathologists and the fact that many men are diagnosed with intermediate-risk Gleason 7 tumors. In this study, we developed a 30-gene mRNA expression signature that strongly differentiated Gleason score ≤6 from ≥8. When applied to men with Gleason score 7, the signature significantly improved prediction of lethal disease—that is, men who molecularly were more similar to Gleason score ≥8 were at a higher risk of dying from prostate cancer. This signature can be used clinically to help confirm the Gleason score and could become a useful additional prognostic tool to improve treatment decision making.
Introduction
One of the strongest predictors of lethal prostate cancer is Gleason score, which measures the extent of differentiation of prostate tumor tissue (1). The risk of prostate cancer–specific mortality with a Gleason score of 6 or below is minimal but increases dramatically with increasing Gleason score (2). In PSA-screened populations, approximately 40% of men diagnosed with prostate cancer have Gleason score 7, a group with intermediate risk of dying from prostate cancer (3). Gleason score 7 tumors may be divided into primary pattern 3 (3 + 4) and primary pattern 4 (4 + 3), with 4 + 3 tumors carrying a 3.1 times higher risk of death (2). However, treatment decisions remain difficult, as most men with Gleason score 7 do not die from their disease, even in the absence of treatment (4). Markers that improve prognosis prediction among these men are scarce, so this situation results in substantial overtreatment (5). Given that Gleason score is such a strong predictor of prognosis, identifying molecular markers that can augment the pathologic classification system could significantly improve risk stratification and treatment selection.
Using a limited gene set, we previously identified a 157-gene expression signature that distinguished high (Gleason score ≥8) from low Gleason (Gleason score ≤6) and significantly improved the prediction of lethal disease among men with Gleason score 7 (6). Since our original publication, this signature has been shown to predict lethal prostate cancer in additional gene expression studies (7–9). Here, in a new genome-wide gene expression dataset, we sought to improve upon this signature.
Patients and Methods
Study population
The men in the current study are participants from two ongoing prospective studies: the Physicians' Health Study (PHS; ref. 10) and the Health Professionals Follow-Up Study (HPFS; ref. 11). We included participants diagnosed with incident, histologically confirmed prostate cancer between 1982 and 2005, verified by systematic medical record review, through which we abstracted clinical data. Follow-up for metastases and mortality is available through March 1, 2011, in the PHS (>99% complete) and December 31, 2014, in the HPFS (>98% complete). The Human Subjects Committees at Partners HealthCare and the Harvard T.H. Chan School of Public Health approved this study.
Tissue collection and Gleason score assessment
The PHS and HPFS Tumor Cohort includes the subset of the men with prostate cancer from the cohorts for whom we have collected archival formalin-fixed paraffin-embedded (FFPE) radical prostatectomy (RP; 95%) and transurethral resection of the prostate (TURP; 5%) specimens (N = 2,200). A study pathologist (M. Fiorentino), blinded to clinical data, rereviewed all hematoxylin and eosin slides to assign Gleason scores, avoiding the problems of the shift in scoring over the past several decades (2, 12) and variability among pathologists (13, 14). After data generation, we returned to and further reviewed the Gleason 7 cases to attempt to identify from which pattern the core had been sampled. We could definitively identify the grade of the core sampled in 84 of the 241 Gleason score 7 cases.
Gene expression profiling
For a subset of the tumor cohort, we undertook whole-transcriptome gene expression profiling for lethal and indolent prostate cancer cases using the Affymetrix Human Gene 1.0 ST microarray. We sampled men using an extreme case design, including 113 men who died of their cancer or developed distant metastases (“lethal” cases) and 291 men who neither died of prostate cancer nor were diagnosed with metastases through 2011 and who lived at least 8 years after cancer diagnosis (“indolent” cases). For a subset (n = 202; 62 lethal and 140 indolent cases), we also profiled adjacent normal tissue.
RNA extraction and amplification and generation of expression profiles
Methods and procedures for RNA extraction and data generation are described elsewhere (15, 16). For the expression profiles generated, we regressed out technical variables and then shifted and normalized the residuals, as described previously (15). We mapped gene names to Affymetrix transcript cluster IDs using the NetAffx annotations as implemented in Bioconductor annotation package pd.hugene.1.0.st.v1, resulting in 20,254 genes. Gene expression data are available through Gene Expression Omnibus (GEO) accession number GSE79021.
Validation of the original 157-gene signature
The original 157-gene signature was developed using gene expression assessed on a custom-made Illumina array with approximately 6,000 genes (6). The signature uses the 157-gene expression values to calculate a continuous risk score between 0 and 1, with higher values, indicating expression patterns more similar to those of Gleason ≥8 tumors and lower values indicating more similarity to Gleason ≤6 tumors. It was built using a nearest shrunken centroids classifier [Prediction Analysis for Microarrays (PAM); ref. 17]. To assess this signature's performance in expression data generated with a different technology, we centered values of gene expression from each array and applied the classifier. Seventy-five PHS participants who were part of our original signature study are also in the current Affymetrix expression data, so we assessed how well the signature's risk scores for these men agreed across both arrays graphically and with Spearman correlation.
Generation of a new Gleason signature and assessment of predictive ability
To build a new signature, we implemented PAM and compared the 57 Gleason score ≤6 cases with the 106 Gleason score ≥8 cases. We restricted the analysis to the 7,880 genes demonstrating moderate to high gene expression (expression level ≥5) in at least one case (114 of the original signature's 157 genes were in this subset of expressed genes). We used 10-fold cross-validation to determine the optimal number of genes to distinguish the two Gleason score groups. We constructed box plots comparing the new signature genes' expression in Gleason ≤6 tumor tissue and Gleason ≥8 tumor tissue, as well as in all normal tissue (n = 202). For each signature gene, we ran t tests (P values calculated using the bootstrap to account for matching between some of the tumor and normal samples) to compare each pair of these groups with Bonferroni-adjusted P values <0.05 considered significant.
When applied to an individual, the signature uses the gene expression values for the 30 genes to calculate a continuous risk score between 0 and 1, with higher values indicating more similarity to Gleason ≥8 tumors and lower values indicating more similarity to Gleason ≤6 tumors. We applied the signature to the Gleason 7 cases and used the area under the receiver operating characteristic (ROC) curve (AUC) to assess how well the signature classified 4 + 3 versus 3 + 4 in these cases. We then assessed how the signature could improve prediction of lethality among men with Gleason 7 when used in conjunction with 3 + 4/4 + 3 status. To minimize bias in estimating this improvement, we restricted to the PHS to estimate the best weights for combining the signature with 3 + 4/4 + 3 status in a logistic regression model for lethality. We then applied this model to the HPFS and compared its ability to predict lethal outcomes with the ability of 3 + 4/4 + 3 status alone. We compared the ROC curves of these two models. As a sensitivity analysis, we repeated this assessment eliminating the participants with TURP tissue, because the TURP and RP procedures could potentially affect the RNA expression, and these procedures may more frequently capture tumors of different zonal origins with different patterns of gene expression (18, 19). We also compared this model with an analogous model built using the original signature.
We assessed the performance of a risk score constructed by adding together the gene expression values for genes positively associated with high score and subtracting the gene expression values for genes negatively associated with high score.
We created violin plots to visually compare the new signature scores among four categories of Gleason 7: Gleason 3 + 4/4 + 3 cross-classified with lethal/indolent outcome. All main analyses described above use the total Gleason score for each person; as a sensitivity analysis, we repeated these analyses among the 84 men where we could definitively identify the grade of the core.
Signature of Gleason score in normal tissue
We applied the new Gleason signature to the gene expression data from adjacent normal prostate tissue to determine whether the signature could differentiate normal tissue from Gleason score ≤6 men (n = 32) and from Gleason score ≥8 men (n = 47). Within the normal tissue, we also attempted to agnostically develop a signature that could separate these two groups using PAM.
Results
Clinical characteristics of the PHS and HPFS participants are presented in Table 1. Because we oversampled for lethal cases, the covariate distributions do not reflect those of the underlying population; for example, there are 14% Gleason score 5–6, 60% Gleason score 7, and 26% Gleason score 8–10.
Clinical characteristics of prostate cancer patients in the PHS and HPFS
Validation of the original 157-gene signature
The original signature's predictions of Gleason score calculated using the original data and the new data are plotted in Supplementary Fig. S1. The original signature showed good concordance (Spearman correlation 0.77; P < 0.0001) among the 75 PHS men with available gene expression data from the original and the current study, giving us confidence that the original signature could be applied to these new data. We then validated the original signature's performance in the HPFS: It successfully classified Gleason score ≤6 and Gleason score ≥8 [AUC = 0.86; 95% confidence interval (CI), 0.78–0.93; P < 0.0001; Supplementary Fig. S2] and, among Gleason score 7 cases, was statistically significant in a logistic regression model predicting lethality beyond 3 + 4/4 + 3 status alone (P = 0.02).
Development of new gene signature to predict Gleason score
We developed a new signature classifying Gleason score ≤6 versus Gleason score ≥8 tumors. The classifier was optimized at 30 genes, five of which were in the original signature (MYBPC1, SERPINA3, DPP4, MYLK, and ASPN). Twenty-nine of the 30 genes are downregulated in Gleason score ≥8 compared with Gleason score ≤6; only ASPN is upregulated in Gleason score ≥8. Box plots of gene expression of the 30-genes in Gleason score ≤6, Gleason score ≥8, and all available adjacent normal tissue are displayed in Fig. 1. The bootstrap P values from t tests comparing these three groups are presented in Supplementary Table S1, and the most significant (Bonferroni-adjusted P < 0.05) differences are indicated in Fig. 1. Twenty-three genes display similar expression values between normal and Gleason score ≤6, with differences only observed in Gleason score ≥8. Five genes (MYBPC1, SERPINA3, SOD2, NCAPD3, and ZFP36) were different in Gleason score ≤6 but similar in normal and Gleason score ≥8; whereas two genes (ANO7 and NR4A1) were different in all three tissues.
Expression levels of signature genes in Gleason score ≤6, Gleason score ≥8, and normal tissue. Expression levels of the 30 genes in the Gleason signature, in Gleason score ≤6 (Gl ≤6, n = 57) and Gleason score ≥8 (Gl ≥8, n = 106), and in all adjacent normal tissue (N, n = 202). The genes are in order of importance in the model (reading left to right, then down). Differences with Bonferroni-corrected bootstrap P value of ≤0.05 are indicated: significant decreases (left to right) are indicated with a “−”; significant increases (left to right) are indicated with a “+.” There are three main “patterns” of significant differences, indicated by the number of stars next to the gene name (0, 1, or 2). Genes without stars are similar between normal and Gl ≤6 but differ in Gl ≥8; genes with one star are significantly higher in Gl ≤6 compared with normal and significantly lower in Gl ≥8 when compared with both Gl ≤6 and normal; and genes with two stars are significantly higher in Gl ≤6 compared with normal and significantly lower in Gl ≥8 compared with Gl ≤6, but not significantly different between normal and Gl ≥8.
Prediction of lethal outcome
To assess whether the signature improves the prediction of lethal outcome among men with Gleason score 7, we fit a logistic regression model with the PHS participants (n = 76; 9 lethal and 67 indolent cases) to find appropriate weights for the signature and 3 + 4/4 + 3 status. On the basis of this model, a case with 4 + 3 has odds of being a lethal case that are 453% higher than a case with 3 + 4, and each 0.20 increase in the Gleason signature score is associated with a 33% increase in the odds of being a lethal case. When applied to the HPFS, the model of 3 + 4/4 + 3 status alone had an AUC for lethal prostate cancer of 0.68 (95% CI, 0.59–0.76). The new signature alone had an AUC for lethal prostate cancer of 0.73 (95% CI, 0.64–0.82). Although not a statistically significant difference (P = 0.23), this signature is a nominally stronger predictor of lethality than 3 + 4/4 + 3 status. The model including the signature and 3 + 4/4 + 3 status had an AUC for lethality of 0.76 (95% CI, 0.67–0.84); this was a statistically significant improvement over 3 + 4/4 + 3 status alone (P < 0.0001). ROC curves for this comparison are shown in Fig. 2. A sensitivity analysis repeating these steps after eliminating Gleason 7 cases with TURP tissue (3 in the PHS, 8 in the HPFS) showed a similar improvement when adding our signature to the model. In addition, the new signature showed a nonsignificant (P = 0.09) improvement over the original signature (AUC for the original signature with 3 + 4/4 + 3 = 0.72; 95% CI, 0.63–0.82).
ROC curves for the 30-gene signature. Curves are compared for the model including 3 + 4/4 + 3 status alone to the model, including both 3 + 4/4 + 3 status and the 30-gene signature (P < 0.0001).
Violin plots are presented in Fig. 3 to display the distributions of the signature predictions among four categories of Gleason score 7: Gleason 3 + 4 with indolent outcome, Gleason 3 + 4 with lethal outcome, Gleason 4 + 3 with indolent outcome, and Gleason 4 + 3 with lethal outcome. These plots illustrate that the signature is higher among the Gleason 4 + 3 cases—indeed, the signature is predictive of 4 + 3 versus 3 + 4 status (AUC = 0.74; 95% CI, 0.67–0.80). They further demonstrate the higher signature values for lethal cases compared with the indolent cases within each subcategory of Gleason score 7. All Gleason score 7 analyses were repeated among the 84 men for whom grade of the core sampled for expression profiling was known; results were similar and are presented in the Supplementary Materials and in Supplementary Fig. S3.
Violin plot of the values of the 30-gene signature. Predicted probability of being high score (y-axis) among four categories of Gleason 7: Gleason 3 + 4 with indolent outcome, Gleason 3 + 4 with lethal outcome, Gleason 4 + 3 with indolent outcome, and Gleason 4 + 3 with lethal outcome.
When we replaced the nearest shrunken centroid classification probability with an additive risk score based on the gene expression values, the performance of the signature was similar. The ROC curves for the nearest shrunken centroids signature and this summed risk score are overlaid in Supplementary Fig. S4, where we see that they are not producing identical results, but that their classification performance is similar. The prediction of lethality among Gleason score 7 cases using this simple summed risk score and 3 + 4/4 + 3 status had an AUC of 0.75 (95% CI, 0.68–0.81).
Signature of Gleason score in normal tissue
We applied the 30-gene signature of Gleason to the normal tissue expression data and found it did not distinguish the men whose tumor was Gleason score ≥8 from those whose tumor was Gleason score ≤6 (AUC = 0.49; 95% CI, 0.37–0.62), and we were not able to identify any gene signature that did. Expression data from normal tissue were, therefore, uninformative in predicting Gleason in the adjacent tumor.
Discussion
In whole-transcriptome expression profiling, we identified a new 30-gene signature that significantly improved predictions of lethal and indolent outcomes among intermediate-risk Gleason score 7 men beyond 3 + 4/4 + 3 status. We believe the predictive ability of this new signature, combined with its small size, makes it a potentially valuable clinical tool. The signature creates a continuous risk score, which can be added to pathologic grading to better inform treatment recommendations. Because the signature is predictive of worse prognosis whether or not 3 + 4/4 + 3 status is included in the model, it may be useful in patients where the Gleason score is uncertain or categorized differently by different pathologists. This signature will also be applicable to patients graded using the recent recategorization of differentiation proposed by the International Society of Urological Pathology, where it would improve the prediction of lethal disease among men in group II (Gleason score 3 + 4) and group III (Gleason score 4 + 3; ref. 20).
We began this study by validating our earlier 157-gene signature in a new study population. This earlier signature was developed using a limited targeted panel; we demonstrate that it can be applied to a new study population with gene expression measured on a different platform and successfully improve outcome prediction among Gleason score 7 cases. This confirmation demonstrates that our approach is robust and effective. By using the same methods with whole-transcriptome gene expression profiling, we identified a new 30-gene model that further improves prediction. The new signature's performance is similar whether we apply it directly as a nearest shrunken centroids classifier or by simply adding the gene expression values together in the direction of association, which further renders the signature easy to apply in new patients.
In the current study, we were not constrained by the targeted panel used previously, and this allowed us to restrict to genes that exhibited higher levels of expression in the tissue. Although genes with lower expression levels may be relevant biologically, they may not be as useful for prediction because of difficulties separating signal from noise and consistently measuring with different platforms. By building the signature with genes with higher expression levels, we are relying on genes more likely to be consistently expressed in prostate tissue, thus producing a more robust model.
A strength of this study is that we could evaluate the signature's performance in adjacent normal tissue. We hypothesized that our signature developed in tumor tissue might predict nearby Gleason when applied to normal tissue; however, we found that this was not possible, and that no signature was capable of doing this. That the signature was unable to distinguish Gleason score ≥8 from Gleason score ≤6 from normal tissue underscores how specific the signature is to the cellular and tissue architecture of tumor grade.
Our analysis examining the signature genes in normal, Gleason score ≤6, and Gleason score ≥8 tissue showed that the majority had similar gene expression values among normal and Gleason score ≤6. This suggests that these expression changes may be more associated with dedifferentiation than simply tumorigenesis. However, seven of the genes are significantly higher in Gleason score ≤6 than normal but are similar to normal or significantly lower than normal in Gleason score ≥8. These genes may play different roles during cancer development and progression.
The five genes selected in both the new and original signatures are compelling candidates for biological involvement with de-differentiation. MYBPC1 (myosin-binding protein c) plays a role in muscle contraction and was associated with Gleason score in the same direction (lower in high score) in a previous study (21). Another myosin-related gene was significantly lower in high score tumors (MYLK, myosin light chain kinase). This gene is differentially expressed in metastatic and nonmetastatic primary prostate cancer (22) and is downregulated in metastases (23). DPP4 has many roles including immune regulation, glucose metabolism, cell adhesion, and bone marrow mobilization. It has been shown to block fibroblast growth factor signaling in prostate tissue (24). In addition, variants within this gene influence PSA levels in serum (25), and the activity of this gene is reduced in the sera of metastatic prostate cancer patients (26). SERPINA3 (serpin peptidase inhibitor, clade A) is a plasma protease inhibitor that can inhibit molecules from immune cells and can form complexes with PSA in prostate tissue (27). Asporin (ASPN) was the one gene in the signature that had higher levels in high Gleason score compared with low Gleason score. ASPN prevents the ossification of periodontal ligaments by preventing the binding of bone-morphogenic protein 2 to its receptor, BMPR1B. Men who carry the ASPN germline D14 allele are at higher risk of developing prostate cancer metastases (28). This gene is highly expressed in cancer-associated fibroblasts (29), and it increases in the stroma of bone metastases in mice and primary prostate tumors in mice and humans (28, 30).
A clear theme across these genes is their possible role in prostate cancer stroma. Our colleagues have recently demonstrated in laser-capture microdissected prostate tissue that there are many more gene expression differences across Gleason score categories in tumor-associated stroma than in the epithelial cells (Tyekucheva and colleagues; unpublished data). This may be reflected in our signature; however, the gene differences could also be due to the varying proportion of epithelial and stromal cells across Gleason score because the tissue was macrodissected. Although microdissected tissue may more cleanly address this issue, a signature in macrodissected tissue has a more realistic application in a clinical setting. From an etiologic standpoint, these results add to the growing literature on the role that stromal genes play with respect to prostate cancer dedifferentiation and progression.
In the future, the genes identified in our signatures could potentially be targets to either inhibit dedifferentiation or maintain differentiation. Currently, our signatures can supplement pathologic score as a means of improving prediction of lethal disease. The next step will be to apply this signature in prostate cancer biopsy samples to assess its concordance across matched biopsy-prostatectomy samples and its performance predicting prognosis in diagnostic tissue. There is substantial current work trying to identify the best candidates for active surveillance, and patients diagnosed with biopsy Gleason 3 + 4 and other favorable risk factors may benefit from an active surveillance approach (31). Low-risk scores predicted by our 30-gene signature could serve as one such favorable risk factor and allow more men to avoid unnecessary treatment.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Disclaimer
The authors assume full responsibility for analyses and interpretation of these data.
Authors' Contributions
Conception and design: J.A. Sinnott, L.A. Mucci, M. Loda, K.L. Penney
Development of methodology: J.A. Sinnott
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Fiorentino, M.J. Stampfer, L.A. Mucci, M. Loda
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.A. Sinnott, S.F. Peisch, S. Tyekucheva, T. Gerke, M. Fiorentino, M.J. Stampfer, L.A. Mucci, M. Loda, K.L. Penney
Writing, review, and/or revision of the manuscript: J.A. Sinnott, S.F. Peisch, S. Tyekucheva, T. Gerke, R. Lis, J.R. Rider, M.J. Stampfer, L.A. Mucci, M. Loda, K.L. Penney
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.A. Sinnott, S.F. Peisch, R. Lis, K.L. Penney
Study supervision: K.L. Penney
Grant Support
The Physicians' Health Study was supported by grants CA34944, CA40360, CA097193, HL26490 and HL34595. The Health Professionals Follow-Up Study was supported by grants CA133891 and UM1CA167552. This study was supported by CA136578, CA141298, CA131945, and P50CA090381. J.A. Sinnott was supported by CA09001and the A. David Mazzone Career Development Award. J.P. Rider, L.A. Mucci, and K.L. Penney were supported by Prostate Cancer Foundation Young Investigator Awards.
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.
Acknowledgments
We are grateful to the participants and staff of the Physicians' Health Study and Health Professionals Follow-Up Study for their valuable contributions, as well as the following state cancer registries for their help: AL, AZ, AR, CA, CO, CT, DE, FL, GA, ID, IL, IN, IA, KY, LA, ME, MD, MA, MI, NE, NH, NJ, NY, NC, ND, OH, OK, OR, PA, RI, SC, TN, TX, VA, WA, and WY.
Footnotes
Note: Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).
- Received May 17, 2016.
- Revision received July 27, 2016.
- Accepted September 15, 2016.
- ©2016 American Association for Cancer Research.
References
Volume 23, Issue 1
