A Tissue Biomarker Panel Predicting Systemic Progression after PSA Recurrence Post-Definitive Prostate Cancer Therapy

Background Many men develop a rising PSA after initial therapy for prostate cancer. While some of these men will develop a local or metastatic recurrence that warrants further therapy, others will have no evidence of disease progression. We hypothesized that an expression biomarker panel can predict which men with a rising PSA would benefit from further therapy. Methodology/Principal Findings A case-control design was used to test the association of gene expression with outcome. Systemic (SYS) progression cases were men post-prostatectomy who developed systemic progression within 5 years after PSA recurrence. PSA progression controls were matched men post-prostatectomy with PSA recurrence but no evidence of clinical progression within 5 years. Using expression arrays optimized for paraffin-embedded tissue RNA, 1021 cancer-related genes were evaluated–including 570 genes implicated in prostate cancer progression. Genes from 8 previously reported marker panels were included. A systemic progression model containing 17 genes was developed. This model generated an AUC of 0.88 (95% CI: 0.84–0.92). Similar AUCs were generated using 3 previously reported panels. In secondary analyses, the model predicted the endpoints of prostate cancer death (in SYS cases) and systemic progression beyond 5 years (in PSA controls) with hazard ratios 2.5 and 4.7, respectively (log-rank p-values of 0.0007 and 0.0005). Genes mapped to 8q24 were significantly enriched in the model. Conclusions/Significance Specific gene expression patterns are significantly associated with systemic progression after PSA recurrence. The measurement of gene expression pattern may be useful for determining which men may benefit from additional therapy after PSA recurrence.


Introduction
The majority of men with prostate cancer are now diagnosed with cancers that have a low risk of cause-specific mortality [1]. These men are usually treated with radical retropubic prostatectomy (RRP), external beam radiotherapy, or interstitial brachytherapy and are then followed by regular serum PSA evaluations. Over the next 5 to 10 year period, 15-30% of these men will develop a rising PSA [2][3][4][5][6], defining a rapidly growing population of major clinical and public health significance. Of this PSA relapse group some men will have local recurrence or have clinically-detectable metastasis, but many will have no other evidence of recurrent prostate cancer other than the rising PSA. The PSA ''doubling time'' has been identified as a potential surrogate for cause-specific mortality, and is used by some clinicians to determine which men with PSA relapse deserve adjuvant hormonal ablation, local radiation therapy, or simple observation [4][5][6]. Biomarkers that predict which of these men would benefit from any additional therapy are needed.
Large scale gene expression studies of prostate cancers of different grade and stage have been performed by several groups [7][8][9][10][11][12][13][14][15][16]. These expression studies have utilized arrays containing probe sets of up to 35,000 genes. While these studies are important for biomarker discovery, several difficulties preclude their translation into a clincial setting. First, it is likely that smaller panels will be used clinically. Second, because the previous studies required frozen material, the number of specimens analyzed was limited. Third, since adverse clinical events in prostate cancer patients require lengthy followup, the testing methods must be applicable to archival paraffinembedded material. Finally, none of the previous studies was focused on the development of a biomarker panel to predict prostate cancer systemic progression in the setting of PSA recurrence.
Using the Mayo Clinic Radical Retropubic Prostatectomy (RRP) Registry, we designed a nested case-control study to test the hypothesis that a limited set of RNA expression biomarkers can predict which men with a rising PSA post-RRP might benefit from additional clinical intervention. The Illumina DASL TM expression microarray platform was selected as the biomarker measurement method, because it measures the expression of gene targets using paraffin tissues [17][18][19]. Using expression data from the literature and derived from our own research program we developed a limited set of expression markers that would likely be altered in association with prostate cancer progression. The panel also included expression biomarkers from several other previously published panels that are associated with surrogates (high Gleason Score, high pathologic stage, or metastatic disease) for prostate cancer aggressiveness [12][13][14][15][16].
We report that the array-based measurements showed excellent correlation with quantitative RT-PCR measurements of paraffinderived RNAs. We also report excellent intra-array, inter-array and within-gene reproducibility. We then describe the testing and validation of a gene expression tissue biomarker panel for the prediction of prostate cancer systemic progression following a rising PSA after radical prostatectomy. We compare the performance of our panel with other previously published panels. Finally, we show that the overexpression of genes mapped to chromosome band 8q24 is associated with prostate cancer systemic progression.

