The population-based Swedish Watchful Waiting Cohort consists of 1256 men with localized prostate cancer. These men had symptoms of benign prostatic hyperplasia (lower urinary tract symptoms) and were subsequently diagnosed with prostate cancer. All men in this study were determined at the time of diagnosis to have clinical stage T1–T2, Mx, N0, according to the 2002 American Joint Commission Committee Tumor-Node-Metastasis staging system as previously described (6
). The prospective follow-up time is now up to 30 years. The regional cohort includes men who were diagnosed at University Hospital in Örebro (1977–1991) (8
) and at four centers in the southeast region of Sweden: Kalmar, Norrköping, Linköping, and Jonköping (1987–1999) (). All patients with prostate cancer were recruited through an informed consent process at the respective institutions. This study is compliant with Karolinska and Örebro Ethical Committees. A subset of men from these cohorts (n = 388) were included in the study. Inclusion criteria required the availability of greater than 90% tumor cells compared with surrounding stroma or benign tissue in the diagnostic Trans Urethral Radical Prostatectomy (TURP) biopsy sample. Samples included were derived from an equal proportion of men who died of prostate cancer or developed metastasis and men who lived a minimum of 10 years without clinical recurrence of their disease. Of these 388 patients, only the 354 with reliable TMPRSS2-ERG fusion results were included in the analyses.
Characteristics of the Swedish Watchful Waiting Cohort (with available tissue blocks, N = 1256) by Center
Physician Health Study (PHS) Prostatectomy Confirmation Cohort
This cohort included 116 US men who were diagnosed with prostate cancer between 1983 and 2003, and were treated by radical prostatectomy as primary therapy (11
) (). The men were participants in an ongoing randomized trial in the primary prevention of cancer and cardiovascular disease. This study was approved by the Harvard School of Public Health Institutional Review Board, and all patients provided written informed consent at time of initial enrollment. Only the 101 patients with reliable TMPRSS2-ERG fusion results were included in the analysis.
Clinical characteristics of Physician Health Study prostate cancer cohort individuals included in the study, 1983–2003*
cDNA-mediated Annealing, Selection, Ligation, and Extension Array Design
We designed a set of four cDNA-mediated annealing, selection, ligation, and extension (DASL) Assay Panels (DAPs) for the discovery of molecular signatures relevant to prostate cancer (Hoshida Y, Setlur SR, Perner S, Camargo A, Gupta S, Moore J, Reich M, Gabriel S, Rubin MA, Golub TR, in preparation). We prioritized informative genes, i.e., genes showing differential expression across samples in previously generated microarray data sets (the datasets are at http://www.broad.mit.edu/cancer/pub/HCC
), which included 24 studies, 2149 samples, and 15 tissue types. The top-ranked transcriptionally informative genes that showed the largest variation in expression across the different datasets comprised genes in most of the known biological pathways. To ensure that prostate cancer–related genes were included in the DAP, we performed a meta-analysis of previous microarray datasets from the Oncomine Database (12
) and included from that a list of genes that were transcriptionally regulated in prostate cancer. The final array consisted of 6144 genes (6K DAP).
Sample Processing and cDNA-mediated Annealing, Selection, Ligation, and Extension
Foci highly enriched for prostate cancer (>90%) were identified by microscopic examination of the tissue sections by the study pathologists (MAR, SP). Three 0.6-μm biopsy cores per patient were taken from these enriched areas and were placed in one well of a 96-well plate for high-throughput RNA extraction. The CyBi-Well liquid handling system (CyBio AG, Jenna, Germany) was used for high-throughput extraction. Cores were first deparaffinized by incubation with 800 μL Citrisolv (Fisher Scientific, USA) at 60°C for 20 minutes and then with 1.2 mL Citrisolv:absolute alcohol (2:1) at room temperature for 10 minutes. Cores were then washed with absolute alcohol, dried at 55°C, and incubated overnight at 45°C in 300 μL lysis buffer (10 mM NaCl, 500 mM Tris pH 7.6, 20 mM EDTA, 1% SDS) containing 1 mg/mL proteinase K (Ambion, Austin, TX). RNA was extracted from the lysate using the TRIzol LS reagent (Invitrogen, Carlsbad, CA). TRIzol LS reagent (900 μL) was added to the cell lysate, followed by 240 μL of chloroform (Sigma-Aldrich, St. Louis, MO). The samples were mixed thoroughly and centrifuged at 4°C, 5600g for 40 minutes (the same centrifugation settings were used for the rest of the protocol). After centrifugation, the aqueous phase was transferred to a new plate, and the RNA was precipitated by incubation with 620 μL of isopropanol (Sigma-Aldrich) at room temperature for 10 minutes. Glycogen (20 μg; Invitrogen) was added as a carrier. The samples were centrifuged as above, and the pellet was washed with 80% ethanol (Sigma-Aldrich), air dried, and dissolved in RNase-free water. The RNA was quantified using a NanoDrop spectrophotometer (NanoDrop technologies, Wilmington, DE).
SYBR green (QIAGEN Inc., Valencia, CA) quantitative polymerase chain reaction (qPCR) assay for a housekeeping gene, ribosomal protein L13a (RPL13A), was used to estimate RNA quality (RNA with crossover threshold, Ct, of less than 30 cycles was considered to be good quality). Primer sequences for RPL13A were as follows: RPL13A-FWD, GTACGCTGTGAAGGCATCAA, and RPL13A-REV, GTTGGTGTTCATCCGCTT (GenBank accession NM_012423.2). DASL expression assay (Illumina Inc., San Diego, CA) was performed using 50 ng of cDNA according to manufacturer’s instructions.
Cell Lines and Transfection
The prostate cancer cell lines NCI-H660, VCaP, PC3, DU145, and 22Rv1 were obtained from American Tissue Culture Collection (ATCC, Manassas, VA). Cells were maintained according to the supplier’s instructions.
VCaP cells were transiently transfected with an ERβ-containing plasmid (kindly provided by M. Lupien) using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instructions. Transfection medium was removed after 6 hours, cells were washed in phosphate-buffered saline (PBS, 137 mM NaCl, 10 mM phosphate, 2.7 mM KCl, pH of 7.4) twice, and phenol red–free DMEM (Cellgro Mediatech, Herndon, VA) supplemented with 5% charcoal/dextran-treated-fetal bovine serum (CDT-FBS) (Invitrogen) was added. ERβ mRNA expression levels were determined after transfection by qPCR using the following primers: ERβ-FWD, AAGAAGATTCCCGGCTTTGT and ERβ-REV, TCTACGCATTTCCCCTCATC (GenBank accession code, NM_001437.2).
NCI-H660 cells were transiently transfected with smart pool small-interfering RNA (siRNA) against ERβ or anti-LUC control siRNA (both from Dharmacon Inc., Chicago, IL), both at a concentration of 25 nM, using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instructions.
17β-Estradiol (E2, Sigma-Aldrich), the ERα agonist propylpyrazole triol (PPT, Tocris Bioscience, Ellisville, MO), and the ERβ agonist diarylpropionitrile (DPN, Tocris Bioscience, Ellisville, MO) were each dissolved in 100% ethanol. Raloxifene, tamoxifen, and fulvestrant (Sigma-Aldrich) were each dissolved in dimethyl sulfoxide (DMSO). All reagents were used at a final concentration of 10nM.
NCI-H660 and VCaP cells were hormone deprived by culture in their respective phenol red–free media (for NCI-H660 without E2 and hydrocortisone) supplemented with 5% CDT-FBS (Invitrogen), for 48 hours (VCaP) or for 72 hours (NCI-H660). Transfected cells were treated with hormones or vehicle 24 hours after transfection.
NCI-H660 cells and transfected VCaP-ERβ cells were treated with the following compounds: E2, DPN, PPT, raloxifene, fulvestrant, or tamoxifen; all at 10 nM final concentration; or vehicle for 12, 24, or 48 hours. Untransfected VCaP cells were treated for 12 or 24 hours with E2, DPN, raloxifene, or fulvestrant (all 10 nM final concentration).
Results were analyzed using GraphPad Prism version 4.00 for Windows (GraphPad Software, San Diego CA).
Determination of TMPRSS2-ERG Fusion Status in Biopsy Samples and Expression in Hormone-Treated Prostate Cancer Cell Lines
TMPRSS2-ERG fusion status was determined by ERG break-apart fluorescence in situ hybridization (FISH) assay (13
)(n=362) and qPCR (for cases not assessable by FISH (n=98). An aliquot of the RNA used for DASL was used for qPCR. cDNA was synthesized as above using the Illumina kit (Illumina Inc., San Diego, CA). The TMPRSS2-ERG fusion product was detected using SYBR green assay (QIAGEN) with TMPRSS2-ERG_f and TMPRSS2-ERG_r primers (GenBank accession code NM_DQ204772.1) (3
). RPL13A was used for normalization. RNA from NCI-H660 cells, which express TMPRSS2-ERG (14
), was used as a positive control and a calibrator for quantification. Relative quantification was carried out using the comparative ΔΔCt
The same protocol was used to quantify the TMPRSS2-ERG fusion product after treatment of NCI-H660 and VCaP prostate cancer cell lines with estrogenic and antiestrogenic compounds. In this case, cDNA was synthesized using Omniscript RT kit (QIAGEN Inc), and a housekeeping gene, hydoxymethylbilane synthase (HMBS) was used for normalization. The primer sequences for HMBS are as follows: HMBS-FWD, CCATCATCCTGGCAACAGCT and HMBS-REV, GCATTCCTCAGGGTGCAGG (Gen Bank accession code NM_000190.3). Two independent experiments were performed, with each sample in triplicate.
Cell Viability Assays
NCI-H660 cells were seeded in 96-well plates (approximately 5 × 103 cells per well) and treated for 8 days with hormone (E2, PPT, or DPN) or vehicle alone as above. Relative cell number was determined before (time 0 used for normalization) and 2, 3, 6, and 8 days after treatment using the Cell Titer-Glo luminescent assay (Promega Corporation, Madison, WI) according to the manufacturer’s instructions. Two independent experiments were performed, both in octuplicate.
Expression of Estrogen Receptor α and Estrogen Receptor β in Prostate Cancer Cell Lines
Quantitative Polymerase Chain Reaction
Total RNA was extracted from the VCaP, NCI-H660, LNCaP, PC3, DU145, and 22Rv1 cell lines. The RNA was reverse transcribed (RT), and 50 μg of the resultant cDNAs was subjected to PCR analysis. RT-PCR was carried out using primers for ERα(16
) (GenBank accession code, NM_000125.2) and ERβ (17
) (GenBank accession code, NM_001437.2). cDNA was synthesized using the Omniscript RT kit (QIAGEN Inc.), and HMBS was used for normalization. Two independent experiments were performed with each sample in triplicate.
Expression of ERαand ERβ in was assessed in NCI-H660 cells and in untransfected VCaP cells and VCaP-ERβ cells (48 hours after transfection). Whole-cell extracts were prepared in RIPA buffer (50 mM Tris pH 7.5, 150 mM NaCl, 2 mM sodium orthovanadate, 0.1% Nonidet P-40, 0.1% Tween 20) with 1x Complete Protease Inhibitor Cocktail (Roche, Indianapolis, IN). Protein concentration was determined using the Bio-Rad DC protein assay (Bio-Rad Laboratories, Hercules, CA). Equal amounts (20 μg) of total protein were loaded on NuPAGE 4–12% Tris-Bis gels (Invitrogen) and transferred to Immobilon-P polyvinylidene fluoride membranes (Millipore, Billerica, MA). Blots were incubated with primary antibodies (mouse monoclonal anti-ERα[1:100, NeoMarkers, Labvision Corporation, Fremont, CA] or mouse monoclonal anti-ERβ [1:200, clone 14C8, GeneTex Inc., San Antonio, TX]), washed three times with PBS containing 0.1% Triton X-100, and then incubated with peroxidase-conjugated anti-mouse secondary antibody (1:8000, Amersham Biosciences, Piscataway, NJ) for 1 hour. Rabbit polyclonal anti-β-Actin (1:1000, Cell Signaling technology, Danvers, MA) was used as a control for protein loading and transfer. Antibody–protein complexes were detected using the ECL Western Blotting Analysis System (Amersham Biosciences, Piscataway, NJ) according to the manufacturer’s instructions. Three independent experiments were carried out.
ChIP was performed as described by Carroll et al. (18
). Briefly, NCI-H660 cells were hormone -deprived by culture for 3 days in phenol red–free medium (Cellgro Mediatech, Herndon, VA) supplemented with 5% CDT-FBS (Invitrogen) lacking E2 and hydrocortisone. Cells were treated with E2 or vehicle for 60 min, and chromatin was crosslinked using 1% formaldehyde (15
). Due to the high homology between ERαand ERβ, sites analyzed included ERαrecruitment sites that had previously been identified upstream of the TMPRSS2
gene (ER3429–ER3433) in an unbiased genome-wide ChIP-Chip experiment in the MCF-7 breast cancer cell line (18
). In addition, the TMPRSS2 promoter (TMPRSS2_prom) was added to this analysis. Crosslinked chromatin was immunoprecipitated with mouse monoclonal anti-ERβ antibody (1:200, clone 14C8, GeneTex Inc.). The precipitated DNA was amplified using primers spanning the putative ER binding sites (Supplementary Tables 2 and 3
, available online). Three independent experiments were performed.
We used several strategies to both identify and evaluate a gene signature of TMPRSS2-ERG prostate cancer. Briefly, we identified candidate genes by t-test statistic via repeated sampling to ensure robustness. We then built different classification models to evaluate prediction performance of the gene signature using both the Swedish and PHS cohorts. Specifically, we first tested the method on the Swedish cohort by means of a hold-out procedure, i.e., two-thirds (n=235) of the samples were used to build the model and one-third (n=119) was used to evaluate it. In each evaluation phase, a different subset of samples from those used to build the models was used, thus ensuring that the final model did not depend on a specific dataset. Second, a model built using the entire Swedish dataset was used to evaluate the PHS dataset.
Gene Selection and Classification Models Built on the Swedish Cohort
A feature selection procedure (identification of candidate genes) was carried out by applying a two-sided t-test on each gene on 10 repeats of 10-fold cross-validation, resulting in 100 t-tests for each gene. A 10-fold cross-validation works as follows: within each iteration, 1/10 of the samples are held out as “test” and a t-test is computed for each gene using the remaining 9/10 of the samples, known as “training” samples. This procedure, called partitioning, is performed to ensure that similar proportions of fusion-positive and -negative cancers are present in the two partitions. Partitioning is then repeated nine times, resulting in 10 disjoint “test” partitions. This approach ensures that each sample is used as a “training” and a “test” set at least once. In addition, to prevent results that are biased to a specific random splitting of the data, the 10-fold cross-validation procedure was repeated another nine times.
Genes were then selected if they were found to be statistically significantly associated with the TMPRSS2-ERG fusion as determined by FISH or qPCR in at least 50% of the iterations. The number of genes selected was not set a priori but rather depended on the data. This procedure was first applied within the Swedish cohort training set (defined as two-thirds of the entire Swedish cohort) using a P value threshold of .0001. The selected genes were then used to build the classification model (see Classification Model) to predict fusion status based on their expression. The resulting classification model was evaluated on the Swedish cohort validation set (one-third of the entire Swedish cohort) to obtain an initial assessment of classifier performance.
To verify whether the molecular signature seen in the Swedish cohort was present in an independent set of cancers (different population), we used the PHS cohort as an additional evaluation dataset. For this procedure, the same feature-selection (i.e selection of genes) was applied within the entire Swedish cohort, using a P value threshold of .00001. The classification model built using this gene signature on the entire Swedish cohort was then evaluated on the independent set of cancers from the PHS cohort. Finally, to determine the sensitivity of the molecular signature model with respect to the number of genes selected, we used the above described iterative procedure to select the top ranking 75%, 50%, 25%, and 10% of genes from the gene signature (ranking is defined by the frequency with which a subset of genes are called “statistically significant” during the iterative procedure) and evaluated the classifier performances on the PHS dataset.
Several classification models based on gene expression values were generated. Support vector machines (SVMs) (19
) using both radial basis function and polynomial kernels (degrees equal to 1, 2, or 3) with different costs (cost equal to 0.01, 0.1, 1, or 10) were used. The area under the receiving operator characteristic (ROC) curve (AUC) was used as a performance measure for SVM classification models. The 10 repeats of 10-fold cross-validation previously described were also used to select the best SVM parameters. Each iteration included genes with a P
value of .00001or less that had been selected by two-sided t-test. The genes selected here were used only for the purpose of selecting the best SVM model. The best SVM parameters were identified as the ones giving highest mean AUC computed on the test sets. The mean AUC was computed on the total number of iterations, namely 100.
The AUC P value was evaluated via a randomization approach, specifically by counting how many times, out of 1000 iterations, randomly obtained class predictions outperformed the actual classification model. A binomial distribution was used to generate random predictions by setting the probability of positive class according to the frequency of the SVM-predicted positive cancers. Finally, the 95% confidence interval [CI] of the AUC was estimated from the 10 repeats of 10-fold cross-validation within the Swedish training set.
Gene Lists Comparison
To further assess the robustness of the gene signature, we compared the gene list obtained with the iterative procedure described in the previous paragraphs with that obtained by means of a standard two-sided t-test controlled with false discovery rate (FDR, Q value<.01) (20
) within the Swedish training set (n=235). In addition, we performed the iterative gene selection on the PHS cohort and measured the overlap with the genes identified in the entire Swedish cohort. A hypergeometric distribution test was computed to assess for the statistical significance of the overlap.
We selected the gene lists based on different thresholds of t-test P
values and on the percentage of genes that were statistically significant throughout the repeated cross-validations. We obtained similar results as those presented in this manuscript. Analysis was performed using R 2.4.0 (21
Computational Functional Analysis
The gene signature identified with the entire Swedish cohort was used for connectivity map (CMAP) (22
) and molecular concepts map (MCM) (23
) analyses. Connectivity map analyzes the association (ie, the positive or negative correlation) between the given test signature and gene expression profiles of cell lines treated with specific concentrations of various drugs. MCM analyzes the association between the given test signature and various gene sets or “molecular concepts”. Fisher exact two-sided test was used for pairwise comparison of the concepts (MCM, P
<.005, odds ratio >1.5).