The Mutational Landscape of Lethal Castrate Resistant Prostate Cancer

Characterization of the prostate cancer transcriptome and genome has identified chromosomal rearrangements and copy number gains/losses, including ETS gene fusions, PTEN loss and androgen receptor (AR) amplification, that drive prostate cancer development and progression to lethal, metastatic castrate resistant prostate cancer (CRPC)1. As less is known about the role of mutations2–4, here we sequenced the exomes of 50 lethal, heavily-pretreated metastatic CRPCs obtained at rapid autopsy (including three different foci from the same patient) and 11 treatment naïve, high-grade localized prostate cancers. We identified low overall mutation rates even in heavily treated CRPC (2.00/Mb) and confirmed the monoclonal origin of lethal CRPC. Integrating exome copy number analysis identified disruptions of CHD1, which define a subtype of ETS fusionnegative prostate cancer. Similarly, we demonstrate that ETS2, which is deleted in ~1/3 of CRPCs (commonly through TMPRSS2:ERG fusions), is also deregulated through mutation. Further, we identified recurrent mutations in multiple chromatin/histone modifying genes, including MLL2 (mutated in 8.6% of prostate cancers), and demonstrate interaction of the MLL complex with AR, which is required for AR-mediated signaling. We also identified novel recurrent mutations in the AR collaborating factor FOXA1, which is mutated in 5 of 147 (3.4%) prostate cancers (both untreated localized prostate cancer and CRPC), and showed that mutated FOXA1 represses androgen signaling and increases tumour growth. Proteins that physically interact with AR, such as the ERG gene fusion product, FOXA1, MLL2, UTX, and ASXL1 were found to be mutated in CRPC. In summary, we describe the mutational landscape of a heavily treated metastatic cancer, identify novel mechanisms of AR signaling deregulated in prostate cancer, and prioritize candidates for future study.


Integration of exome sequencing with transcriptome sequencing of prostate cancer cell lines
As transcriptome sequencing has also been used to discover recurrent mutations in cancer 17,18 , we analyzed the transcriptome of 11 prostate cancer cell lines (primarily CRPC, Supplementary Table 11 Table 12). Given the lack of normal genomic DNA from these cell lines, germline and somatic variants cannot be distinguished. Thus, we only considered variants fulfilling one of three high stringency filters as likely somatic mutations: 1) deleterious variants affecting a gene harboring a somatic mutation in our study (Supplementary Table 13), 2) variants affecting the same nucleotide as a somatic mutation in our study (Supplementary Table   14 Table 15).
This integrative approach identified additional variants in TP53, AR and APC, supporting the utility of the analysis. A TP53 R248W variant, present in WA10 and previously reported as somatic 19 , was identified in the VCaP cell line, while previously reported P223L and V274F somatic variants were identified in DU-145 16 , with a V274G variant present in WA37. A confirmed somatic TP53 variant R175H was identified in both WA30 and LAPC-4, consistent with previous reports (Supplementary Table 15) 19 . Finally, a Y234H confirmed somatic variant (predicted to be probably damaging) was also present in C4-2B (Supplementary Table 15).
Importantly, integrating transcriptome sequencing data also identified recurrent variants in genes not previously identified as being mutated in prostate cancer, including STAG2, MLL3, CNOT1, FAM123B (WTX) and FOXA1 (Supplementary Tables 13-15). WA32 harbored a R370W somatic mutation in STAG2, and a R370G variant was identified by transcriptome sequencing in LNCaP; mutations in STAG2 have recently been identified as causing aneuploidy across cancer types 24 . WA56 and WA50 harbored a frameshifting indel and a probably damaging Significantly mutated genes and pathways in prostate cancer Grasso et al. 8 As described in the main text, we identified three genes as significantly mutated that did not have a previously described role in prostate cancer: MLL2, CDK12, and OR5L1. MLL2 encodes a H3K4-specific histone methyltransferase 25 that is recurrently mutated in multiple cancers including diffuse large B-cell lymphoma 26 , urothelial carcinoma 27 and medulloblastoma 28 , and is a direct coactivator of the estrogen receptor 29 . CDK12, which encodes a transcription elongation-associated C-terminal repeat domain (CTD) kinase 30 , was recently identified as one of nine significantly mutated genes in ovarian serous carcinoma 31 , and silencing of CDK12 has previously been shown to cause resistance to tamoxifen and estrogen deprivation in ERdependent breast cancer models 32 , suggesting a potential role in endocrine resistance in CRPC.
OR5L1 is an olfactory gene that exhibits a higher than average mutation rate as a result of its late replication, arguing against a role in cancer 33 .
Previously, through single gene and focused panel resequencing studies of prostate cancer, mutations in candidates including AR, TP53, CHEK2, KLF6, EPHB2, ZFHX3 (ATBF1), NCOA2, PLXNB1, SPTA1, and SPOP have been reported 2,16,[34][35][36][37][38][39][40] . Through our exome sequencing, we identified recurrent mutations in several genes previously reported to be recurrently mutated in prostate cancer, including AR, TP53, and ZFHX3 (each of which was significantly mutated), as well as SPOP; however no mutations were identified in CHEK2, KLF6 or NCOA2 (previously reported to be mutated in prostate cancer). Importantly, 61 (100%), 51 (84%) and 60 (98%) of the 61 samples had at least 70% of bases with sufficient coverage to call somatic mutations for CHEK2, KLF6 and NCOA2, respectively (Supplementary Table 6), suggesting that the lack of identified mutations is unlikely to be due to inadequate sequencing, and instead suggests that mutations in these genes may be rare, present in a small population of tumour cells or negatively selected for in CRPC. Grasso et al. 9 We also identified 88 significantly mutated canonical pathways out of 880 considered (see Methods & Supplementary Table 17), including 49 with substantial contributions from the nine significantly mutated genes. For example, we identified the 'WNT signaling' KEGG pathway to be significantly mutated (57 somatic mutations, 38 samples, q-value = 1E-6). Half of these mutations occurred in genes other than TP53 and APC, including three missense mutations in CTNNB1 and a splice site mutation in MYC. Additionally, WA57 harbored concurrent nonsense mutation (W509*) and high-level copy loss in SMAD4, a gene which has recently been described as controlling lethal metastasis in CRPC 41 .
We next utilized the matched somatic point mutation and exome copy number data to identify altered subnetworks in a large protein-protein interaction network using HOTNET 42 . Table 18) as significantly mutated in CRPC, including a PTEN interaction network, which was altered in 81% of samples (Supplementary Fig. 11). While 48% of CRPC samples have PTEN mutations, 33% of CRPC samples have mutations in a protein that directly interacts with PTEN, indicating an even broader role for PTEN in prostate cancer pathogenesis and suggesting that mutational status of numerous genes may be required for stratification of therapies targeting the PTEN pathway. For example, we identified a probably damaging R215W mutation in WA57 of  Table 4), suggesting potential roles in prostate cancer progression. Grasso et al. 10 In addition to those significantly mutated genes and pathways, we also identified recurrent mutations in intriguing candidates and pathways. For example, we identified three CRPC samples with mutations in FRY (R100C in WA32, I1480T in WA56, S2510N in WA57), the homologue of the Drosophila gene Furry that encodes a microtubule binding protein required for precise chromosome alignment 43 . Mutations in FRY may promote chromosomal instability in CRPC or result from selection during treatment with docetaxel (a microtubule binding agent), a standard therapy for men with CRPC. In addition, we identified a KRAS G12V mutation in WA42 (ETS -), consistent with our previous reports of rare RAF and RAS family aberrations in ETSprostate cancers 44,45 .

