Multivariate analysis.

We used a forward stepwise selection process to define an optimal model for prediction of prostate cancer outcome, yielding a four-variable model containing TOP2A, CDH10, and the fusion and aneuploidy status variables. At each iteration, all remaining variables were tested for their additional predictive power given the variables selected thus far (Appendix Tables A1, A2, and A3, online only). The first variable was TOP2A. With TOP2A in the model, BIRC5 and RRM2 were not significant predictors. The predictive performance of the model did not change significantly if either BIRC5 or RRM2 were selected as the first variable of the model. When any one of TOP2A, BIRC5, and RRM2 was selected as the first variable in the model, the other two genes lost significance, demonstrating these genes are interchangeable and suggesting they are related. The next most significant variable added to the model was CDH10. After adding CDH10 and TOP2A to the model, the significance of the fusion status increased and became the most significant remaining parameter (see Appendix for additional information). Adding fusion status to TOP2A and CHD10 rendered the remaining genes insignificant. Each successive step of this modeling improved the corresponding area under the curve (AUC) in the ROC plot by at least 0.02 (Fig 1). The AUC with TOP2A alone was 0.71— it increased to 0.74 with CDH10, and increased to 0.79 with the predicted fusion parameter (Fig 1). Aneuploidy remained a significant variable in this analysis, so it was added to the final model (Table 3) consisting of predicted fusion, TOP2A, CHD10, and aneuploidy which resulted to an AUC of 0.81 (95% CI, 0.74 to 0.88). The change in AUCs at each step was also tested for significance using the integrated discrimination improvement (IDI) method proposed by Pensina (Appendix) which agreed with the stepwise result. Further confirmation for the significance of the final model is provided by permutation tests and by the penalized fit to the original data, which also showed only minor changes in the coefficients. The score derived from the regression model is displayed in Figure 1B, showing the overall separation of the cases and the controls.

Independent validation using frozen tumor tissue from a separate cohort of 57 high-risk patients showed a slight decrease in the AUC to 0.79 (95% CI, 0.66 to 0.92; Fig 2A). The score derived from the regression model in this independent validation is displayed in Figure 2B.

It is important to note that patients were not matched for adjuvant postoperative therapy (hormone and radiation); however, nearly equal numbers of cases and controls received adjuvant therapy. In addition, stratified analysis indicated no association of adjuvant therapy with the final model, the treatment modalities were not significant predictors of outcome and the coefficients and P values of the other variables in the model changed only slightly.

Model building case-control design.

To assess the potential utility of selected biomarkers and to facilitate the construction of a multivariable model on the prediction of systemic progression, a case-control study design was used. Men with systemic progression within 5 years of RRP were identified, and matched with controls known to be free of systemic progression at 7 years. Subjects were also matched on GS, TNM stage, margin status, and preoperative serum prostate-specific antigen (PSA). Systemic progression was defined as the development of metastatic disease as determined by clinical, radiologic (bone or computed tomography scan), or pathologic (biopsy) evaluation. PSA progression alone was not sufficient for declaration of a systemic progression. Patients undergoing RRP at the Mayo Clinic are evaluated postoperatively at least quarterly for 1 year, semiannually for 1 year, and annually thereafter. Digital rectal examination and serum PSA are evaluated at each visit. If patients had an abnormal elevation in serum PSA postoperatively, a radioactive bone scan, plain radiograph in the presence of an abnormal bone scan, and/or computed tomography were performed. Patients that did not return to the Mayo Clinic were mailed kits for blood submission and serum PSA testing, and additional medical information was obtained from the local physicians as needed (Blute ML, et al: J Urol 165:119-125, 2001).