Gene Selection and Array Design for the DASL TM Assay
Two Illumina DASL expression microarrays were utilized for the experiments. The standard commercially available Illumina DASL expression microarray (Cancer Panel TM v1) containing 502 oncogenes, tumor suppressor genes and genes in their associated pathways. Seventy-eight of the targets on the commercial array have been associated with prostate cancer progression.
A custom Illumina DASL TM expression microarray containing 526 gene targets for RNAs, including genes whose expression is altered in association with prostate cancer progression. Probes for the custom DASLH panel were designed and synthesized by Illumina, Inc. (San Diego, CA).
Four different sets of prostate cancer aggressiveness genes were included in the study (if the genes were not present on the Cancer Panel v1 array, they were included in the design of the custom array):

1)
Markers of prostate cancer aggressiveness identified by a Mayo/University of Minnesota Partnership [20]: The expression profiles of 100 laser-capture microdissected prostate cancer lesions and matched normal and BPH control lesions were analyzed using Affymetrix HG-U133 Plus 2.0 microarrays. Ranked lists of significantly over-and under-expressed genes comparing 10 Gleason 5 and 7 metastatic lesions to 31 Gleason 3 cancer lesions were generated. The top 500 genes on this list were compared to lists generated from prior expression microarray studies and other marker studies of prostate cancer (see 2-4 next). After this analysis there was space for 204 novel targets with potential association with aggressive prostate cancer on the custom array.

3)
Previously published markers associated with prostate cancer aggressiveness (e.g. PSMA, PSCA, Cav-1): Expression microarray data has also been published. This literature was evaluated for additional tissue biomarkers. For example, at the time of array design we were able to identify 13 high quality expression microarray studies of prostate cancer aggressiveness (See Tables S1 and S2 for full reference list). In addition, among the 13 reports, 5 papers presented 8 expression biomarker panels to predict prostate cancer aggressiveness [12][13][14][15][16]. When appropriate probes suitable for the DASL chemistry could be designed for these panels they were included on the custom array. We also identified 12 articles reviewing genes associated with prostate cancer. These criteria resulted in the selection of 150 genes.

4)
Markers derived from Mayo SPORE research (including genes and ESTs mapped to 8q24). Ninety-three additional biomarkers were identified (see Tables S1 and S2).
The custom array also included probe sets for 45 genes that were not expected to differ between case and control groups based on Mayo/University of Minnesota Partnership data. Thirty-eight of these genes were also present on the commercial array (see Tables S1 and S2).
After enumerating the potentially prostate cancer relevant genes on the commercially available cancer panel 570 potentially prostate cancer relevant genes and 451 other cancer-related genes were evaluated across both arrays.