Identification of potential drivers by combined copy number and mutation analysis
To identify potential drivers, we considered genes with recurrent high-level gains or losses present in peaks of global copy number change, and compared results to mutated genes ( Supplementary Fig. 12). Using this approach, AR had the maximum copy number sum (57), with 25 samples showing high-level copy number gain (all CRPC). Likewise PTEN had the minimum copy number sum (-64), with 25 samples showing high-level copy number loss. Both genes also harbored recurrent somatic mutations (Figure 1), supporting the validity of this approach.
Finally, by aCGH of our matched cohort, we confirmed focal deletions in WA7, WA19 and WA10, and identified three additional localized prostate cancers with focal, high copy number loss of CHD1 (Supplementary Fig. 13a).
To further explore the association of CHD1 focal deletion/mutation and ETS status in prostate cancer, we analyzed the association between CHD1 and ETS status in prostate cancer using three prostate cancer aCGH studies (totaling 331 additional cancers) using Oncomine Powertools (http://powertools.oncomine.com). Importantly, each study showed a peak of copy number loss on 5q21, and in each study, all cancers with focal deletions of CHD1 were ETS -(15 of 331 total, 4%) (Supplementary Table 19, Supplementary Fig. 14 48,49 . Additionally, we also identified other tumours with focal deletions involving other genes at 5q21, including PJA2 (high-level copy loss in T65 and T53, and Y505C in WA53) suggesting the existence of other potential drivers at 5q21 (Supplementary Fig. 13a,c&d). Hence, in summary, our integrated analysis identifies deletion or mutation of CHD1 as defining a novel subtype (CHD1 -) of ETSprostate cancer.

ETS2 is both deleted and mutated in prostate cancer
Our copy number profiling data generated here demonstrates that multiple ERG rearrangement positive CRPCs show focal deletions extending telomeric from ERG ( Supplementary Figure 15a), consistent with previous observations in the LuCap35 xenograft and the NCI-H660 prostate cancer cell line (small cell ERG + ) 15,50 . Importantly, WA31 (ERG + through insertion) shows a focal, high copy number loss of ETS2, and our gene expression data demonstrates decreased ETS2 expression in localized cancer and CRPC, with the lowest expression in WA31 (Supplementary Fig. 15b).
Effects on cell invasion were even more pronounced (Fig. 2d, 51,52 , ETS2 may directly compete with other ETS transcription factors for binding to target genes and further investigation will be needed to clarify the role of ETS2 in ETS + and ETSprostate cancers. These results are consistent with distinct ETS genes having oncogenic and potential tumour suppressive roles in prostate cancer 53 .

Identification of chromatin/histone modifying genes mutated in CRPC that interact with the androgen receptor
In addition to MLL2 and CHD1, our integrated analysis identified mutations and copy number aberrations in multiple other genes involved in chromatin/histone modification (Fig. 1)., including MLL2, which was the 7th ranked significantly mutated gene in our data set. The MLL genes (MLL, MLL2 and others) encode histone methyltransferases that function in multi-protein complexes that mediate H3K4 methylation required for epigenetic transcriptional activation 25

3).
Additional deregulated epigenetic modifiers identified in our analysis included the polycomb group gene ASXL2 which was the 17th ranked significantly mutated gene in our data set (p=3.4E-4) and was mutated in 4 samples, with 3 samples harboring nonsense mutations (Y1163* in WA31, Q1104* in WA56 and Q172* in WA23) (Figure 1). We also identified single samples with nonsense mutations in ASXL1 (P749fs in WA52) and ASXL3 (L2240V and R2248* in WA22). Interestingly, ASXL1 is recurrently mutated in myeloid disorders, predominantly through frameshift mutations in the last exon 54 , the same exon affected by the P749fs mutation observed in WA52. Similarly, although UTX (KDM6A), which encodes a histone H3K27 demethylase that complexes with MLL3 25 , is located in a broad region of copy number gain on chr X, it is located at a local copy number minimum, and two samples (WA28 and WA40) show focal high copy loss (Fig. 1). UTX has been shown to be mutated in a number of cancers including renal carcinoma and urothelial carcinoma 25,55,56 . This was of interest, as we previously showed that the histone H3K27 methyltransferase EZH2 is overexpressed in the majority of CRPCs 57 . Grasso et al. 15 As described in the main text, we identified a somatic 2 bp insertion in FOXA1 (S453fs) in the localized prostate cancer sample T12, and 340fs and P358fs indels in DU-145 and LAPC-4 (identified by transcriptome sequencing), respectively. Screening 101 localized and 46 CRPCs Thus, we hypothesized that mutations in FOXA1 may affect proliferation and/or AR signaling, and generated stable LNCaP cells expressing empty vector (LNCaP vector), wild type FOXA1 (FOXA1 wt) and the five FOXA1 mutants observed in our clinical samples as N-terminal FLAG fusions. Western blot and QPCR analyses confirmed equivalent levels of expression of each FOXA1 construct ( Fig. 4b and Supplementary Fig. 17a). The S453fs insertion allele encodes a protein with a predicted molecular weight 49 kDa, similar to wild type FOXA1 (49.2kDa). Grasso et al. 16 In LNCaP cells grown in the presence of 10nM DHT, all FOXA1 mutants, as well as FOXA1 wt, showed significantly increased cell proliferation compared to LNCaP vector (p=0.006 for FOXA1 F400I, p<0.001 for all comparisons to LNCaP vector), while only FOXA1 L388M showed significantly increased growth compared to FOXA1 wt (p=0.005). Expression of FOXA1 wt or mutants had no significant effect on LNCaP proliferation in the absence of androgen (Supplementary Fig. 17b).

Identification and characterization of recurrent mutations in FOXA1
Given the role of FOXA1 as a cofactor for AR signaling, and the reported ability of

Supplementary Discussion
Future studies are needed to understand the temporal development of aberrations identified in this study, to delineate their occurrence during prostate cancer development and progression (including therapy and resistance mechanisms). Importantly, our results provide insight into the ability of sequencing based studies to identify potential therapeutic targets in patients with advanced cancers 62  Similarly, identification of rare lesions such as high-level somatic copy loss of BRCA2 ( Supplementary Fig. 4), high-level copy gain of CDK4 (Supplementary Fig. 7b), or potentially activating mutations of RET (R873W in WA47) may provide rationale for targeted clinical trials or treatment with novel agents. Grasso et al. 19 Methods

Tissue samples and cell lines
Prostate tissues were from the radical prostatectomy series at the University of Michigan and from the Rapid Autopsy Program 64  All CPRC specimens (WA2-WA60) and paired normal tissues were obtained at rapid autopsy from men who died of lethal castrate resistant metastatic disease. Our rapid autopsy protocol has been described in detail previously 64 . Briefly, at the time of autopsy, portions of all cancerous tissue grossly identified, as well as uninvolved organs, were processed by routine formalin fixation and paraffin embedding (FFPE). Corresponding samples were also snap frozen in OCT or chunks. Hematoxylin and eosin (H&E) stained FFPE and frozen sections were reviewed by study pathologists (R.M., L.P.K. and SAT) to identify blocks with highest tumour content (tumour) or lack of cancer or high grade prostatic intraepithelial neoplasia (normal). For each frozen block used, a level was taken for H&E staining, consecutive 3 x 10um sections were cut for DNA isolation, a level was taken for H&E staining, consecutive 3 x 10um sections were cut for RNA isolation, and a final level was taken for H&E staining. All H&E stained levels were reviewed to confirm tumour/normal content before DNA/RNA isolation.
All prostatectomy specimens (N1-N29 and T1-T97) were obtained from treatment naïve men at the time of prostatectomy, where fresh tissue was obtained as part of our standardized procurement protocol. Intervening sections not embedded for routine histological assessment from prostatectomies were quartered and snap frozen in OCT. Evaluation of frozen sections to Grasso et al. 20 identify blocks with the highest tumour content (tumour) or lack of cancer and unremarkable morphology (normal) and RNA/DNA isolation were performed as described above.
The immortalized prostate cancer cell lines 22Rv1, C4-2B, CWR22, DU-145, LAPC-4, LNCaP, MDA-PCa-2B, NCI-H660, PC3, VCaP and WPE1-NB26 (Supplementary Table 11 Table 1) were generated using the Illumina Paired-End Genomic DNA Sample Prep Kit, following the manufacturers' instructions. 3 ug of each genomic DNA was sheared using a Covaris S2 to a peak target size of 250 bp. Fragmented DNA was concentrated using AMPure XP beads (Beckman Coulter), and DNA ends were repaired using T4 DNA polymerase, Klenow polymerase, and T4 polynucleotide kinase. 3' A-tailing with exo-minus Klenow polymerase was followed by ligation of Illumina paired-end adapters to the genomic DNA fragments. The adapter-ligated libraries were electrophoresed on 3% Nusieve 3:1 (Lonza) agarose gels and fragments between 300 to 350 bp were recovered using QIAEX II gel extraction reagents (Qiagen). Recovered DNA was then amplified using Illumina PE1.0 and PE2.0 primers for 9 cycles. The amplified libraries were purified using AMPure XP beads and the DNA concentration was determined using a Nanodrop spectrophotometer. 1 mg of the libraries were hybridized to the Agilent biotinylated SureSelect Capture Library at 65°C for 72 hr or to the Grasso et al. 22 Roche EZ Exome capture library at 47°C for 72 hr following the manuufacturer's protocol. The targeted exon fragments were captured on Dynal M-280 streptavidin beads (Invitrogen), washed, eluted, and enriched by amplification with the Illumina PE1.0 / PE2.0 primers for 8 additional cycles. After purification of the PCR products with AMPure XP beads, the quality and quantity of the resulting exome libraries were analyzed using an Agilent Bioanalyzer.

Somatic point mutation identification by exome capture sequencing
All captured DNA libraries were sequenced with the Illumina GAII Genome Analyzer or the Illumina HiSeq in paired end mode, yielding 80 base pairs from the final library fragments.
The reads that passed the chastity filter of Illumina BaseCall software were used for subsequent analysis. Next, matepairs were pooled and then mapped as single reads to the reference human genome (NCBI build 36.1, hg18), excluding unordered sequence and alternate haplotypes, using Bowtie 65 , keeping unique best hits, and allowing up to two mismatched bases per read. Reads in the tumour that mapped to another location in the genome with three mismatches were excluded from further consideration. Likely PCR duplicates were removed by removing reads that have the same match interval on the genomic sequence. Individual basecalls with Phred quality less than Q20 were excluded from further consideration.
A mismatched base (SNV) was identified as a somatic mutation only when 1) it had six reads of support (this cut-off was selected based on Sanger validation rates in T12 [ Supplementary Fig. 1]), 2) it was in at least 10% of the coverage at that position in the tumour, 3) it was observed on both strands, 4) there was 8X coverage in the matched normal, and 5) it did not occur in the matched normal sample in more than two reads and 2% of the coverage (to ensure that somatic variants are not filtered out due to tumour contamination in the normal, we Grasso et al. 23 retained variants present in 2-4% of the coverage in the matched normal, if they were in at least 20% of the coverage in the tumour). SNVs were excluded from further consideration as somatic mutations if 1) they did not fall within 50 bases of a target region, 2) they occurred in any two matched normal samples in at least two reads and 2% of the coverage, or 3) they occurred in another tumour and its matched normal sample in two reads and 4% of the coverage.