Based on this data, an initial set of 200 samples, comprising matched cases and controls, was defined. Tissues were acquired from the tissue bank and reviewed by two pathologists blinded to the case-control status. Of the 200 samples, 76 cases and 81 controls passed the pathology review and had sufficient tissue of the highest Gleason pattern cancer available for experimental analysis. The clinical and pathologic features of the selected subjects are listed in Table A2. A number of cases and controls were not included due to an absence of tumor block or insufficient tumor for analysis. In addition, some controls were excluded because of inaccurate Gleason grading; often the block containing the highest grade of the tumor was exhausted by previous studies. This occurred for both controls and cases. The controls were selected based on the surgical pathology report, and all materials were reviewed by a urologic pathologist to ensure appropriate Gleason grading when matching cases and controls. All GS6 samples were excluded from this study, hence cases and controls were GS ≥7. Due to the exclusion of several subjects following review by the pathologists, exact case-control matching was not preserved. However, the retained cases and controls remained balanced on age, preoperative serum PSA, TNM stage, GS, margin status, and adjuvant (< 90 days after surgery) therapy, and an unmatched analysis was used. Based on the final numbers the study had approximately 80% power to detect (α = .05, two sided) a mean difference in gene expression between cases and controls equivalent to 0.45 standard deviations. For present/absent gene expression, a 22-percentage point difference (eg, 40% v 62%) between cases and controls could be detected.

Validation case-control design.

An independent nested case-control design was conducted to assess the model identified in the model building component of the analysis. Twenty-one systemic progression cases were identified from 2000 to 2005 inclusive who had both fresh frozen and paraffin-embedded tissue. (Fresh frozen specimens have been archived since 1999.) For each case, two potential controls were identified from those who were at risk but who had not yet progressed at the time of the case's event (time of systemic progression or death) and matched on GS, pathologic stage, age, and date of surgery. The clinical and pathologic features of the selected subjects are listed in Table A3. The total number of controls with sufficient RNA was 36. Because of the risk-set sampling design the average follow-up of controls in this set was much less than the first case-control set of patients. Nevertheless, the risk-set sampling design was necessary as there is not yet sufficient follow-up in this cohort for a 7 year rule. We would expect that the consequences of this design are that it would still validate the model from stage 1, but that the model would discriminate slightly less well.

Processing of FFPE samples.

In all experiments, the processing of samples was randomized to prevent processing biases. Each case was reviewed by a pathologist and tumor was identified on H-E stained sections. Subsequent sections (10 μm) from each case were prepared under RNase free condition and de-paraffinized with Xylene. Identified tumor areas were scraped into 2 mL tubes containing digestion buffer (RecoverAll kit; Ambion, Austin, TX). Total RNA was isolated according to the RecoverAll RNA isolation procedure. The isolated RNA was treated with DNase using Ambion Turbo DNA free kit, according to manufacturer's instructions (Ambion). The amount of RNA in each case was measured by the Quant-iT RiboGreen kit (Invitrogen Carlsbad, CA). Reverse transcription (RT) was performed using Superscript III First Strand Synthesis system (Invitrogen) and 500 ng of RNA from each case in a 40 μL reaction volume.

Processing of frozen samples for validation.

Tissue was cut by our tissue-processing core facility and was kept frozen. Sections were stored on slides at −80 for less than 1 week. Slides were placed directly from −80% to 75% ethanol (ETOH) for 30 seconds × 2 to remove OTC, then placed in 95% ETOH for 30 seconds and 100% ETOH for 30 seconds and air dried. Tissue was scraped directly into lysis solution from the RNeasy mini or midi kit (Qiagen, Valencia, CA) and processed immediately according to manufacturer's instructions. Reverse transcription and quantitative polymerase chain reaction (qPCR) were performed as described for FFPE tissue except 200 ng of RNA was used in the RT reaction and 1 μL cDNA was used in subsequent qPCR (5 ng RNA equivalent).

qPCR was performed on each sample by adding 12.5-ng total RNA equivalent of cDNA to a 20 μL reaction volume for each gene of interest using SYBR green PCR Master Mix (Applied Biosystems ABI, Foster City, CA) on the ABI 7900HT real-time PCR machine using the manufacturer's default cycling conditions. Primers for qPCR were designed using Primer Express software (ABI) to amplify a 70 to 85 base pair fragment from the Affymetrix target sequence for the gene of interest. The primer pair concentrations (0.15 or 0.2 μmol/L final) were optimized by generating standard curves using a pool of prostate cDNA from normal and tumor tissue. To check for genomic DNA, No RT samples were run in a qPCR reaction and those with cycle threshold (Ct) less than 35 for GapDH were considered contaminated with DNA and were reprocessed. In the analysis of data, undetermined values for Cts were replaced with a Ct of 40, which was the maximum cycle number in the qPCR experiments.