Design of Nested Case-Control Study
For this study we sampled individuals from the Mayo Clinic RRP Registry. The registry consists of a population of men who received prostatectomy as their first treatment for prostate cancer at the Mayo Clinic (For a current description and use of the registry; see reference [23]). As systemic progression is relatively infrequent, we designed a case-control study nested within a cohort of men with a rising PSA. Between 1987-2001, inclusive, 9,989 previously-untreated men had RRP at Mayo. On follow-up, 2,131 developed a rising PSA (.30 days after RRP) in the absence of concurrent clinical recurrence. PSA rise was defined as a followup PSA . = 0.20 ng/ml, with the next PSA at least 0.05 ng/ml higher or the initiation of treatment for PSA recurrence (for patients whose follow-up PSA was high enough to warrant treatment). This group of 2,131 men comprises the underlying cohort from which SYS cases and PSA controls were selected.
Within 5 years of PSA rise, 213 men developed systemic progression (SYS cases), defined as a positive bone scan or CT scan. Of these, 100 men succumbed to a prostate cancer-specific death, 37 died from other causes and 76 remain at risk.
PSA recurrence controls (213) were selected from those men without systemic progression within 5 years after the PSA rise and were matched (1:1) on birth year, calendar year of PSA rise and initial diagnostic pathologic Gleason score (, = 6, 7+). Twenty of these men developed systemic progression greater than 5 years after initial PSA rise and 9 succumbed to a prostate cancer-specific death.
A set of 213 No Evidence of Disease (NED) Progression controls were also selected from the Mayo Clinic RRP Registry of 9,989 men and used for some comparisons. These controls had RRP from 1987-1998 with no evidence of PSA rise within 7 years of RRP. The median (25 th , 75 th percentile) follow-up from RRP was 11.3 (9.3, 13.8) years. The NED controls were matched to the systemic progression cases on birth-year, calendar year of RRP and initial diagnostic Gleason Score. Computerized optimal matching was performed to minimize the total ''distance'' between cases and controls in terms of the sum of the absolute difference in the matching factors [24].
The current study was approved by the Institutional Review Board of Mayo Clinic.

Block Identification, RNA Isolation, and Expression Analysis
The list of 639 cases and controls was randomized. An attempt was made to identify all available blocks (including apparently normal and abnormal lymph nodes) from the randomized list of 639 eligible cases and controls. Maintaining the randomization, each available block was assessed for tissue content by pathology review and the block containing the dominant Gleason pattern cancer was selected for RNA isolation.
Four freshly cut 10mm sections of FFPE tissue were deparaffinized and the Gleason dominant cancer focus was macrodissected. RNA was extracted using the High Pure RNA Paraffin Kit from Roche (Indianapolis, IN). RNA was quantified using ND-1000 spectrophotometer from NanoDrop Technologies (Wilmington, DE). The RNAs, including intra-plate and inter-plate replicates, were distributed on 96-well plates in the randomized order for DASL analysis.
RNA samples were processed, hybridized to Sentrix Universal 96-Arrays, scanned using BeadArray Reader, and data initially processed in BeadStudio according to the manufacturer's instructions. Microarray data is available from the GEO database (http:// www.ncbi.nlm.nih.gov/geo/ accession number GSE10645).
To evaluate the accuracy of the gene expression levels defined by the DASL technology, we performed quantitative SYBR Green RT-PCR reactions for 9 selected ''target'' genes (CDH1, MUC1, VEGF, IGFBP3, ERG, TPD52, YWHAZ, FAM13C1, and PAGE4) and 4 commonly-used endogeneous control genes (GAPDH, B2M, PPIA and RPL13a) in 384-well plates, with the use of Prism 7900HT instruments (Applied Biosystems, Foster City, CA). 210 RNA samples with abundant RNA from the group of total 639 patients were analyzed. Because of RNA shortage, only 77 samples were analyzed for PAGE4. mRNA was reversetranscribed with SuperScript III First Strand Synthesis SuperMix (Invitrogen, Carlsbad, CA) using random hexamers. For each of the nine genes studied, the cycle threshold (Ct) was determined in triplicate and the expression was normalized relative to the set of four reference genes.

Pathology Review
The Gleason score in the Mayo Clinic RRP Registry was defined as the initial Gleason score. Since there have been changes in pathologic interpretation of the Gleason score over time, a single pathologist (JCC) reviewed the Gleason score of each of the blocks selected for expression analysis. This clinical variable was defined as the revised Gleason score.