Identification of coding indels in exome capture data
The methodology for identifying indels in exome capture data was adapted from 66 with minor modifications. Reads for which Bowtie was unsuccessful in identifying an ungapped alignment were converted to fasta format and mapped to the target regions, padded by 200 bases on either side, with cross_match (v0.990329, http://www.phrap.org), using parameters -gap_ext -1 -bandwidth 10 -minmatch 20 -maxmatch 24. Output options were -tags -discrep_listsalignments. Alignments with an indel were then filtered for those that: 1) had a score at least 40 more than the next best alignment, 2) mapped at least 75 bases of the read, and 3) had two or fewer substitutions in addition to the indel. Reads from filtered alignments that mapped to the negative strand were then reverse-complemented and, together with the rest of the filtered reads, remapped with cross_match using the same parameters (to reduce ambiguity in called indel positions due to different read orientations). After the second mapping, alignments were refiltered using criteria 1-3. Reads that had redundant start sites were removed as likely PCR duplicates, after which the number of reads mapping to either the reference or the non-reference allele was counted for each. An indel was called if there were at least six non-reference allele reads making up at least 10% of all reads at that genomic position. Indels were reported with respect to genomic coordinates. For insertions, the position reported is the last base before the Grasso et al. 24 insertion. For deletions, the position reported is the first deleted base. Indel somatic mutation candidates were excluded from further consideration if 1) they did not occur on both strands, 2) they did not fall within 50 bases of a target region, 3) there wasn't 8X coverage in the matched normal at that position, 4) they occurred in the matched normal sample in more than 2 reads and 4% of the coverage, 5) they occurred in any two matched normal samples, or 6) they occurred in any single matched normal sample in more than 2 reads.

Annotation
We annotated the resulting somatic mutations using CCDS transcripts wherever possible.
If no CCDS transcript was available, we used the coding regions of RefSeq transcripts. HUGO gene names were used. The impact of coding nonsynonymous amino acid substitutions on the structure and function of a protein was assessed using PolyPhen-2 67 . We also assessed whether the somatic variant was previously reported in dbSNP or COSMIC v56 68 .

Calculation of somatic mutation rates
The somatic mutation rate was calculated as described 2 . We identified a base as "covered", if there was at least 14X total coverage after PCR duplicate removal in the tumour and 8X total coverage after PCR duplicate removal in the matched normal sample. We considered only  In order to filter out regions of possible amplification or deletion, we used exon coverage ratios to infer copy number changes, following the approach of 14 . Resulting SNV candidates were not used for estimation of tumour content if the segmented log-ratio exceeded 0.25 in absolute value. Candidates on the X and Y chromosomes were also eliminated because they were unlikely to exist in 2-copy genomic regions.
Using this set of candidates, we fit a binomial mixture model with two components using the R package flexmix, version 2.2-8 69  sample. Thus, the estimated binomial proportion as obtained from the mixture model was doubled to obtain an estimate of tumour content in each sample.

Determination of significantly mutated genes and pathways
The determination of significantly mutated genes and pathways was done as described 2,70 using methodology based on that of Getz et al. 71 and Ding et al. 72 . Before doing the calculations, we selected one of the three samples derived from distinct metastatic sites from the same individual (WA43) for inclusion in the sample set in order to ensure that the requirement of independence was met for the set of considered mutations. We selected WA43-44, because it contained all of the recurrent somatic mutations that occurred in WA43-27 or WA43-71, along with additional recurrent mutations not contained in the other two. We also excluded hypermutated sample WA16. In this approach, significantly mutated genes are identified based on the observed number of mutations for each sequence context-based mutation class (CpG, other C:G, A:T, and indels), the sample-specific and class-specific background mutation rates, and the number of covered bases per gene (Supplementary Table 6 Benjamini-Hochberg procedure for controlling False Discovery Rate (FDR).
We repeated this analysis to consider significantly mutated pathways, considering a list of 880 gene sets corresponding to the set of canonical pathways used in Gene Set Enrichment Analysis (GSEA). For this analysis, we tabulated the number of mutations and the number of covered bases in all component genes of each gene set and the total number of covered bases in the set.
As in the single-gene analysis, the mutation counts were broken down into the context-based mutation classes (CpG, other C:G, and A:T), and then the P-value and subsequent q-value were calculated.

Sanger sequencing to validate somatic point mutations and indels
Various genomic locations nominated for somatic point mutations and indels were amplified from whole genome amplified DNA 73 from corresponding matched normal-tumour tissue pairs or cell lines. Briefly, fifty ng of input genomic DNA was subjected to fragmentation, library preparation and amplification steps using Genomeplex-Complete Whole Genome Amplification Kit (Sigma-Aldrich) according to manufacturer's instructions. The final whole genome amplified DNA was purified by AMPure XP beads (Beckman-Coulter) and quantified by a Nanodrop spectrophotometer (Thermo Scientific). Fifty ng of DNA were used as template in PCR amplifications with HotStar Taq DNA polymerase (Qiagen) with the suggested initial denaturation and cycling conditions. Primer sequences were as described 9 and are based on human hg18, March 2006 assembly. Primers for FOXA1 can be found in Supplementary Table   20. The PCR products were subjected to Sanger sequencing by the University of Michigan DNA Sequencing Core after treatment with ExoSAP-IT (GE Healthcare) and sequences were analyzed using MacVecotr software (MacVector).

Exome Copy Number Analysis
Copy number aberrations were quantified and reported for each gene as the segmented normalized log2-transformed exon coverage ratios between each tumour sample and its matched normal (Supplementary Table 7) as described 14 .
We used sample-specific cutoffs, based on estimated tumour content, to define regions of gain and loss, as follows. For a sample with tumour percentage P, genomic regions with N copies in cancer cells and 2 copies in normal cells would be predicted to give log-ratios centered at log2(N*P + 2*(1-P))-1. For each sample, using its estimated tumour content we computed the predicted locations of these N-copy peaks in the distribution of log-ratios, and chose cutoffs to fall between these predicted peaks.
To define high-level gains (i.e., greater than 3 copies), we computed the weighted average of the 3-copy and 4-copy predicted peaks with weights 0.25 and 0.75, respectively.
Similarly, to define low-level gains (i.e., greater than 2 copies), we computed the weighted average of the 2-copy and 3-copy predicted peaks, using the same weights. These weighted averages were used as cut-offs to define high-level gain and low-level gain, respectively. Next, the negatives of the cutoffs for high-level gain and low-level gain were used as the cut-offs for high-level loss (two-copy loss) and low-level loss (single-copy loss), respectively. Histograms of the distributions of segmented log2 copy number ratios were then examined ( Supplementary   Fig. 7) and cutoffs adjusted manually in cases in which this algorithmic approach appeared to misclassify large numbers of genomic regions (due to lower tumour content, multiple clones, severe aneuploidy, etc. All cutoffs with estimated tumour percentages are given in Supplementary Table 8. The resulting copy number alterations were reported for all sixty one Grasso et al. 29 prostate cancer tumours with +2 representing high-level (> 1copy) gain, +1 representing 1 copy gain, 0 representing no change, -1 representing 1 copy loss, and -2 representing high-level (>1 copy) loss (Supplementary Table 9). To identify potential drivers, we summed all called copy number alterations across all samples and identified genes with the maximum number of high level gains or losses occurring in peaks summed copy number gain or loss, respectively. For all analyses and visualizations, WA43-27 and WA43-71 were excluded.

Identification of Significantly Mutated Protein-Protein Interaction Subnetworks
We used HotNet 42 to find subnetworks of a large protein-protein interaction network containing a significant number of mutations and copy number alterations (CNAs). The input to HotNet is a dataset of matched somatic mutations and copy number alterations for a set of tumour samples. The output of HotNet is a list of subnetworks, each containing at least n genes.
HotNet employs a two-stage statistical test to assess the significance of the output. In the first stage the p-value for the number of subnetworks in the list is computed. In the second stage the false discovery rate (FDR) of the list of subnetworks is estimated. At the end, the significance of each individual subnetwork in the list is assessed by comparison to known pathways and protein complexes.
We analyzed the combined somatic mutations (Supplementary Table 4) and CNAs (generated from exome data [Supplementary Table 9]), considering only high-level (> 1 copy) gains and two-copy losses, for the 47 metastatic samples (we did not include hyper-mutated sample WA16 and only included one of the three metastatic sites from the same individual, WA43-44). We removed CNAs for which the sign of aberration was not consistent in at least 90% of the altered samples. We used the interaction network derived from the Human Protein Grasso et al. 30 Reference Database (HPRD) 74 . For the statistical test, we generated random aberrations as follows. We simulated mutations using the estimated background mutation rate (1.97x10 -6 ). We simulated CNAs from the observed distribution of CNA lengths, permuting their positions. In this way artifacts resulting from functionally related genes that are both neighbors on the network and close enough on the genome (and thus affected by the same CNA) are minimized.
We also discarded subnetworks reported by HotNet that contain 3 or more genes in the same CNA in one or more samples. Moreover, for subnetworks with two genes g1, g2 in the same CNA in 1 or more samples, we removed the genes that were not reported when alterations in either g1 or g2 are removed.
Using the approach above, HotNet identifies 28 candidate subnetworks containing at least 10 genes (p < 0.01) with FDR = 0.32. A total of 24 subnetworks remained after CNA filtering.
We then compared those 24 subnetworks with known pathways in the KEGG database and with protein complexes from PINdb 75 . Of the 24 subnetworks, 14 had statistically significant (p<0.05 after Bonferroni correction) overlap with at least one KEGG pathway or protein complex (Supplementary Table 18).

RNA isolation and cDNA synthesis
Total RNA was isolated from frozen prostate tissue samples (for gene expression analysis and qPCR) and cell lines (for transcriptome sequencing, qPCR/expression profiling from cell lines) using Trizol (or Qizol [Qiagen]) and an RNeasy Kit (Invitrogen) with DNase I digestion according to the manufacturer's instructions. RNA integrity was verified on an Agilent Bioanalyzer 2100 (Agilent Technologies). cDNA was synthesized from total RNA using Grasso et al. 31 Superscript III (Invitrogen) and random primers (Invitrogen).

Transcriptome library preparation and sequencing
Next generation RNA sequencing was performed on 11 prostate cell lines according to Illumina's protocol using 2ug of total RNA. RNA integrity was measured using an Agilent 2100 Bioanalyzer, and only samples with a RIN score >7.0 were advanced for library generation.
PolyA+ RNA was selected for using Sera-Mag oligo(dT) beads (Thermo Scientific) and

Somatic point mutation identification in transcriptome sequence data
Transcriptome short reads were trimmed to remove the first two bases and as many bases as necessary to ensure the read length was less than 40bp. Trimmed short read sequences were mapped to the reference human genome (NCBI build 36.1, hg18), excluding unordered sequence and alternate haplotypes, and the 2008 Illumina splice junction set using Bowtie 65 in single read Grasso et al. 32 mode keeping unique best hits and allowing up to two mismatched bases. Matepairs from paired end runs were pooled and treated as single reads. Likely PCR duplicates were removed by removing reads that have the same match interval on the genomic sequence or an exon junction.
Individual basecalls with Phred quality less than Q20 were excluded from further consideration.
A mismatched base (SNV) was identified as a candidate somatic mutation when it had three reads of support and was in at least 10% of the coverage at that position in the tumour. Less stringent criteria were applied for nominating candidate somatic mutations in the transcriptome as compared to the exome capture data, since only variants in the transcriptome recurrent to known somatic mutations were further considered (see below). SNVs were excluded from further consideration as recurrent somatic mutations if 1) they occurred in any two matched normal exomes in at least two reads and 2% of the coverage, or 2) they occurred in another tumour exome and its matched normal exome in two reads and 4% of the coverage.

Identification of coding indels in transcriptome data
The methodology for identifying indels in transcriptome data was adapted from 66 . Reads for which Bowtie was unsuccessful in identifying an ungapped alignment were converted to fasta format and mapped to the set of full-length CCDS transcripts, padded by 32 genomic bases on either side, with cross_match (v0.990329, http://www.phrap.org), using parameters -gap_ext -1 -bandwidth 10 -minscore 24 -minmatch 16 -maxmatch 24. Output options were -tags -discrep_lists -alignments. Alignments with an indel were then filtered for those that: 1) had a score at least 20 more than the next best alignment; and 2) had two or fewer substitutions in addition to the indel. Reads from filtered alignments that mapped to the negative strand were then reverse-complemented and, together with the rest of the filtered reads, remapped with Grasso et al. 33 cross_match using the same parameters (to reduce ambiguity in called indel positions due to different read orientations). After the second mapping, alignments were re-filtered using criteria 1) and 2). Reads that had redundant start sites were removed as likely PCR duplicates, after which the number of reads mapping to either the reference or the non-reference allele were counted for each. An indel was called if there were at least four non-reference allele reads making up at least 10% of all reads at that transcript position. Indels were reported with respect to genomic coordinates. For insertions, the position reported is the last base before the insertion.
For deletions, the position reported is the first deleted base. Indel somatic mutation candidates were excluded from further consideration if they were present in dbSNP132, or if they occurred in a single read in any two matched normal exome samples or in a single matched normal exome sample with two or more reads. Identified indel variants are given in Supplementary Table 13.

Identification of transcriptome somatic SNVs recurrent to known somatic variants
The somatic mutations identified in the exome data in this study (excluding the eight that did not validate by Sanger sequencing) were combined with the confirmed somatic variants in COSMIC v56 to yield a comprehensive somatic mutation dataset. A transcriptome SNV was considered recurrent to a known somatic variant, if it resulted in the same nucleotide change, amino acid change, or if it disrupted the same amino acid. Identified variants recurrent to our exome data are given in Supplementary

Gene Expression Microarray Analysis
Gene expression microarray analysis of the same prostate tissue samples subjected to aCGH (Supplementary Table 10) was performed using Agilent Whole Human 44k element arrays (1x44k or 4x44k format) as described 76 . RNA from indicated prostate samples were labeled with Cy-5 (test channel) and hybridized against Cy-3 (reference channel) labeled pooled benign prostate RNA (Clontech). Arrays were scanned using an Agilent Model G2505B scanner, and data was extracted using Agilent Feature Extraction software. Control probes were removed from all arrays and the LogRatio for all probes, which were used for subsequent analysis, were converted to log (2). Although the 1x44k and 4x44k arrays have the same probes, the 4x44k arrays have 10 replicates of some probes. Thus, to generate a final data set, the median value of replicated probes was used for 4x44k arrays. The final data set (including benign prostate, localized prostate cancer and CRPC) was uploaded into a custom instance of Oncomine for automated analysis. In Oncomine, the dataset was median centered (per array) prior to indicated analyses. Copy number and gene expression data is also available from GEO (GSE35988). were considered ETS + , samples with CHD1 outlier under-expression were considered CHD1and samples with outlier over-expression of SPINK1 were considered SPINK1 + (see Supplementary Table 19).

ETS2
Full length wild type ETS2 with N-terminal HA-tag was PCR amplified and cloned into pCR8/GW/TOPO vector (Invitrogen). ETS2 R437C was generated using the Quick changemutagenesis kit (Stratagene). ETS2 wildtype and R437C were transferred into pLenti4-V5 DEST Grasso et al. 37 vector (Invitrogen). After confirmation of the insert sequence, lentiviruses were generated by the University of Michigan Vector Core. VCaP cells were infected and stably expressing ETS2 wild type, ETS2 R437C mutant and lacZ control were generated by selection with Zeocin (Invitrogen). ETS2 expression was confirmed by qPCR for ETS2 expression and western blotting with anti-HA antibody as above. For proliferation assays, 50,000 cells (n=4) were plated per well in 24-well poly-lysine coated plates, and cells were harvested and counted at the indicated time points by Coulter counter (Beckman Coulter, Fullerton, CA). For in vitro migration and invasion, 2.0x10 5 cells (migration n=8; invasion n=12) were placed in the top chamber with a noncoated membrane or Matrigel coated membrane, respectively (24-well insert; pore size 8µm; BD Biosciences). In both the assays, cells were plated in medium without serum, and medium supplemented with 10% serum was used as a chemoattractant in the lower chamber. Cells were incubated for 48hr and cells that did not migrate or invade through the pores were gently removed with a cotton swab. Cells on the lower surface of the membrane were stained with crystal violet and counted.

FOXA1
FOXA1 wildtype and FOXA1 mutants (S453fs, G87R, L388M, L455M and F400I) were cloned and inserted into pCDH (System Biosciences), which has been modified to express an Nterminal FLAG tag and puromycin resistance. Lentiviruses were generated in 293FT cells using the ViraPower Lentiviral Expression System (Invitrogen). LNCaP cells were infected with the generated viruses (or empty control virus) and stable pooled populations were selected with puromycin. Expression was confirmed by western blotting with anti-FLAG antibody (Sigma) or qPCR for FOXA1 expression as above, and FOXA1 primers are in Supplementary Table 20. qPCR experiments were performed in triplicate, and FOXA1 expression was normalized to GAPDH. For proliferation, cells were starved for 24 hours in phenol red-free media media with Grasso et al. 40 1% charcoal-dextran stripped serum, and grown in media with 1% charcoal-dextran stripped serum +/-10nM DHT. Relative cell numbers were measured in quadruplicate by WST-1 assays at indicated time points following the manufacturer's instructions (Roche).
Gene expression microarray analysis was performed as above, using LnCAP cells expressing empty vector, FOXA1 wild type, or FOXA1 mutants as just described. Cells were starved in 1% charcoal stripped media for 24 hours and treated with 10 nM DHT or vehicle for 48 hours. RNA was isolated using Qiazol. All samples were hybridized against vehicle stimulated vector control in duplicate. Probes not passing filtering in both duplicate hybridizations were excluded from further analysis, and remaining probes were averaged. For each set of FOXA1 wildtype or mutant hybridization (duplicates of DHT and vehicle treated), DHT vs. vehicle stimulated ratios for each probe passing filtering on all four arrays were computed. Probes were filtered to include only those with average LogRatio (converted to log base 2) of >1 or <-1 in the DHT vs. vehicle stimulated pair. Clustering of probes using centroidlinkage clustering was performed using Cluster 3.0 and heatmaps were generated using JavaTreeview.
In parallel, FOXA1 wildtype and FOXA1 mutant (S453fs; resulting from chr14:37130381insCC observed in T12) ORFs were generated by gene synthesis (Blue Heron) and cloned into the pLL_IRES_GFP lentival vector. Lentiviruses (and pLL_IRES_GFP expressing LACZ as control) were generated by the University of Michigan Vector Core. LNCaP cells were transduced in the presence of 4µg/mL polybrene (Sigma). After 72 hours, GFP+ cells were sorted at the University of Michigan flow cytometry core. Cells were genotyped to confirm identify. GFP fluorescence was monitored every other day. Soft agar colony forming assays were performed as described 51 , except colonies were counted and photographed without staining. Grasso et al. 41 For xenograft experiments, four week-old male SCID C.B17 mice were procured from a mice breeding colony at University of Michigan (maintained by K.J.P.). Mice were anesthetized using a cocktail of xylazine (80 mg/kg IP) and ketamine (10 mg/kg IP) for chemical restraint.

Indicated LNCaP cells (2 million cells per implantation site) as above (or parental LNCaP cells)
were suspended in 100ul of 1X PBA with 20% high concentration Matrigel (BD Biosciences).
Cells were implanted subcutaneously on both sides into the flank region. Tumour growth was monitored bi-weekly by using digital calipers in LNCaP FOXA1 wildtype (n=9), LNCaP FOXA1 S453fs (n=10) and parental LNCaP (n=6) groups. Tumour volumes were calculated using the formula (π/6) (L×W2), where L= length of tumour and W= width. All procedures involving mice were approved by the University Committee on Use and Care of Animals (UCUCA) at the University of Michigan and conform to their relevant regulatory standards.

DLX1
For DLX1 immunoblotting, prostate tissues were homogenized in NP40 lysis buffer containing 50 mm Tris-HCl (pH 7.4), 1% NP40 (Sigma), and complete proteinase inhibitor mixture (Roche). Western blotting with ten micrograms of each protein extract was performed as above. Transferred membrane was incubated for 1 h in blocking buffer and over-night with anti-DLX1 rabbit polyclonal antibody (PTG laboratory, #13046-1-AP, 1:1000 dilution). After washing three times with TBS-T buffer, the membrane was incubated with horseradish peroxidase-linked donkey anti-rabbit IgG antibody (GE Healthcare, 1:5,000) for 1 h at room temperature prior to visualization by enhanced chemiluminescence (GE Healthcare). To monitor equal loading, the membrane was re-probed with anti-β-Actin mouse monoclonal antibody (1:30,000 dilution; Sigma, # A5316). Grasso et al. 42 qPCR was performed on 10 benign prostate tissues (included in gene-expression profiling), 55 localized prostate cancers (including 32 samples subjected to gene-expression profiling) and 7 CRPCs (including 6 samples subjected to gene-expression profiling) as above.
The amount of DLX1 in each sample was normalized to the average of GAPDH and HMBS for each sample. Primers for DLX1 are given in Supplementary Table 20; GAPDH and HMBS primers were as described 87 . All oligonucleotide primers were synthesized by Integrated DNA Technologies.   T8  T12  T32  T90  T91  T92  T93  T94  T95  T96  T97  Genome wide copy number plots by exome sequencing for two hypermutated CRPC samples (WA16 and WA48). For each sample, the Log2 copy number ratio between tumor and matched normal is plotted for each targeted exon and then ordered by genomic location. Genes with 1 copy gain/loss are indicated by red and blue points, respectively, and those with > 1 copy gain/loss are indicated by orange and cyan points, respectively. The location of focal high level deletions of MSH2 in WA16 and BRCA2 in WA48 are indicated.      Figure 7. Genome wide copy number analysis by exome sequencing and identification of 1 copy and >1 copy gains/losses. For every sample, segmented normalized log2-transformed exon coverage ratios between each tumor sample and its matched normal were computed. Sample-specific cutoffs were generated, based on estimated tumor content, to define regions of 1 copy gain and loss, and >1 copy gain/loss. Manual adjustment of cutoffs was performed for samples where the algorithmic approach appeared to misclassify large numbers of genomic regions. a. Distribution histogram of all Log2 copy number ratios (tumor to normal) for each targeted exon in WA15. The automated cutoffs for no change (gray), 1 copy loss (blue) and gain (red), and > 1 copy loss (cyan) and gain (orange) are indicated. Manual adjustment was not performed for this sample. b. Genome wide copy number aberrations for WA15. The Log2 copy number ratio (tumor to normal) for each targeted exon in WA15, ordered by genomic location is shown. Using the cutoffs in a, genes with called 1 copy gain/loss are indicated by red and blue points, respectively, and those with >1 copy gain/loss are indicated by orange and cyan points, respectively. The location of a focal amplification including CDK4 on chr12 is indicated.          Figure S12 b c Supplementary Figure 12. Identification of high level, focal copy number aberrations in prostate cancer. a. Genome wide copy number analysis of each sample was performed using exome sequencing. For all genes, the sum of copy number calls (+/-1: one copy gain or loss, respectively; +/-2: high level copy gain/loss, respectively) across all samples is plotted and ordered by genome location. b. As in a, but only the sum of high level copy gains/losses (+/-2) is plotted. c.    T65  T56  T53  T45  T43  T42  WA7  WA47  WA42 N9  N21  N16  N22  N23  N11  N7  N26  N19  N5  N2  N29  N10  N20  N27  N6  N8  N15  N17  N1  N18  N24  N12  N4  N28  N13  N25  N14  WA4  WA39  WA13  T19  T62  T8  T60  WA53  T54  WA2  T11  T17  WA28  WA55  T50  T5  T59  T44  T75  T68  T51  T83  T69  WA6  T70  T12  T46  WA54  T82  WA31  T64  T41  WA18  WA37  T73  T48  WA40  WA16  T63  T85  T39  WA22  WA24  WA23  T20  WA26  T57  T52  T3  T9a  T10  WA7  T45  WA10  T42  T43  T65  WA19  T53  T37  T58  T55  T26  T61  T67  WA32  T40  WA5  WA20  T21  T6  WA33  WA14  WA47  T56  T27  T49  T29  T32  T25  WA35  T66  T31  T24  T47  WA3  WA25  T1  T7  WA46  WA29  WA11  WA42  WA30 ETS +  Genome wide copy number analysis of high-grade localized prostate cancer and castrate resistant prostate cancer by exome sequencing identified a peak of copy number loss on chr 5q21 centered on CHD1. A subset of samples used for exome sequencing, and additional benign prostate tissue (N1-29) and localized prostate cancers were used for matched aCGH and gene expression profiling. a. Genome wide analysis by aCGH identified a similar peak of copy number loss on 5q21 (upper panel, sum log2 copy number across all samples plotted) centered on CHD1. The expanded view is as in Figure 2a, except the area (absolute Log2 ratio) and color intensity (Log 2 ratio; copy number loss in blue) of each box are proportional to binned copy number for that gene according to the legend. ETSand ETS + samples are indicated in black or gray type, respectively. Samples with focal deletions of CHD1 (CHD1 -) or other genes within 5q21 (5q21 -) by aCGH are indicated with green or red background, respectively, according to the legend. The adjoining plot shows the genome wide copy number plot for T56, which harbors a focal, high level deletion on 5q21 including CHD1. b. Co-expression of CHD1 and ETS family members. Heatmap of CHD1, ETS genes (ERG, ETV1, ETV5) and SPINK1 gene expression. Samples are stratified by benign prostate tissue and prostate cancer (including localized and CRPC). ETS and CHD1 status was determined, with black and green indicating ETS + and CHD1 -, respectively. c. Genome wide copy number plot for T65, which shows focal, high level deletion of 5q21, including PJA2, but not CHD1. d. Expression of PJA2 stratified by benign prostate tissues (orange), localized prostate cancers (cyan) and CRPCs (black). T65 is indicated in red.    Only samples with at least one gene in the region with a called copy number gain/loss are shown, and missing boxes indicate that gene has no called copy number gain/loss. Mutations in ETS2 are indicated according to the legend. b. Gene expression profiles from benign prostate tissues (n=28, orange), localized prostate cancer (PC, n=59, cyan), and metastatic castrate resistant prostate cancer (CRPC, n=35, black), including samples subjected to exome sequencing, were loaded into Oncomine for automated analysis. Expression of ETS2 is shown, including for sample WA31, which harbors a focal, high copy loss of ETS2. c. VCaP prostate cancer cells (ERG + ) stably expressing wild type (wt) ETS2 (black) or ETS2 R437C (yellow) with N-terminal HA tag, or LACZ as control (purple), were generated using lentiviruses (see Fig. 2). qPCR for ETS2 expression was performed for each stable line, and the amount of ETS2 was normalized to GAPDH. Normalized ETS2 expression is plotted relative to LACZ control. Mean of normalized ETS2 expression + S.E. (n=3) are plotted. The inset shows western blotting of the same samples with anti-HA to confirm protein expression.   ASH2L and MLL. a. ). a. Reverse immunoprecipitation using anti-FOXA1 (positive control), an antibody against MLL2, two anti-ASH2L antibodies, an antibody against MLL, or IgG control, with Western blotting for androgen receptor (AR). 1% whole lysate was used as control. b. VCaP cells were treated with siRNAs against ASH2L or MLL (or non-targeting as control). qPCR for ASH2L or MLL expression (normalized to the average of ACTB and GAPDH) confirmed knockdown (n=3, + S.E.), prior to androgen stimulation experiments (see Fig. 3b).   . 4). qPCR for FOXA1 expression was performed for each stable line, and the amount of FOXA1 was normalized to GAPDH. Normalized FOXA1 expression is plotted relative to vector control. Mean of normalized FOXA1 expression + S.E. (n=3) are plotted. b. Cell proliferation in 1% charcoaldextran stripped serum was measured by WST-1 colorimetric assay (absorbance at 450nM) at the indicated time points. c. Soft agar colony forming assays using LNCaP cells stably expressing LACZ (control, purple), or Nterminally HA-tagged FOXA1 wild type (wt, black) or FOXA1 S453fs (green) generated through lentiviral infection. Representative photographs and quantification of colonies formed in the absence (white) or presence of 1nm R1881 are shown. For b and c, mean + S.E. (n=3) are plotted; * indicates p=0.05 from two tailed t-test.