All samples were analyzed in duplicate and all studies were carried out under Mayo institutional review board–approved protocols. Primers used for the quantitative RT-PCR expression analysis of the genes in the final model are: Forward reverse: GAPDH, CATGGCCTCCAAGGAGTAAGAC TCTCTTCCTCTTGTGCTCTTGCT; RPS 28, GCTGCTCGCTGGGTCTTG GGAGCAGATTGTGACAGACCATT; ERG, GCTGCCACAATCAGAAATCA TCGCGACTCAAAGGAAAACT; TOP2A, TGGCTGCCTCTGAGTCTGAA AGTCTTCTGCAATCCAGTCCTCTT; CDH10, GAACAGGATAGTTCTCCCTTAAGCA CAAGGGCAGGACATGTACCTAAC; ETV1, TGTTTTTGCTTTGCATTTGG TCCCCATTTACTCATGGTTTTT; ETV4, GCAGATCCCCACTGTCCTAC CCACTTTTCCTTCCCAATGA.

Data analysis.

All quantitative PCR measurements were normalized by subtracting the number of cycles measured for a candidate gene from the average number of cycles measured for GAPDH and RPS28 in the same tissue sample. This normalization method inverts the amplification values, such that higher values correspond to higher expression. The full set of un-normalized data was used to assess the variance of the assay using a Bland-Altman plot of the average of a pair of replicates versus their difference (not shown). Assay variability was close to constant for values below 33, from which point it increased linearly. The inverse of the estimated variance was used to create a weighted average for each replicate pair. In samples where one replicate measure returned no value, the other replicate value was used. In samples where both replicate measures returned no value, the measurement was treated as missing.

Ploidy analysis.

Flow cytometry was performed as previously described (Zanetta G, et al: Am J Obstet Gynecol 175:1217-1225, 1996). The nuclear content of 10,000 nuclei was measured with a FACScan (Becton Dickinson, Sunnyvale, CA) flow cytometer. Cell cycle evaluation of the DNA histogram was performed with a Modfit 5.2 (Verity Software, Topsham, ME) computerized software program. Tumors with only one identifiable gap0-gap1 peak were classified as DNA diploid (2n). Tumor samples that contained a significant increase in the 4n peak (> 9% of nuclei) and an identifiable 8n population were categorized as DNA tetraploid. Tumor DNA content was classified as DNA aneuploid if a separate, identifiable gap0-gap1 peak was present. All DNA histograms were analyzed and classified without knowledge of the clinicopathologic features or patient outcome.

Role of ploidy.

In this case-control study, cases and controls were not matched for aneuploidy as this is not a standard pathologic assessment of prostate cancer in RRP specimens, and its association with prostate cancer outcome is debated. Ploidy generally correlates with other prognostic factors, such as GS, tumor stage, and tumor volume: low-stage tumors are usually diploid, and high-stage tumors are usually nondiploid, and therefore some investigators have shown that there is no association of tumor ploidy with outcome in a multivariate model (Humphrey PA, et al: Am J Surg Pathol 15:1165-1170, 1991). However, in a Mayo Clinic study, when patients with tumor that had spread beyond the prostate were examined, tumor ploidy was found to be significantly associated with progression (Robertson CN, et al: Acta Oncol 30:205-207, 1991; Winkler HZ, et al: Mayo Clin Proc 63:103-112, 1988). In the Mayo Clinic practice, cancer ploidy is assessed on every RRP specimen, and therefore we investigated the significance of ploidy in the molecular model. The ploidy status parameter is often defined as 0 for diploid tumors, 1 for tetraploid tumors, and 2 for anauploid. However, we decided that it would be difficult to separate tetraploid from proliferation so we defined a binary variable only for aneuploidy (1 for aneuploid and 0 for everything else). This variable is still a significant predictor in univariate and multivariate analysis and entered into the model after the fusion-status variable. As measurement of aneuploidy is not a standard of care for prostate cancer we decided to report the model with and without the aneuploidy in this manuscript (Fig 1 in the manuscript).

Fusion status analysis.