Statistical Methodology
Collection of gene expression data was attempted for the 623 patients as described in Results. Of these, there were 596 (n SYS = 200, n PSA = 201, n NED = 195) patients for whom data was collected, the rest having failed one or both expression panels as described in Results. To assure selection of similar training and validation sets, 100 case-control-control cohorts comprised of 133 randomly chosen SYS patients (two-thirds of 200 for training) along with their matched PSA and NED controls were selected as a proposed training set. The remaining cases and controls were treated as a proposed validation set. The clinical variables were tested for independence between the proposed training and validation sets separately within the SYS cases and the PSA controls. Discrete clinical factors (pathologic stage, hormonal treatment adjuvant to RRP, radiation treatment adjuvant to RRP, hormonal treatment adjuvant to PSA recurrence, and radiation therapy adjuvant to PSA recurrence) were tested using Chi-square analysis. Continuous clinical variables (Gleason score (revised), age at PSA recurrence, first rising PSA value, second rising PSA value, and PSA slope) were tested using Wilcoxon rank sum. Six of the one hundred randomly sampled sets failed to show dependency for any of the clinical variables at the 0.2 level, and the first of these was chosen as the training set: 391 patients (n SYS = 133, n PSA = 133, n NED = 125). This reserved 205 patients for the validation set (n SYS = 67, n PSA = 68, n NED = 70).
The raw data from BeadStudio was normalized using cyclic loess (fastlo) [25]. The training data were analyzed using random forests [26] using R Version 2.3.1 (http://www.rproject.org) and randomForest version 4.5-16 (http://stat-www. berkeley.edu/users/breiman/RandomForests). The data were analyzed by panel (Cancer, Custom and Merged, where Merged was the Cancer and Custom data treated as a single array). By testing the ntree parameter of the randomForest function we determined that 4000 random forests were sufficient to generate a stable list of markers. The top markers as sorted for significance by the randomForest program were combined with various combinations of clinical variables using logistic regression R program (glm() with family = binary (a logistic model), where glm refers to generalized linear model). The resulting scoring function was then analyzed using Receiver Operating Characteristic (ROC) methods and the cut-off was chosen that assumed an equal penalty for false positives and false negatives. A review of the models permitted a subset of markers to be identified, and a subset of supporting clinical data identified. The number of features in the model was determined by leave 1/3 out Monte Carlo Cross Validation (MCCV) using 100 iterations. The number of features was selected to maximize AUC and minimize random variation in the model. The final model was then applied to the 391 patient training set and the reserved 205 patient validation set. For comparison, other previously reported gene expression models were also tested against the training and validation sets [12][13][14][15][16].
We compared the previously reported models for their classification of patients into the known PSA recurrence control and SYS progression case groups. We used the Cramér's Vstatistic [27] to compare models.

Study Design/Paraffin Block Recovery/RNA Isolation and Expression Panel Success
Briefly, a nested case-control study was performed using the large, well-defined cohort of men with rising PSA following RRP( Figure S1). SYS cases were 213 men who developed systemic progression between 90 days and 5.0 years following the PSA rise. PSA controls were a random sample of 213 men who were 5 years post-RRP with PSA recurrence but with no evidence of further clinical progression. NED controls were a random sample of 213 men who were 7 years post-RRP without PSA rise (the comparison of PSA controls with NED controls-will be presented in a subsequent paper). SYS cases and PSA controls were matched (1:1) for birth year, calendar year of PSA rise, and initial pathologic Gleason score (, = 6 vs. . = 7). The list of eligible cases and controls was randomizeed for the blind ascertainment of blocks, isolation of RNA and performance of the expression array experiments. Table 1 summarizes the distribution of clinical parameters between the SYS cases and the PSA and NED control groups. There was no significant difference between the groups for the matching variables (there was no significant difference in initial diagnostic Gleason score when the , = 6 and .7 groups-the matching criteria-were compared). Comparison of the initial diagnostic Gleason score to the revised Gleason scores revealed that Gleason scores have increased over time. In addition, the proportion of Gleason 8-10 tumors increased comparing NED controls to PSA controls, and PSA controls to SYS cases. The revised Gleason score was used in all the biomarker modeling.
All paraffin-embedded blocks from eligible men were identified and each block was surveyed for the tissue present (primary and secondary Gleason cancer regions, normal and metastatic lymph nodes, etc.). We macrodissected the dominant Gleason pattern region and attempted to isolate RNA. Illumina Cancer Panel TM and custom prostate cancer panel DASL array analyses were then performed on all RNA specimens. The Methods section and Tables S1 & S2 describe the composition of the Cancer Panel and the design of the custom panel. Table 2 summarizes the final block availability, the RNA isolation success rate and the success rates of the expression array analyses. Of the 639 eligible patients, blocks were available on 623 (97.5%). RNA was isolated and DASL assays successfully performed on a high proportion of patients/specimens: usable RNA was prepared from all 623 blocks, and the Cancer Panel and custom panel DASL arrays were both successful (after repeating some specimens-see below) on 596 RNA specimens (95.7% of RNAs; 93.3% of design patients). Only 9 (1.4%) RNA specimens failed both panels. The primary reason for these failures was poor RNA quality-as measured by qRT-PCR of the RPL13a gene expression [19]. Of the 1246 initial samples run on both panels, 87 (7.0%) specimens failed. Those specimens for which there was residual RNA were repeated with a success rate of 77.2% (61 of 79 samples).

Expression Analysis Reproducibility
Replicate analysis results ( Figure S2), RT-PCR comparisons ( Figure S3) and inter-and intra-panel gene expression comparisons are described in Results S1.

Specific Gene Expression Results Comparing the Systemic Progression Cohorts with the PSA Recurrence and No Evidence of Progression Cohorts
Univariate Analyses by gene. Because the DASL assay appeared to generate precise and reproducible results, the array data was examined for genes whose expression was significantly altered when the SYS cases were compared with the PSA controls. For this initial analysis, the DASL gene expression value was determined to be the average of up-to-three probes for each gene on each array. Upon univariate analysis (two-tail t-test) of the probe-averaged and fastlo normalized data [25], 68 genes were highly significantly over-or under-expressed in the SYS cases versus PSA controls (p,9.73610 27 , Bonferroni correction for p,0.001) ( Table 3). One hundred twenty-six genes were significantly over-or under-expressed in the SYS cases versus the PSA controls (p,4.86610 25 , Bonferroni correction for p,0.05). Table S3 provides the complete gene list ordered by pvalue. Figure 1 illustrates nine genes with significantly different expression in the SYS cases and PSA controls.

Systemic Progression Prediction Model Development and
Testing on a Training set. A statistical model to predict systemic progression (with and without clinical variables) using a training set was developed using random forests [21] and logistic regression as described in Methods. The training data were analyzed by panel (cancer, custom and merged), by gene (the average expression for all gene-specific probes), and by individual probes. Table 4 lists the 15 genes and 2 individual probes selected for the final model. Table 5 and Figure 2A summarize the areas under the curve (AUCs) for three clinical models, the final 17 gene/probe model and the combined clinical probe models. The variables in the clinical models (Table 6) were based on available clinical information. Clinical model A included revised Gleason score and pathologic stage (information available immediately after RRP). The addition of diagnostic PSA and age at surgery did not significantly add to the AUC and was left out of this model (data not shown). Clinical model B added age at surgery, preoperative PSA value, and any adjuvant or hormonal therapy within 90 days after RRP (information available after RRP but before PSA recurrence). Clinical model C added age at PSA recurrence, the second PSA level at time of PSA recurrence, and the PSA slope (information available at the time of PSA recurrence).
Using The arrays were selected to include probe sets for several previously published prostate aggressiveness models [12][13][14][15][16]. Table 5 summarizes the AUCs for array expression results for these biomarker models. Figure 2C illustrates the AUCs for four of these models with the appropriate comparison with clinical model C and with the 17 gene/probe model. Each of these models generated AUCs that were smaller than the model we developed. However several of the models generated AUCs (e.g. We also compared the models for their classification of patients into the known PSA recurrence control and SYS progression case groups. Table S4 summarizes    model, both with clinical data. Most of the models classified the same patients into the known groups (e.g. classifying a patient in the PSA control group as a PSA recurrence and a patient in the SYS case group as a systemic progression). They also tended to incorrectly classify the same patients (e.g classifying a patient in the PSA control group as a systemic progression and vice versa). The 17 gene/probe model correctly classified 5-15 more patients into their known category (PSA controls or SYS cases) compared to the other models (data not shown).

Secondary Analyses
Exploratory Survival studies. As noted above, the 17 gene/ probe model and the previously reported models each classified some of the SYS cases in the good outcome category (e.g. to be PSA recurrences, not systemic progressors) and some of the PSA controls in the poor outcome category (e.g. to go on to systemic progression). We were curious to know whether these apparently false classifications had any biologic or clinical relevance.
Seventeen men in the PSA control group (who had both array and clinical model C data) went on to have systemic progression beyond 5 years at the time of last follow-up. Of these 17 patients, 9 were predicted to have a poor outcome by the 17 gene/probe model. Of the 179 patients who did not have any systemic progression, 38 were classified in the poor outcome category by the model (p value = 0.0066, Fisher exact test). Figure 3A illustrates the systemic progression-free survival for the good and poor outcome groups in the PSA controls. PSA controls with a tumor classified as having a poor outcome had significantly increased risk for developing systemic progression beyond 5 years (log rank p-value = 0.00050) (HR = 4.7, 95% CI: 1.8-12.1). Ninety-three men in the SYS case group (who also had array and clinical model C data) went on to prostate cancer death at the time of last follow-up. Of these 93 patients, 78 were predicted to have a poor outcome by the 17 gene/probe model. Of the 98 patients who did not suffer a prostate cancer death, 61 were classified in the poor outcome category by the model (p value = 0.0008, chi-square test). Figure 3B illustrates the prostate cancer-specific overall survival for the good and poor outcome groups in the SYS cases. SYS cases with a tumor classified as having a poor outcome had significantly increased risk for suffering a prostate cancer-specific death (HR = 2.5, 95% CI:  Figure 3C illustrates the prostate cancer-specific overall survival for the good and poor outcome groups in the SYS cases. SYS cases whose tumor classified as having a poor outcome had significantly increased hazard of suffering a prostate cancerspecific death (HR = 2.3, 95% CI: 1.3-4.2). The median survival from first positive bone scan or CT was 3.1 years (95% CI: 2.5-4.3) in the group classified as having a poor outcome and 8.6 years (95% CI: 8.3-') in the group classified as having a good outcome (log rank p-value = 0.0033). Gini Decrease for a variable is the average (over all random forest trees) decrease in node impurities from recursive partitioning splits on that variable. For classification, the node impurity is measured by the Gini index. The Gini index is the weighted average of the impurity in each branch, with impurity being the proportion of incorrectly classified samples in that branch. The larger the Gini decrease, the fewer the misclassification impurities. ** Genes mapped to 8q24 *** Single probes for C8orf53 and CDK10 were selected. The Mean Gini Decrease for these probes are derived from an independent random forest analysis of the all probes separately. doi:10.1371/journal.pone.0002318.t004 Table 5   Exploratory 8q24 Studies. Because of recent tumor chromosome dosage and germ line association studies, the custom array included 82 8q genes on the custom array. Fourteen 8q genes were within the top 68 genes based upon univariate analysis (Table 3). Compared to the proportion of 8q genes on both arrays the prevalence of 8q genes is non random (p = 0.003, Fisher exact test). Twelve additional 8q genes were within the top 126 genes. The prevalence of 26 8q genes in the top 126 is statistically significant (p = 1.56610 25 , Fisher exact test). Chromosome band 8q24.1 has the greatest overrepresentation of genes in the top 68 gene and 126 gene lists (11 genes, p = 6.35610 27 and 19 genes, p = 9.34610 212 , Fisher exact test). Of the 17 genes/probes in our final model, five map to 8q24 (p = 0.0043, Fisher exact test)(see Table 4).
Exploratory ets Transcription Factor Studies. Alterations of several ets-family oncogenes are associated with the development of prostate cancer [28][29][30]. We included oligonucleotide probe sets for the three major members of the ets family involved in prostate cancer: ERG, ETV1, and ETV4, as well as their translocation partner TMPRSS2. Figure 4 summarizes the expression results for these genes for the SYS cases and the PSA and NED controls. Several observations can be made: 1) With only 8 exceptions ERG, ETV1 and ETV4 overexpression are mutually exclusive; i.e. the overexpression of each generally occurs in different tumors. 2) Different probe sets for ERG give nearly identical expression results (see Figure S4A). 3) The prevalence of ERG overexpression was 50.0%, 52.2% and 53.8% in the SYS cases, PSA controls and NED controls, respectively. There was no significant difference in the mean expression and the prevalence of ERG overexpression between the three cohorts (see Figure 4). 4) The prevalence of ETV1 overexpression was 11.5%, 6.5% and 5.1% in the SYS cases, PSA controls and NED controls, respectively (see Figure 4). The prevalence of ETV1 overexpression was significantly higher in SYS cases (p = 0.043, chi-square test). 5) The prevalance of ETV4 overexpression ranged from 2.5%-5.5% among the three groups and was not significantly different. 6) None of the genes were selected by the formal statistical modeling (see Table 4). In fact, the 17 gene/ probe model predicted similar rates of progression in ERG+ and ERG-patients (data not shown).
Exploratory Pathway Analysis. We used the 461 genes from both cancer and custom panels that are potentially differentially expressed between SYS cases and PSA controls (p#0.05) as the focus genes for Ingenuity Pathway Analysis (IPA, Ingenuity Systems Inc. Redwood City, CA). IPA identified 101 canonical pathways that are associated with the focus genes, 51 of which are overrepresented with p#0.05 (see Table S5). However, because we measured a limited number of genes on both DASL panels, the p values from IPA analysis may not accurately quantify the degree of overrepresentation of focus genes in each pathway.
We then performed Gene Set Enrichment Analysis (GSEA) [31], on chromosome 8 genes grouped by map location. Genes mapped to 8q24.1 had a significant p value (p = 0.0002) with a FDR q value = 0.001 (see Table S6).

Discussion
Patients with a rising PSA following definitive therapy comprise a heterogeneous cohort; a significant number develop metastasis, followed by hormone refractory prostate cancer. Of these, a substantial number, but not all, will die of the disease. PSA failure following RRP or radiation therapy is associated with a 15% to 25% five year prostate cancer death rate [32,33]. Androgen deprivation has been increasingly used in all stages of prostate cancer to improve mortality rates [34,35] or to facilitate prostate cytoreduction [36,37]. Two recent studies described the natural history of progression after PSA elevation following RRP or radiation therapy [32,38]. They identified PSA doubling time as a potential surrogate for prostate mortality. In three retrospective studies early androgen deprivation in patients with biochemical failure and short (,12 months) PSA doubling time after prostatectomy improved survival [39][40][41]. We hypothesized that additional biomarkers beyond PSA doubling time could help predict which men with a rising PSA post-RRP might suffer systemic progression. Such a panel could be incorporated into future prospective clinical trials in the setting of PSA progression.
Using an array methodology optimized for RNA from paraffinembedded tissues and a rigorous statistical modeling algorithm, we developed a 17 gene/probe tissue gene expression model to predict the likelihood of systemic progression in men with a rising PSA post-RRP. In a training set the 17 gene/probe model was significantly better than the use of clinical variables alone. While accuracy decreased when the 17 gene/probe model was further tested with a reserved validation set, the performance of the 17 gene/probe model with clinical model C was better than the clinical model alone.
The reduction in AUC between the training and validation sets was in part due to the overfitting inherent in these types of analyses. Since we maximized the AUC on the training set, the validation set AUC would be predicted to be lower. Another cause of the reduction in AUC could be a relative lack of precision of the Illumina DASL technology. Except for the poor correlation of the DASL and RT-PCR measurements for genes with low DCt values, all of the intra-plate, inter-plate, intra-gene and inter-gene reproducibility analyses suggested that the DASL chemistry was very precise. We observed greater coefficients of variation in our replicate RT-PCR measurements than in the DASL measurements (data not shown).
Perhaps the best explanation for the validation set AUC reduction is that prostate cancer is genetically heterogeneous. This heterogeneity can result in a reduction in a validation set AUC even for relatively large datasets. This hypothesis predicts that several different models could be developed from the same dataset. One of the advantages of the Illumina DASL platform is its ability to analyze up to 1536 probes (or 512 genes if thee probes are selected per gene) on a single array. We included probes from eight previously reported gene expression panels associated with prostate cancer aggressiveness [12][13][14][15][16]. The models showed strong correlation with each other and generally predicted the same patients to be PSA recurrences or systemic progressors. A recent comparison of several breast cancer gene expression models by Fan et al. [42] also showed high correlation between models. The implication is that several independent biomarker panels can be developed to classify cancer patients with similar clinical or biologic endpoints. After the formal analysis was completed we began secondary exploratory studies of the whole dataset. We hypothesized that systemic progression beyond 5 years in the PSA controls might be predicted by the models. When the 17 gene/probe model score (and three of the other models) predicted a poorer outcome, there was a high likelihood that a PSA control would have a positive bone scan or CT beyond 5 years. There was also evidence supporting a second hypothesis that prostate cancer-specific death in the SYS cases might be predicted by the models. When the 17 gene/probe model score (and three of the other models) predicted a poorer outcome, the median overall survival of a SYS case from a positive bone scan or CT was 2.8 years (compared to 8.6 years with a better model score). These secondary analyses imply that tissue expression biomarker panels may have utility for the stratification of patients for interventions at the time of PSA recurrence as well as for systemic progression. Importantly, the expression data was collected on primary tumor specimens resected several years before the occurrence of the clinical events.
Overrepresentation of 8q24 is associated with clinically aggressive prostate cancer (for example, references [43][44][45]). Furthermore, tumor overexpression of genes on chromosome 8 (and from 8q24) is also reproducibly associated with prostate cancer progression [46][47][48][49]. We recently mapped the region of 8q24 overrepresentation, and it involves a ,5 Mb region surrounding c-Myc [50]. Our secondary analyses demonstrated that 8q24 genes were significantly overrepresented in the top 68 and 126 genes by t-test and in the final 17 gene/probe model. Each of the genes exhibited similar magnitudes of overexpression in the SYS cases suggesting an association with chromosomal dosage. A common germline polymorphism mapped near the c-Myc gene on 8q24, has recently been associated with prostate cancer development [51]. This finding has been replicated by at least five different groups, with the further suggestion that at least two germline haplotypes near c-Myc are associated with prostate cancer development [52][53][54][55][56]. It is not known if the men who inherit the at-risk haplotype(s) have a clinically more aggressive prostate cancer, a poorer prognosis, somatic (tumor) overrepresentation of 8q24, or overexpression of 8q24 genes.
Alterations of several ets-family oncogenes are associated with the development of prostate cancer [28][29][30]. Our panel(s) included probe sets for three members of the ets family involved in prostate cancer; ERG, ETV1, and ETV4, as well as their translocation partner, TMPRSS2. As a group, these genes are over-expressed in approximately 62% of prostate cancers. This overall prevalence was nearly identical in our three case and control groups. In addition, with the possible exception of ETV1, whose prevalence of overexpression was about 2-fold higher in the SYS cases, none of the genes seemed to be associated with systemic progression of prostate cancer. It has been recently reported that ERG fusion is associated with lethal prostate cancer in Scandinavian men treated with watchful waiting [30]. However, the prevalence of fusion (and presumably ERG overexpression) in that study was only 15%; far lower than in our dataset and other reports [28,29]. These differences are likely a result of the types of prostate cancer (and clinical outcomes) diagnosed where PSA screening is common (North America) and uncommon (Scandinavia).
We conclude that the measurement of gene expression patterns may be useful for determining which men are likely to benefit from additional therapy following PSA recurrence. These measurements should be included in prospective evaluation of various therapeutic interventions when PSA rises following definitive treatment of prostate cancer. Figure S1 Summary of the nested case-control study design.