The presence of fusions between TMPRSS2 and ERG, ETV1, or ETV4 was assigned based the expression values of ERG, ETV1, and ETV4 genes. The thresholds for status assignments were determined for each gene by analyzing the normalized Ct expression values of these genes. The fusion-status parameter was set to 1 when one of the ETS factors (ERG, ETV1, or ETV4) was above their specified threshold and to −1 when all of them were below the thresholds. If the RT-PCR failed for one of these three genes or the values were to close to call that fusion-status was set to 0. The thresholds were set by plotting histograms of the expression values for all cases for each of these genes and finding two cutoff points (Fig A2). Samples where the expression was greater than the upper threshold in any of the three genes (ERG > −4.8; ETV1 > −5.5; ETV4 > −6) were designated fusion positive. Samples where the expression was less than the lower threshold for all three genes (ERG < −5.4; ETV1 < −6.8; ETV4 < −5.9) were assigned fusion negative. The same procedure was implemented for the frozen samples and the thresholds were −5, −3, and −7, respectively.

Methods.

In the cohort study, the Cox proportional hazards model was used to evaluate predictors of time to systemic progression. The concordance statistic (C), which has an interpretation similar to receiver operator characteristic curve area was used to compare models. In the case-control study the association of gene expression with systemic progression was assessed using univariate and multivariate logistic regression with the dependent variable (0 = control, 1 = case) and the normalized genes as the independent predictors. Computations were done using the S-Plus statistical package (Insightful Corp, Seattle WA) and the open-source R statistical package. The clustering analysis was performed using clustering methods in R.

Multivariate analysis.

The P values for the univariate and the stepwise logistic regression analysis can be followed in Table A1. In iteration 1 the model was adjusted for TOP2A (model 1). In iteration 2 the model was adjusted for TOP2A and CDH10 (model 2). In iteration 3 the model was adjusted for TOP2A, CDH10, and pred_fusion (model 3). The final model included TOP2A, CDH10, and pred_fusion and aneuploidy (model 4).

The change in area under the curve from model 1 (more than 2), model 2 (more than 3), and model 3 (more than 4), was tested for significance using the integrated discrimination improvement statistic proposed by Pensina et al (Pensina MJ, et al: Statist Med 27:157-172, 2008). The P values were .006, .001, and .025 respectively for each step. We also applied bias reduced regression on the final model using the plr function in R. The penalized fit showed only minor changes in the coefficients (shrinkage was chosen to give 3.5 effective df).

Importance of fusions in the multivariate analysis.

The ETS-fusion status became a strong predictor of outcome after CDH10 entered in the model (Table A1). Logistic regression showed a strong association of CDH10 expression with the fusion status in all three data sets. As we examined the association of the ETS fusion with CDH10, we observed that CDH10 expression was variably overexpressed in much higher levels in ETS fusion–positive samples than in ETS-fusion–negative samples. The pred_fusion status variable entered the model to compensate for this association. At the time we do not know the biologic explanation for the association of CDH10 and ETS transcription factors. We are in the process of performing studies to better understand this association. Of more clinical significance we observed that if we apply our model in a fusion dependent manner, we see an improvement in the model's predictive value in fusion-positive cases to 0.82 area under the curve. Our data suggests that the presence or absence of the fusion is associated with unique gene expression profiles of prognostic markers, and further studies are underway. For example, COL2A1 and KHDRBS3 genes, that have a strong association with the fusion, also became strong predictors when the CDH10 entered into the model presumably imposing a similar adjustment.

Independent validation.

The concordance index (area under the curve) for the validation set is 0.79 (0.66,0.92) when using the coefficients derived from the training set. We also attempted to fit the four variables in this set as well. The concordance index (area under the curve) of the panel fit for the validation set was 0.82 (95% CI, 0.7 to 0.94). The model parameters are as follows: variable original model coefficient (95% CI) refit coefficient: TOP2A, 0.71 (95% CI, 0.38 to 1.04) 0.81; CDH10, 0.62 (95% CI, 0.93 to 0.31) 0.78; pred_fusion, 0.91 (95% CI, 0.36 to 1.46) 1.07; aneuploidy 1.34 (95% CI, 0.14 to 2.54).

As shown earlier, the model coefficients for the individual variables were also consistent in the validation set with the exception of aneuploidy. In the validation data set only three cases were aneuploid, hence the coefficient is not estimatable. If there was serious overfitting in the training set we would expect that the coefficients in the validation set to be shrunk toward 0. Further confirmation is provided by the penalized fit to the original data, which also showed only minor changes in the coefficients (shrinkage was chosen to give 3.5 effective df).