Purpose:

Cross-resistance renders multiple lines of androgen receptor (AR) signaling inhibitors increasingly futile in metastatic castration-resistant prostate cancer (mCRPC). We sought to determine acquired genomic contributors to cross-resistance.

Experimental Design:

We collected 458 serial plasma cell-free DNA samples at baseline and progression timepoints from 202 patients with mCRPC receiving sequential AR signaling inhibitors (abiraterone and enzalutamide) in a randomized phase II clinical trial (NCT02125357). We utilized deep targeted and whole-exome sequencing to compare baseline and posttreatment somatic genomic profiles in circulating tumor DNA (ctDNA).

Results:

Patient ctDNA abundance was correlated across plasma collections and independently prognostic for sequential therapy response and overall survival. Most driver alterations in established prostate cancer genes were consistently detected in ctDNA over time. However, shifts in somatic populations after treatment were identified in 53% of patients, particularly after strong treatment responses. Treatment-associated changes converged upon the AR gene, with an average 50% increase in AR copy number, changes in AR mutation frequencies, and a 2.5-fold increase in the proportion of patients carrying AR ligand binding domain truncating rearrangements.

Conclusions:

Our data show that the dominant AR genotype continues to evolve during sequential lines of AR inhibition and drives acquired resistance in patients with mCRPC.

This article is featured in Highlights of This Issue, p. 4457

Translational Relevance

Using serial plasma cell-free DNA (cfDNA) collections, we show that resistance against the widely used androgen receptor (AR) signaling inhibitors abiraterone and enzalutamide in metastatic castration-resistant prostate cancer (mCRPC) is associated with Darwinian selection for highly aggressive AR genotypes, including increased copy number and increased frequency of ligand binding domain (LBD) truncating rearrangements. No recurrent treatment-driven changes were observed in 71 other prostate cancer–associated genes. Our results support the continued development of novel strategies targeting the AR, as it remains a relevant clinical target even after multiple lines of AR signaling inhibition. Plasma circulating tumor DNA fraction was the strongest predictor of treatment outcomes, surpassing other clinical characteristics and underscoring its potential as a prognostic biomarker. Our study demonstrates the value and limitations of serial plasma cfDNA biopsies in the study of mCRPC treatment resistance, and establishes a framework for future investigations.

Metastatic prostate cancer is a lethal disease, but can be transiently controlled by inhibition of androgen receptor (AR) signaling (1). This inhibition is achieved clinically via suppression of androgen synthesis or direct antagonism of the AR ligand binding domain (LBD). To overcome these treatments, prostate cancers develop augmented AR signaling through mutations altering AR ligand binding affinity, AR amplification, or LBD truncation (2–4).

Several potent AR signaling inhibitors have clinical activity in AR-augmented prostate cancer (5), but cross-resistance is common (6, 7). Adaptive mechanisms can drive renewed resistance (8, 9), but most posttreatment tumors retain AR signaling, suggesting that the pathway remains an oncogenic driver (10). The extent to which the AR gene and the wider somatic genome continue to evolve and drive acquired resistance over sequential lines of AR-targeted therapy remains unknown, largely because of the difficulty of serially biopsying a treatment-standardized cohort of patients.

Recently, we reported primary outcomes from a phase II trial (NCT02125357) of sequential AR inhibitor treatment in metastatic prostate cancer that has already progressed on conventional androgen deprivation therapy (a clinical state described as metastatic castration-resistant prostate cancer, mCRPC; ref. 7). We analyzed pretreatment plasma cell-free DNA (cfDNA) samples collected from the 202 trial patients and demonstrated that circulating tumor DNA (ctDNA) fraction and key genomic alterations were independently prognostic for time to progression (TTP; ref. 11). Our prior work was limited to one time point, had short follow-up, and did not address acquired resistance. Here, we leverage this clinical trial population to present the first comprehensive analysis of longitudinal cfDNA collections (n = 458) from a large standardized mCRPC cohort. We show that most somatic driver alterations are consistently detected in ctDNA over sequential clinical progressions. However, temporal shifts in cancer cell populations represented in ctDNA were common, and the AR locus was highly plastic, evolving toward aggressive genotypes associated with increased AR activation and ligand independence. Collectively, this work provides insight into real-time evolution of the metastatic niche during treatment and proposes a genomic explanation for resistance to the most commonly used mCRPC therapies.

Patient sample cohort and clinical endpoints

Blood samples were collected as part of an open-label randomized phase II crossover trial (NCT02125357) conducted in six centers in British Columbia, Canada (7). In the published trial, patients with newly diagnosed mCRPC were randomly assigned 1:1 to receive abiraterone 1,000 mg orally daily + prednisone 5 mg twice daily followed by enzalutamide 160 mg orally daily (arm A), or the opposite sequence of enzalutamide followed by abiraterone + prednisone (arm B). Patients were to crossover from the first treatment assigned (first-line treatment) to the second treatment (second-line treatment) at the first occurrence of either confirmed PSA progression, radiographic progression as assessed by CT scan (RECIST 1.1 criteria) and/or bone scan (PCWG2 criteria; ref. 12), unequivocal clinical disease progression, or unacceptable toxicity. The study was designed and conducted in accordance with Good Clinical Practice and the Declaration of Helsinki. All patients provided written informed consent, including for collection and analysis of peripheral whole blood samples. Approval for this study was granted by the University of British Columbia Research Ethics Board (Number H14-00738 and H18-00944). Patients were required to have histologically proven prostatic adenocarcinoma with metastatic disease on CT scan, MRI, and/or bone scan, and a rising PSA (PSA progression per PCWG2 criteria; ref. 12) with castrate levels of testosterone (≤1.7 nmol/L) with ongoing medical castration or previous bilateral orchiectomy. Prior use of CYP17A1 inhibitors or novel AR inhibitors was prohibited, whereas prior use of docetaxel for castration-sensitive disease was allowed. Exclusion criteria included contraindications to abiraterone and enzalutamide per the manufacturer's label and Eastern Cooperative Oncology Group performance status of > 2. The presence of visceral metastasis and pain requiring opioid analgesia were allowed. Patients were required to meet the same eligibility criteria to crossover to second-line treatment.

PSA response was defined as ≥30% PSA decline from baseline confirmed on repeat measurement ≥28 days later. Time to PSA progression (on first- or second-line therapy) was defined as the time from start of therapy to confirmed PSA progression. Confirmed PSA progression was defined as an increase of 2 μg/L and 25% from nadir confirmed by subsequent rising PSA ≥ 28 days later. For patients with no PSA decline, PSA progression was defined as an increase of 2 μg/L and 25% from baseline after ≥12 weeks of treatment. Patients without PSA progression were right-censored at the last treatment date. Overall survival was defined as the time from start of first-line therapy to time of death from any cause, or last follow-up (censored).

TTP (on first- or second-line therapy) was defined as the time from start of therapy to confirmed PSA progression, radiographic progression (PCWG2 criteria), clinical progression, or death from prostate cancer, whichever occurred first. If no events occurred, patients were right-censored at the last treatment date.

All Cox regression analyses, associated confidence intervals (CI), and Kaplan–Meier curves were calculated using R (version 3.6.0) with the “survival” package (version 2.44.1.1). CIs for PSA response rates, Pearson's χ2 tests, rank-sum tests, and Fisher exact tests were calculated using Julia (version 1.1.0) with the “HypothesisTests” package (version 0.8.0).

Blood collection, DNA isolation, and sequencing

Whole blood was collected in either 4 × 6 mL EDTA tubes (2014–2016) or 2 × 9 mL Streck Cell-Free DNA BCT tubes (2016–2019). Samples collected in EDTA tubes were kept on ice or at 4°C until processing, which occurred within 1 hour. For EDTA tubes, whole blood was aliquoted first and remaining samples were then centrifuged at 1,600 rpm for 10 minutes at 4°C, after which plasma was transferred to a new tube and spun for 10 additional minutes at 1,600 rpm. Samples collected in Streck tubes were maintained at room temperature prior to and during processing. Streck tubes were centrifuged at 1,600 rcf for 15 minutes, buffy coat was aliquoted, and plasma was transferred to a new tube and spun for 10 additional minutes at 5,000 rcf (or maximum attainable speed).

Plasma aliquots were stored at −80°C prior to DNA extraction. cfDNA was extracted from up to 6 mL of plasma with the QIAGEN Circulating Nucleic Acids kit. We included a 1-hour lysis incubation at 60°C, and cfDNA was eluted in 60 μL water. After extraction, cfDNA was quantified with the Qubit 2.0 Fluorometer and Qubit dsDNA HS Assay Kit, or the Quantus Fluorometer and QuantiFluor ONE dsDNA system. Buffy coat fraction aliquots from centrifuged whole blood samples (Streck) and whole blood aliquots (EDTA) were stored at −80°C prior to DNA extraction. White blood cell (WBC) DNA was extracted from the buffy coat fraction using the QIAGEN DNeasy Blood and Tissue kit or Promega Maxwell RSC Blood DNA kit and Maxwell RSC system, as per the manufacturer's instructions. WBC DNA was quantified with a NanoDrop spectrophotometer following extraction.

Sequencing libraries were prepared from a minimum of 10 ng of plasma DNA input per sample. WBC DNA was sheared into ∼180 bp fragments (with a Covaris focused ultrasonicator or via enzymatic digestion) prior to A-tailing, end repair, Illumina-compatible adapter ligation, and PCR amplification. Library quantification was carried out with a NanoDrop spectrophotometer, and each library was run on an ethidium bromide gel to confirm success. Purified sample libraries were multiplexed to obtain single pools with a combined mass of 1 to 2 μg. For targeted sequencing with the 72-gene or AR locus panels, library pools were hybridized to custom-designed Roche NimbleGen SeqCap EZ Choice capture probe sets. For whole-exome sequencing (WES), library pools were hybridized to the SeqCap EZ MedExome panel. SeqCap EZ HyperCap protocol was followed for hybridization and subsequent wash, recovery, and amplification of the captured regions. Agencourt AMPure beads enabled purification of library pools prior to quantitation with the Qubit 2.0 or Quantus Fluorometer. Pools were diluted and sequenced on an Illumina MiSeq (V3 600 cycle kit) or HiSeq 2500 (V4 250 cycle kit) machine. All deidentified targeted and whole-exome DNA sequencing data have been deposited in the European Genome-phenome Archive (EGA) database under the accession code EGAS00001005318 and are available under standard EGA-controlled release.

Sequence alignment

For both targeted sequencing and WES data, paired-end reads were aligned against an hg38 reference genome using Bowtie-2.3.0 (13). PCR and optical duplicate reads were marked using samblaster-0.1.24 (14) and omitted from analysis. Adapters in read 3′ ends were trimmed in paired-end mode using cutadapt-1.11 (15). Low-quality read tails (smoothed baseq < 30) were trimmed using an in-house algorithm. Per-base read coverages in target regions were quantified using bedtools-2.25.0 (16). cfDNA and WBC sample identities were verified based on SNP genotype profiles.

Identification of mutations

In targeted sequencing data (72-gene panel), somatic mutations (single-nucleotide variants and indels) were called in cfDNA samples by searching for variants with at least 8 supporting reads and a mutant allele fraction (MAF) of at least 0.5%. The allele fraction was required to be at least 20 times higher than the average allele fraction of the same allele in all WBC samples (n = 202), and at least 3 times higher than the allele fraction in the patient's WBC sample. The paired WBC sample had to have at least 20 reads covering the position. The comparisons to patient-matched and cohort-wide WBC samples were necessary to help eliminate (from somatic mutations calls) germline SNPs, stereotypical artifacts arising from sequencing and/or alignment error, as well as mutations of plausible clonal-hematopoietic origin. For base substitutions, the average mapping quality of mutation supporting reads was required to be at least 10, and the average distance of the mutant allele from the nearest read end was required to be at least 15 bases. Our mutation analysis pipeline has been previously validated in matched tissue cfDNA cohorts and in dilution series experiments (11, 17). Protein-level consequences of variants were predicted using ANNOVAR, and all somatic variants were inspected with the Integrative Genomics Viewer for alignment errors.

Somatic mutations were called in WES data using the same pipeline but requiring an allele fraction ≥ 5% to account for the lower sequencing depth. The allele fraction was required to be at least 25 times higher than the average allele fraction across all WBC samples, and at least 10 times higher than the allele fraction in the paired WBC sample.

Germline variants were called in targeted 72-gene panel WBC samples by searching for variants with an alternate allele fraction of at least 15%, and at least 8 supporting reads. Germline variants with population allele frequency 0.5% or higher in the KAVIAR or ExAC databases were discarded. Protein-level consequences of variants were predicted using ANNOVAR. Variant pathogenicity was assessed using ClinVar annotations (18) and based on our previously published protocols (19).

Assessment of mutation redetection rate

Differences in ctDNA fraction and sequencing depth between serial samples can affect sensitivity of mutation detection. To reduce the likelihood of false negatives in comparing mutation presence/absence across samples, mutations independently called in one cfDNA sample were considered not amenable for detection in a later cfDNA timepoint if the expected allele fraction MAFearlier * (Clater / Cearlier) was below 2%, where Clater and Cearlier were ctDNA fractions of the later and earlier samples, respectively. In addition, if the probability of observing eight or more mutant reads in the later sample (based on coverage and expected allele fraction) was less than 90%, the mutation was considered not amenable for detection.

ctDNA fraction estimation

For cfDNA samples sequenced with the 72-gene panel, ctDNA fractions were estimated based on the allele fractions of autosomal somatic mutations in nonamplified regions. MAF in autosomes is elevated when a mutation is accompanied with loss of the other allele (i.e., LOH). In this situation, the MAF and ctDNA fraction are related as ctDNA fraction = 2 / (1 / MAF + 1). Because deletions are not detectable at low ctDNA fractions and LOH is very common in frequently mutated genes such as TP53, we conservatively assumed all mutations to have associated LOH. To account for sampling noise, we modeled the mutant read count as arising from a binomial distribution, and conservatively calculated what the true MAF would be if the highest observed MAF was a 95% quantile outlier. After calculating a ctDNA fraction estimate for each somatic mutation, the highest estimate was adopted as the estimate for the sample, as this mutation was the most likely to be truncal to the metastatic lineage. A diploid tumor model was assumed, because our targeted panel lacks sufficient genome breadth to reliably identify whole genome duplication events. These three assumptions are potential sources of error in calculating ctDNA purity from targeted data. Nonetheless, here and in our prior prostate cancer work (11), we have shown that ctDNA fraction estimates derived from targeted 72-gene panel data are highly correlated with orthogonal approaches based on whole-exome copy-number profiles. Our methodology for approximating ctDNA fraction has also been applied to targeted cfDNA sequencing data from bladder cancer (20).

For samples with WES, ctDNA fractions were estimated based on genome-wide copy-number and heterozygous SNP allele fractions (Supplementary Figs. S1 and S2). Specifically, we divided the genome into segments based on coverage logratio discontinuities identified in any sample from a patient (see Patients and Methods section “Copy-number segmentation”). We then tested combinations of ctDNA fraction and tumor copy-number levels to calculate expected autosomal logratios LR = log2[(F * C + (1 – F) * 2) / 2] (where F = ctDNA fraction and C = number of mutant allele copies), and selected the model that minimized discrepancy between the expected and observed exome-wide autosomal logratios. Candidate models were rejected if any regions had a negative copy number, or if the model suggested a zero copy status for a large genomic region (which is highly biologically improbable). The best model was selected manually to reduce error arising from aneuploidy and potential polyclonality captured by ctDNA; prior experience with automated tools has suggested low accuracy for cfDNA whole-exome data. We corroborated our WES ctDNA fraction estimates using ASCAT (version 2.5.2; ref. 21) and achieved a Pearson correlation of 0.88 between the manual and automated estimates for the 48% of samples where ASCAT was able to find a model.

To model the accuracy of our ctDNA fraction estimates, we compared estimates from the mutation- and copy-number–based procedures using a Bland–Altman plot, and found that a model C = Ctrue * 2N(0, 0.12) + N(0, 0.06) composed of multiplicative noise ∼ 2N(0, 0.12) and additive noise ∼ N(0, 0.06) represented the data well (Supplementary Fig. S3).

Copy-number analysis

Read counts inside panel target regions were calculated using bedtools-2.25.0 (16). A set of 138 reference cfDNA patient samples with no detectable somatic mutations or amplifications from targeted sequencing were used to generate an aggregate reference coverage profile. To do this, we first calculated coverage logratios Lk = log2(Rk,s / Rk,1) for each reference cfDNA sample against an index sample (arbitrarily selected as the first reference sample), where Rk,s represented the read count of capture region k in sample s. We then fitted a Loess regression curve to a two-dimensional scatter plot of (Gk, Lk) pairs—where Gk was the GC fraction of region k—and subtracted the Loess curve from the logratio values to correct for GC bias and differences in global sequencing depth. By inverting the aforementioned formula, we then converted each sample's GC-corrected logratios back to read counts. Finally, we constructed the aggregate reference coverage profile by calculating the median positional coverage of all 138 GC-corrected reference samples. The cfDNA samples used as negative controls were sequenced with the same panel, and to a similar depth.

To calculate coverage logratios for patient cfDNA samples, each sample's coverage profile was compared against the reference profile, and the resulting logratios were GC corrected using Loess regression. A gene's coverage logratio was calculated as the median logratio of all target regions inside the gene.

For samples sequenced with the 72-gene panel, gene deletions and gains were called at coverage logratios of ≤ –0.3 and ≥ 0.3, respectively. These thresholds were determined by examining the distribution of coverage logratios and heterozygous SNP allele fractions in negative control cfDNA samples from healthy volunteers and patients with ctDNA-negative prostate cancer, as previously described (11).

A gene's copy number in the cancer cells contributing to ctDNA was calculated as P * [1 + (2L − 1) / F], where L = the gene's coverage logratio, P = the gene's copy number in normal cells, and F = ctDNA fraction. To derive a CI for a gene's copy number, we assumed that the gene's coverage logratio was normally distributed as L ∼ N(μ, σ2), and used cfDNA samples sequenced with both the 72-gene and whole-exome panel to estimate σ2 = Var(L72GLWES) / 2, where L72G and LWES were coverage logratios measured with the two panels. For simplicity, this calculation assumed that coverage logratios from the two panels were identically distributed.

Copy-number segmentation

To identify temporal copy-number changes between same-patient samples sequenced with the whole-exome or AR panels, we used a segmentation procedure incorporating information from all consecutive samples. Each chromosome was represented as a sequence of coverage logratios L1, …, Ln ordered by genomic position. To identify a breakpoint a, we searched for integers a and b defining a region La, …, Lb such that the statistical difference relative to region L1, …, La-1 was maximized. The difference was defined as Z = (μ1..a − μa..b) / sqrt(1 / a + 1 / (b − a)), where μa..b = mean (La, …, Lb). For each chromosome, the breakpoint position a that produced the highest Z score in any sample from the patient was identified, and the chromosome was split into two segments if the Z score exceeded a specified threshold. This process was repeated recursively for each resulting segment.

To account for variable levels of coverage logratio noise between samples, the Z scores calculated for each sample were divided by a sample noise score N, defined as the average difference between consecutive logratios.

Inference of temporal shifts in genomic populations represented in ctDNA

Somatic genomic profiles from targeted sequencing approaches such as our 72-gene panel do not enable robust characterization of subclonal cancer cell populations, as insufficient somatic mutations are captured. However, a shift in the relative fractions of cancer cell populations between two samples can be inferred even from sparse evidence. Here, we considered a population shift to have occurred between two cfDNA samples sequenced with the 72-gene panel if at least one somatic mutation pair showed a significant change in allele fraction relative to each other, or if at least one gene showed a significant change in its ctDNA fraction adjusted copy number.

To identify a change in the relative allele fraction of two mutations A and B between two cfDNA samples, we started with the null hypothesis p1A / p1B = p2A / p2B, where p1A was the true allele fraction of mutation A in sample 1. This null hypothesis states that the relative allele fraction of mutations A and B did not change between samples. We then ran a million rounds of simulation generating allele fractions s1A ∼ Bin(n1A, o1A) / n1A, s1B ∼ Bin(n1B, o1B) / n1B, s2A ∼ Bin(n2A, cs2A) / n2A, and s2B ∼ Bin(n2B, cs2B) / n2B, where o1A = observed allele fraction of mutation A in sample 1, n1A = observed total read depth at the position of mutation A in sample 1, and c = observed change in average allele fraction between the samples (o2A + o2B) / (o1A + o1B). We then tested how often the observed relative allele fraction o2A / o2B was more extreme than the simulated s2A / s2B. This gave us the probability of the observed change in relative allele fraction occurring due to sampling noise. A population shift was called if this P value was less than 0.00001 for any mutation pair. Samples with ctDNA < 5% were excluded, based on a previous dilution series experiment where mutation detection sensitivity in this cohort approached 100% for samples with ctDNA ≥ 5% (11). If the probability parameter for any binomial distribution would have been zero, it was instead replaced with 1 / read depth, in accordance with Cromwell's rule.

To identify a change in the copy number of gene G between two cfDNA samples, we started with the null hypothesis that the gene had the same copy number (C1 + C2) / 2 in both samples. We then ran a million rounds of simulation where we perturbed the observed ctDNA fractions with noise (Supplementary Fig. S3), then calculated coverage logratios by inverting the C = P * [1 + (2LR − 1) / F] formula, then perturbed the logratios with noise ∼ N(0, S * σG), converted back to copy numbers, and stored the absolute difference of the two simulated copy-number measurements. The amount of noise added to logratios was based on σG = standard deviation of the gene's coverage logratios in control cfDNA samples with zero ctDNA, and S = sample noise multiplier determined by comparing coverage logratios of adjacent genomic regions. A copy-number change was called if less than 1% of simulations (P < 0.01) yielded a difference stronger than the observed difference |C1 - C2|.

In pairs of samples with WES data, we used a different method to identify mutation-based evidence for a shift in cancer cell populations, because the number of tested hypotheses was too large when testing all mutation pairs. We considered a mutation to show evidence for population shift between two WES samples if the ctDNA fraction adjusted MAF was different (Bonferroni adjusted P < 0.00001) between two cfDNA timepoints according to Fisher exact test, and the MAF was higher than 10% in one sample. To adjust for ctDNA fraction in the two samples, the number of mutant reads in each timepoint was divided by C / A, where C = ctDNA fraction of the sample and A = average ctDNA fraction of the two samples. Unlike the method used for the 72-gene panel, this method relies on accurate knowledge of the ctDNA fraction for the two samples.

Identification of structural rearrangements

Somatic structural rearrangements were detected using Breakfast (github.com/annalam/breakfast), with the –merge-duplicates option. Briefly, breakfast extracts 30 bp anchors from the 5′ and 3′ end of unaligned reads and aligns them against the genome. If both anchors align to the genome, an optimal breakpoint inside the original read is sought such that the halves align to their respective anchor regions with the fewest mismatches. Split-read alignments supporting an identical junction structure were grouped together. Somatic rearrangements were required to be supported by at least four unique junction-spanning reads, and to have no supporting reads in any of the 202 WBC samples. All rearrangement junction sequences were inspected using UCSC BLAT, and were discarded if alternative genomic alignments in homologous regions were detected. A rearrangement was classified as a deletion if the 5′ anchor was located upstream of the 3′ anchor in the same strand of the same chromosome, and as a tandem duplication if the 5′ anchor was located downstream of the 3′ anchor. A rearrangement was classified as gene inactivating if it was a deletion or tandem duplication overlapping the gene's coding region, or an inversion or translocation anywhere between the transcription start site and stop codon.

Classification of AR gene structural rearrangements

A rearrangement was classified as an AR gene structural rearrangement (AR-GSR) if at least one breakpoint fell inside the AR gene body (chrX: 67,544,623–67,730,619 in hg38). An AR-GSR was considered LBD truncating if it was (Supplementary Fig. S4):

  • 1) A deletion that fully or partially overlaps the LBD in AR exons 5 to 8, but does not overlap AR exons 1 to 3

  • 2) A tandem duplication that fully overlaps AR exons 1 to 3, but does not fully overlap the LBD in AR exons 5 to 8

  • 3) A tandem duplication that starts downstream of AR exon 3 and ends before the stop codon in exon 8

  • 4) An inversion or translocation with at least one breakpoint located between the end of AR exon 3 and stop codon in exon 8, with a junction flank on the centromeric side containing AR exons 1 to 3

To determine which AR-GSRs showed evidence for an increase in their allele fraction between consecutive cfDNA samples, we calculated the number of unique rearrangement supporting reads R1 and R2 in each sample, and divided them by the median AR sequencing depth A1 and A2 in those samples to account for differences in global sequencing depth between samples.

Longitudinal ctDNA quantity and patient prognosis

In our published trial, we randomized 202 patients with mCRPC to either abiraterone (an inhibitor of androgen synthesis) or enzalutamide (an AR LBD antagonist), followed by crossover to the other AR signaling inhibitor at progression (7). Patients provided serial plasma cfDNA samples before protocol treatment (baseline; representing clinically progressing mCRPC) and at progression on first- and second-line therapy with either inhibitor (Fig. 1A; Supplementary Table S1; Supplementary Fig. S5). In total, 458 samples were collected and subjected to initial targeted sequencing (median read depth 724x) with a panel encompassing 72 prostate cancer genes to estimate ctDNA quantity and to inform on DNA-level temporal changes in driver genes (Supplementary Tables S2 and S3).

Figure 1.

CtDNA fraction consistency across serial clinical progressions. A, Study design. All plasma cfDNA collections across three timepoints were subjected to initial targeted sequencing (72-gene panel). Samples with high ctDNA fraction as a proportion of total cfDNA were additionally subjected to WES (panel) and/or sequencing of all AR introns (AR panel). B, Concordance of somatic MAFs detected by the targeted (72-gene) and exome panels. Each dot represents a cfDNA sample sequenced with both approaches (134 samples in total). C, Validation of ctDNA fraction estimates from targeted sequencing (72-gene panel), based on WES copy-number profiles. D, Correlation of ctDNA fractions in same-patient consecutive cfDNA sample pairs, shown as a mirrored barplot. Color indicates the number of days between the sample collections. Only samples collected after progression are shown (i.e., 21 samples improperly collected during effective treatment were excluded). E, Relative ctDNA fraction change between consecutive cfDNA collection timepoints is associated with relative serum PSA change between the timepoints. Boxes represent increasing quartiles of relative PSA change. Consecutive timepoints where neither sample had detectable ctDNA were excluded. See Supplementary Fig. S9 for more detail. F, Kaplan–Meier curve showing the impact of baseline ctDNA fraction on overall survival. G and H, High ctDNA fraction at first progression and short first-line TTP are associated with short TTP on second-line enzalutamide treatment.

Figure 1.

CtDNA fraction consistency across serial clinical progressions. A, Study design. All plasma cfDNA collections across three timepoints were subjected to initial targeted sequencing (72-gene panel). Samples with high ctDNA fraction as a proportion of total cfDNA were additionally subjected to WES (panel) and/or sequencing of all AR introns (AR panel). B, Concordance of somatic MAFs detected by the targeted (72-gene) and exome panels. Each dot represents a cfDNA sample sequenced with both approaches (134 samples in total). C, Validation of ctDNA fraction estimates from targeted sequencing (72-gene panel), based on WES copy-number profiles. D, Correlation of ctDNA fractions in same-patient consecutive cfDNA sample pairs, shown as a mirrored barplot. Color indicates the number of days between the sample collections. Only samples collected after progression are shown (i.e., 21 samples improperly collected during effective treatment were excluded). E, Relative ctDNA fraction change between consecutive cfDNA collection timepoints is associated with relative serum PSA change between the timepoints. Boxes represent increasing quartiles of relative PSA change. Consecutive timepoints where neither sample had detectable ctDNA were excluded. See Supplementary Fig. S9 for more detail. F, Kaplan–Meier curve showing the impact of baseline ctDNA fraction on overall survival. G and H, High ctDNA fraction at first progression and short first-line TTP are associated with short TTP on second-line enzalutamide treatment.

Close modal

Somatic mutations with an allele fraction above 1% were detected in 119 of 202 (59%) baseline samples, 82 of 158 (52%) first progression samples, and 67 of 98 (68%) second progression samples (Supplementary Table S4). Using the highest allele fraction mutation in nonamplified regions, estimated ctDNA fractions (as a proportion of total cfDNA) ranged between undetectable and 89%, and were strongly correlated with clinical markers of tumor burden (Supplementary Table S5; Supplementary Fig. S6). Select samples (mainly those with high ctDNA fractions, but including some low ctDNA fraction controls) were additionally subjected to targeted sequencing of the entire AR gene and flanking regions (n = 175; median unique read depth 1,093x) and/or WES (n = 134; median depth 144x; Fig. 1A; Supplementary Table S3). Technical concordance was high between targeted sequencing and WES (Fig. 1B; Supplementary Fig. S7). ctDNA fractions estimated based on whole-exome somatic profiles were consistent with mutation-based estimates from targeted sequencing (Fig. 1C; Supplementary Fig. S8).

ctDNA fractions between consecutive cfDNA collections were correlated (Fig. 1D). The correlation was stronger when timepoints were separated by < 180 days (the median collection interval; Pearson r = 0.75 vs. 0.49). Only 21 of 63 (33%) patients with multiple samples and ctDNA fraction < 2% at baseline had ctDNA ≥ 2% in a subsequent progression sample. Changes in ctDNA fraction between collections were correlated with changes in serum PSA, a clinical surrogate for depth of therapy response (ref. 22; Fig. 1E; Supplementary Fig. S9). Eleven cfDNA samples were collected too early (before progression) due to clinical circumstances, and expectedly had lower ctDNA fractions than all other samples (median undetectable vs. 3.9%, P = 0.008, rank-sum test). Collectively, these results suggest that repeat ctDNA testing for genomic characterization may have limited utility in patients with consistently low tumor burden.

Patients with baseline ctDNA ≥ 5% (the cohort median) exhibited shorter overall survival (median 16.6 months vs. not reached; HR, 4.57; 95% CI, 3.01–6.92; P < 0.001) and time from start of first-line therapy to second progression (median 7.8 vs. 19.1 months; HR, 2.63; 95% CI, 1.88–3.67; P < 0.001) than patients with low or undetectable ctDNA. This association was quantitative, with higher ctDNA fractions associated with progressively worse overall survival and TTP (Fig. 1F; Supplementary Fig. S10). Importantly, baseline ctDNA fraction achieved higher prognostic accuracy than a multivariate model combining all routine clinical prognostic factors (Supplementary Fig. S11). In our trial, almost no patients benefited from second-line abiraterone, but 51% of patients experienced a PSA response of 30% or greater to second-line enzalutamide (Supplementary Fig. S5; ref. 7). Two characteristics showed the strongest association with TTP on second-line enzalutamide: ctDNA fraction at crossover and TTP on first-line abiraterone (Fig. 1G and H; Supplementary Table S6).

Consistent detection of truncal driver alterations across serial clinical progressions

To compare driver gene alteration status over time, we focused on 164 same-patient sample pairs (across 66 patients) with sufficient ctDNA for robust genotyping (ctDNA ≥ 5%; see Patients and Methods for justification). In these pairs, ctDNA fraction adjusted MAFs and copy numbers were strongly correlated (Fig. 2A and B). We observed high temporal concordance for the copy number, rearrangement, and mutation status of key prostate cancer driver genes associated with tumor initiation and development, including SPOP, FOXA1, PTEN, and TP53 (Fig. 2C and D; Supplementary Table S7). Across all regions captured via targeted sequencing, 78% of somatic mutations detected in an earlier sample were independently redetected in the subsequent sample. The rate of redetection increased to 91% after excluding discordances explained by insufficient sequencing depth and/or reduced ctDNA fraction in the later sample (Supplementary Fig. S12). Mutations estimated to be present in over 50% of ctDNA-shedding cells were redetected more often than less dominant mutations (93% vs. 76%, P < 0.0001, Fisher exact test). Together, these data demonstrate that pretreatment and postprogression ctDNA populations share a large common reservoir of cancer driver alterations.

Figure 2.

Consistent detection of truncal driver alterations across timepoints. A, Relationship between somatic MAFs in earlier versus later cfDNA samples, highlighting that most mutations are present in similar abundance relative to total ctDNA fraction. B, Correlation of gene copy numbers in same-patient cfDNA sample pairs. Each dot represents the copy number of one gene in one sample pair. Only sample pairs where both samples had at least 30% ctDNA are shown. C, Correlation of the copy-number status of key driver genes between consecutive cfDNA samples from the same patient. Each dot represents a consecutive sample pair. D, Allele fraction of somatic mutations in key driver genes in consecutive cfDNA samples collected at clinical progressions. Each cluster of bars represents consecutive cfDNA samples from one patient. Patient identifiers are shown underneath. To adjust for changes in overall ctDNA fraction, allele fractions were divided by ctDNA fraction/2. Error bars indicate 95% CI based on binomial distribution.

Figure 2.

Consistent detection of truncal driver alterations across timepoints. A, Relationship between somatic MAFs in earlier versus later cfDNA samples, highlighting that most mutations are present in similar abundance relative to total ctDNA fraction. B, Correlation of gene copy numbers in same-patient cfDNA sample pairs. Each dot represents the copy number of one gene in one sample pair. Only sample pairs where both samples had at least 30% ctDNA are shown. C, Correlation of the copy-number status of key driver genes between consecutive cfDNA samples from the same patient. Each dot represents a consecutive sample pair. D, Allele fraction of somatic mutations in key driver genes in consecutive cfDNA samples collected at clinical progressions. Each cluster of bars represents consecutive cfDNA samples from one patient. Patient identifiers are shown underneath. To adjust for changes in overall ctDNA fraction, allele fractions were divided by ctDNA fraction/2. Error bars indicate 95% CI based on binomial distribution.

Close modal

Evidence for shift in tumor populations after treatment

Because rapid autopsy studies have revealed interlesion genomic heterogeneity in mCRPC (23), we hypothesized that treatment may drive temporal changes in the proportions of genomically distinct cancer cell populations represented in ctDNA. Therefore, we searched for quantitative differences in somatic MAFs and gene copy numbers in the 164 sample pairs with ctDNA ≥ 5% (see Patients and Methods). Using targeted sequencing alone (72-gene panel), 36% of patients (24/66) displayed evidence for a temporal shift in detected cancer cell populations (Supplementary Table S8). MAF profiles changed significantly in 21 patients (Fig. 3A; Supplementary Table S9), and copy-number profiles shifted in 17 patients (with an overlap of 12 patients; Fig. 3B, Supplementary Fig. S13).

Figure 3.

Evidence for temporal shifts in cancer cell populations represented in ctDNA. A, Temporal shifts in driver gene mutation profiles, indicated by changes in somatic MAFs between timepoints (targeted sequencing; 72-gene panel). All 21 patients with significant changes are shown (see Supplementary Table S9). Bars represent somatic mutations, are grouped by patient, and are shown in the same order for the two timepoints. Note that allele fractions are not adjusted for ctDNA fraction in this plot. Only patients with ctDNA ≥ 5% in at least two samples were evaluated. B, Changes in gene copy number between timepoints (targeted sequencing; 72-gene panel). Five representative cases are shown, from a total of 17 patients where temporal copy-number differences were found (see Supplementary Fig. S13 for all cases). Horizontal and vertical error bars indicate 80% CIs. Genes showing a significant copy-number difference between timepoints are highlighted in red. C, Validation of shifts in cancer cell populations based on exome-wide somatic mutation profiles from serial cfDNA collections. Three representative cases are shown (see Supplementary Figs. S14 and S17 for all cases). MAFs are shown as dots. Mutations indicating population shift between timepoints are indicated with red ellipses. Protein-altering mutations are shown in black, and silent mutations are shown in gray. D, Somatic MAFs as detected in three serial timepoints from patient 180 (harboring a germline BRCA2 mutation and LOH). Although the timepoints share a small set of apparently truncal mutations, most somatic mutations change in relative prevalence over time, supporting a marked shift in which population B arises during first-line treatment and population A diminishes during both treatments. E, Validation of population shifts based on whole-exome copy-number profiles. Two example cases are shown (see Supplementary Figs. S15 and S18 for all cases). Genomic regions with temporal copy-number differences are indicated with red arrows.

Figure 3.

Evidence for temporal shifts in cancer cell populations represented in ctDNA. A, Temporal shifts in driver gene mutation profiles, indicated by changes in somatic MAFs between timepoints (targeted sequencing; 72-gene panel). All 21 patients with significant changes are shown (see Supplementary Table S9). Bars represent somatic mutations, are grouped by patient, and are shown in the same order for the two timepoints. Note that allele fractions are not adjusted for ctDNA fraction in this plot. Only patients with ctDNA ≥ 5% in at least two samples were evaluated. B, Changes in gene copy number between timepoints (targeted sequencing; 72-gene panel). Five representative cases are shown, from a total of 17 patients where temporal copy-number differences were found (see Supplementary Fig. S13 for all cases). Horizontal and vertical error bars indicate 80% CIs. Genes showing a significant copy-number difference between timepoints are highlighted in red. C, Validation of shifts in cancer cell populations based on exome-wide somatic mutation profiles from serial cfDNA collections. Three representative cases are shown (see Supplementary Figs. S14 and S17 for all cases). MAFs are shown as dots. Mutations indicating population shift between timepoints are indicated with red ellipses. Protein-altering mutations are shown in black, and silent mutations are shown in gray. D, Somatic MAFs as detected in three serial timepoints from patient 180 (harboring a germline BRCA2 mutation and LOH). Although the timepoints share a small set of apparently truncal mutations, most somatic mutations change in relative prevalence over time, supporting a marked shift in which population B arises during first-line treatment and population A diminishes during both treatments. E, Validation of population shifts based on whole-exome copy-number profiles. Two example cases are shown (see Supplementary Figs. S15 and S18 for all cases). Genomic regions with temporal copy-number differences are indicated with red arrows.

Close modal

To validate across a broader genomic context, we analyzed serial cfDNA WES data from 36 patients, focusing on patient-matched sample pairs where high ctDNA fraction (≥30%) enabled more sensitive comparisons. Among patients with evidence for population shifts from initial targeted sequencing, exome-wide changes in mutation and copy-number profiles were confirmed in all evaluable patients (Fig. 3C–E; Supplementary Figs. S14 and S15). In one extreme case (patient 180), we observed an almost complete switch in ctDNA populations between baseline and second progression, and a mixture of the two distantly related populations in the first progression sample (Fig. 3D). Whole-exome analysis supported the absence of ctDNA changes in patients classified as temporally similar based on targeted sequencing data (Supplementary Figs. S16–S18). One patient (026) was excluded from analysis because his mCRPC baseline and “progression” samples shared no mutations, and he was diagnosed with concurrent metastatic bladder cancer (Supplementary Fig. S19).

With the exception of the AR (addressed below), the only recurrent posttreatment changes detected by our 72-gene panel were WNT pathway changes (n = 3 patients), PI3K pathway alterations (n = 3), and cell-cycle gene changes (n = 4, including two patients that acquired a CCND1 amplification at progression; Supplementary Fig. S20).

Dynamic evolution of AR copy number and ligand affinity

Missense mutations in the AR LBD result in differential ligand affinity, and prior work has suggested genotype switching in response to enzalutamide or abiraterone (24, 25). Accordingly, we observed high variation between collection timepoints for somatic AR LBD mutations (Supplementary Fig. S21). We detected a marked increase in AR L702H mutations previously linked to abiraterone resistance through glucocorticoid agonism (glucocorticoids are administered concurrently with abiraterone; Fig. 4A; refs. 24, 26). Interestingly, although 6 patients acquired AR L702H mutations after exposure to abiraterone plus glucocorticoids, 2 patients (035 and 163) acquired L702H mutations after enzalutamide, without accompanying glucocorticoids (Supplementary Fig. S22). Both patients harbored additional AR LBD mutations, and patient 163 harbored a mismatch repair defect and hypermutation (27). The AR F877 L mutation linked with enzalutamide resistance in model systems was not identified in any samples (28, 29).

Figure 4.

Treatment-driven selection of aggressive AR gene alterations. A, Frequency of AR LBD mutations at three cfDNA collection timepoints (1 = baseline, 2 = first progression, 3 = second progression). Mutations are grouped by affected amino acid (codon). B, Change in AR gene copy number between consecutive cfDNA samples. Copy number was defined as the average number of AR copies in cancer cells that released ctDNA into the blood (see Patients and Methods). Sample pairs showing a statistically significant increase in AR copy number are highlighted in red, and those showing a decrease are highlighted in blue. Only samples with ctDNA fraction ≥ 5% were included in analysis. Error bars indicate 80% CIs. C, Box plot showing distribution of AR copy numbers across the three timepoints. Red horizontal lines indicate median. Red diamonds indicate mean. D, Association between AR gene copy number in baseline ctDNA and patient overall survival. Patients where AR copy number could not be assessed due to low ctDNA (<5%) are plotted separately. E, Average copy-number change between consecutive cfDNA samples in a 24 Mb region surrounding AR, based on 30 kb resolution targeted sequencing of AR and its flanking regions. F, Change in AR amplicon structure in three patients that showed a change in AR copy number between cfDNA timepoints. Patient 162 had an estimated ctDNA fraction of 17% at progression (based on 7 somatic mutations that were shared with the baseline sample), but showed no evidence for the AR-amplified cancer cell population detected at baseline. Amplicon structure changes were observed in 17 patients in total (see Supplementary Fig. S23). G, Correlation between coverage logratios measured inside the AR gene and at the AR enhancer, in all 175 cfDNA samples (black dots) sequenced with the AR panel. Some samples where the AR and enhancer copy number differs are highlighted with colored ellipses with a labeled patient identifier. Ellipses containing multiple samples represent samples from the same patient.

Figure 4.

Treatment-driven selection of aggressive AR gene alterations. A, Frequency of AR LBD mutations at three cfDNA collection timepoints (1 = baseline, 2 = first progression, 3 = second progression). Mutations are grouped by affected amino acid (codon). B, Change in AR gene copy number between consecutive cfDNA samples. Copy number was defined as the average number of AR copies in cancer cells that released ctDNA into the blood (see Patients and Methods). Sample pairs showing a statistically significant increase in AR copy number are highlighted in red, and those showing a decrease are highlighted in blue. Only samples with ctDNA fraction ≥ 5% were included in analysis. Error bars indicate 80% CIs. C, Box plot showing distribution of AR copy numbers across the three timepoints. Red horizontal lines indicate median. Red diamonds indicate mean. D, Association between AR gene copy number in baseline ctDNA and patient overall survival. Patients where AR copy number could not be assessed due to low ctDNA (<5%) are plotted separately. E, Average copy-number change between consecutive cfDNA samples in a 24 Mb region surrounding AR, based on 30 kb resolution targeted sequencing of AR and its flanking regions. F, Change in AR amplicon structure in three patients that showed a change in AR copy number between cfDNA timepoints. Patient 162 had an estimated ctDNA fraction of 17% at progression (based on 7 somatic mutations that were shared with the baseline sample), but showed no evidence for the AR-amplified cancer cell population detected at baseline. Amplicon structure changes were observed in 17 patients in total (see Supplementary Fig. S23). G, Correlation between coverage logratios measured inside the AR gene and at the AR enhancer, in all 175 cfDNA samples (black dots) sequenced with the AR panel. Some samples where the AR and enhancer copy number differs are highlighted with colored ellipses with a labeled patient identifier. Ellipses containing multiple samples represent samples from the same patient.

Close modal

Mutations in AR amino acids 742 to 743 were found in 6 patients, all of whom had previously received bicalutamide, a first-generation AR signaling inhibitor. These mutations are a known resistance mechanism to bicalutamide, and showed decreased frequency after enzalutamide and/or abiraterone treatment (Fig. 4A), consistent with preclinical evidence of enzalutamide activity against these mutants (30).

Somatic AR amplifications were detected in 62 of 119 (52%) patients at baseline and 74 of 149 (50%) progression samples (excluding samples with undetectable ctDNA; Supplementary Fig. S21). Given the switch in AR mutation genotypes, we hypothesized that AR amplifications detected before and after treatment may constitute different events in some patients. Nine of 66 patients with sufficient ctDNA in two or more samples to quantify AR copy number displayed a significant temporal increase in AR copy number per average ctDNA-shedding tumor cell (Fig. 4B). Three patients displayed a significant decrease in AR copy number. Consistent with these observations, across all evaluable patients, mean AR copy number increased between the three collection timepoints (7.8 vs. 10.1 vs. 12.3 copies; Fig. 4C). Because an increased AR copy number can result in elevated AR protein expression and greater ligand sensitivity (31), populations with higher AR copy number may have greater fitness during treatment. Accordingly, higher baseline AR copy number was associated with shorter overall survival, even in multivariate models incorporating ctDNA fraction (Fig. 4D; Supplementary Table S6).

To explore the genomic context of AR copy-number increase, we analyzed a 30 Mbp window around AR in 175 cfDNA samples (see Patients and Methods). Consistent with our initial observation of temporally increasing AR gene copy number, a spatial profile of average chromosome X copy-number change between consecutive timepoints revealed a focal amplification centered at the AR and its upstream region (Fig. 4E). All evaluable patients with AR copy-number change also displayed a change in chromosome X amplicon structure (Fig. 4F; Supplementary Fig. S23). The observed changes in amplicon boundaries suggest that AR copy-number increases are due to treatment-induced selection for a distinct cell population with a different pattern of AR gene copy gain. Interestingly, an additional 9 patients displayed a change in AR amplicon structure without significant change in bulk AR copy number, consistent with the notion that patients harbor multiple similar AR locus genotypes.

An enhancer region 640 kb upstream of AR was recently reported to contribute to increased AR expression through focal tandem duplications in patients with mCRPC (32, 33). In our cohort, most AR amplifications also encompassed the enhancer (Fig. 4G). Indeed, 15 patients exhibited AR locus copy gains where the amplification peak began at the AR enhancer (Supplementary Fig. S24). Only 7 patients showed an AR amplification that excluded the enhancer. Fourteen patients showed evidence for a focal amplification of the AR enhancer. Together, these data suggest that potent AR signaling inhibition precipitates population restructuring that strongly favors lineages with high-level AR gene and enhancer amplification.

Emergence of AR truncating genomic rearrangements

AR proteins with a truncated LBD can facilitate target gene transcription independent of ligand availability, and can drive resistance to approved AR signaling inhibitors in preclinical models (34). LBD-truncated AR transcripts can be produced via aberrant splicing, but recent work suggests that genomic breakpoints falling within the AR intragenic region can also give rise to truncated transcripts (3, 35). Here, we identified AR-GSRs in 69 of 175 (39%) cfDNA samples subjected to targeted sequencing of the entire AR locus (Supplementary Table S10). All detected AR-GSRs were unique to a single patient. Of the 170 distinct AR-GSRs identified across all timepoints, 36% were intragenic, 50% juxtaposed AR with another chromosome X region, and 14% were interchromosomal. Copy-number profiles inside AR were consistent with detected AR-GSR breakpoints: among 8 patients with clear evidence for an intragenic AR copy-number change, a concordant AR-GSR was identified in all except one (077, where the expected breakpoint was inside a LINE element and therefore unresolvable; Supplementary Fig. S25). These data are consistent with the AR being one of the most frequently rearranged genes in mCRPC (4, 35).

Based on structural analysis (see Patients and Methods; Supplementary Fig. S4), 34 unique AR-GSRs had potential to produce a truncated AR isoform lacking part or all of the LBD (exons 5–8) while preserving the DNA binding domain (exons 1–3; Supplementary Table S10). The number of detected AR-GSRs was strongly correlated with AR copy number (Fig. 5A). cfDNA samples with amplified AR were 9 times more likely to carry AR-GSRs (61% vs. 7%, P < 0.001, Fisher exact test), including LBD truncating rearrangements (26% vs. 3%, P < 0.001, Fisher exact test).

Figure 5.

Evolution of AR rearrangements during abiraterone and enzalutamide treatment. A, Number of detected AR-GSRs in cfDNA samples, ordered by number of AR gene copies. Rearrangements predicted to truncate the AR LBD are shown in red color, and other AR-GSRs in gray. Inset plot shows the average number of AR-GSRs in samples grouped by AR copy number. B, Percentage of patients carrying LBD truncating (red) and other AR-GSRs (gray) at the three cfDNA timepoints (baseline, first progression, and second progression). C, Relative change in the number of rearrangement supporting reads (divided by the median AR gene coverage) between two consecutive cfDNA timepoints, grouped by type of rearrangement. Values above 100 indicate increases from zero, and values below 0.001 indicate decreases to zero. D, Average sequencing depth (gray silhouette) and locations of rearrangement breakpoints (black dots) across the entire AR region (hg38 coordinates chrX: 67,544,032–67,730,619) in the 175 cfDNA samples sequenced with the AR panel. E, Structure of LBD truncating AR rearrangements (left side), and the number of unique reads supporting each rearrangement at the three cfDNA timepoints (right side).

Figure 5.

Evolution of AR rearrangements during abiraterone and enzalutamide treatment. A, Number of detected AR-GSRs in cfDNA samples, ordered by number of AR gene copies. Rearrangements predicted to truncate the AR LBD are shown in red color, and other AR-GSRs in gray. Inset plot shows the average number of AR-GSRs in samples grouped by AR copy number. B, Percentage of patients carrying LBD truncating (red) and other AR-GSRs (gray) at the three cfDNA timepoints (baseline, first progression, and second progression). C, Relative change in the number of rearrangement supporting reads (divided by the median AR gene coverage) between two consecutive cfDNA timepoints, grouped by type of rearrangement. Values above 100 indicate increases from zero, and values below 0.001 indicate decreases to zero. D, Average sequencing depth (gray silhouette) and locations of rearrangement breakpoints (black dots) across the entire AR region (hg38 coordinates chrX: 67,544,032–67,730,619) in the 175 cfDNA samples sequenced with the AR panel. E, Structure of LBD truncating AR rearrangements (left side), and the number of unique reads supporting each rearrangement at the three cfDNA timepoints (right side).

Close modal

Recent studies have linked LBD truncating AR-GSRs to de novo treatment resistance (11, 35, 36), but it is unknown whether AR signaling inhibitors actively select for these genotypes in responding patients. Here, the fraction of patients carrying AR LBD truncating rearrangements increased successively with progression on AR signaling inhibitors (11% at baseline vs. 20% at first progression vs. 25% at second progression; Fig. 5B). Comparing coverage-normalized read counts at consecutive timepoints, we observed that although 72% of all AR-GSRs increased in allele fraction over time, LBD truncating rearrangements were under stronger positive selection than other AR-GSRs (P = 0.003, rank-sum test; Fig. 5C). Eight patients displayed new LBD truncating rearrangements in progression samples, with no supporting evidence at baseline despite a median ctDNA fraction of 41% and median sequencing depth of 1,388x at the breakpoints. Conversely, no LBD truncating rearrangements decreased to undetectable levels over time, compared with 13% of other AR-GSRs (P = 0.0036, Fisher exact test; Fig. 5C).

Intriguingly, we observed an enrichment for intragenic deletions removing AR exons 5 to 7 while preserving the final exon 8 (Fig. 5D and E). In patients 071 and 110, multiple independent deletions matching this pattern were detected solely in progression samples. Because deletions of exons 5 to 7 juxtapose frame-incompatible exons 4 and 8 (which if spliced together would yield a frameshifted AR isoform), their apparent positive selection may suggest a benefit for preservation of the 3′-untranslated region (and polyadenylation site) or adjacent regulons.

Genomic change is associated with treatment response

Combining results from all analysis approaches, treatment-induced genomic changes were detected in 53% of evaluable patients (Supplementary Table S8). The vast majority (86%) of these patients exhibited AR locus changes, suggesting that differential selective pressure on distinct AR genotypes may be driving the overall population shifts (Fig. 6A). In support of this concept, 7 patients exhibited temporal changes only at the AR locus. Patients with primary resistance against first-line treatment (TTP < 3 months) did not show systematic increases in AR copy number (Fig. 6B), suggesting that most changes were occurring in the context of disease that was initially sensitive to treatment. Consistent with this, patients with evidence for temporal changes in ctDNA populations were more likely to have experienced a deep PSA response to intervening therapy than other patients (median PSA response 76% vs. 28%, P = 0.028; Fig. 6C). Indeed, all 8 patients with greater than 90% PSA decline and ctDNA fraction ≥ 10% in serial cfDNA collections showed population shifts.

Figure 6.

Relationship of treatment response to temporal shifts in somatic populations and AR genotype. A, Oncoprint summarizing the AR gene locus changes during treatment in 66 evaluable patients (those with pairs of cfDNA samples with ctDNA fraction ≥ 5%). Six different types of somatic AR alteration are shown. Red color indicates an increase in the number of AR copies, AR mutations, or AR LBD truncating rearrangements after treatment. Blue color indicates a respective decrease. B,AR copy-number change between baseline and first progression cfDNA collection timepoints in patients with a sustained response to first-line treatment (TTP > 6 months), and patients with evidence for primary resistance (TTP < 3 months). Only patients with ctDNA ≥ 5% at both timepoints are shown. C, Best PSA decline during treatment in 66 patients, grouped according to whether a temporal shift in cancer cell populations was observed. Only treatments received between cfDNA collections with ctDNA fraction ≥ 5% were considered when determining best PSA decline.

Figure 6.

Relationship of treatment response to temporal shifts in somatic populations and AR genotype. A, Oncoprint summarizing the AR gene locus changes during treatment in 66 evaluable patients (those with pairs of cfDNA samples with ctDNA fraction ≥ 5%). Six different types of somatic AR alteration are shown. Red color indicates an increase in the number of AR copies, AR mutations, or AR LBD truncating rearrangements after treatment. Blue color indicates a respective decrease. B,AR copy-number change between baseline and first progression cfDNA collection timepoints in patients with a sustained response to first-line treatment (TTP > 6 months), and patients with evidence for primary resistance (TTP < 3 months). Only patients with ctDNA ≥ 5% at both timepoints are shown. C, Best PSA decline during treatment in 66 patients, grouped according to whether a temporal shift in cancer cell populations was observed. Only treatments received between cfDNA collections with ctDNA fraction ≥ 5% were considered when determining best PSA decline.

Close modal

This study provides insights into mCRPC genomic evolution in living patients as their disease responds and progresses through standard-of-care treatment with abiraterone and enzalutamide. Importantly, observed shifts in the predominant plasma ctDNA profile and AR genotype over sequential clinical progressions propose a Darwinian explanation for acquired resistance to these widely used AR signaling inhibitors.

Somatic profiles changed over time in most patients where AR signaling inhibitors elicited strong responses. Differences in the relative abundance of mutations, copy numbers, and rearrangements as a proportion of total ctDNA supported a model in which the metastatic niche harbors several phylogenetically related competing populations. This is consistent with genomic analyses of metastatic lesions obtained at autopsy, which show interlesion heterogeneity at the time of patient death (23).

Very few temporal changes were detected in TP53, FOXA1, or SPOP, supportive of prior work indicating that somatic alterations in these genes are early events truncal to the metastatic lineage (23, 37, 38). However, we observed a striking convergence of genomic changes upon AR. Compared with baseline mCRPC, the temporal changes detected here in AR after progression on abiraterone and/or enzalutamide reflect a shift toward even higher levels of AR augmentation. Although the concept of an AR genotype that changes during mCRPC treatment is well accepted for AR LBD mutations (39), previous efforts to understand the clinical relevance of AR amplification have not explored whether amplifications detected at different timepoints represent the same event (24, 25, 40). Here, the mean AR copy number in ctDNA increased by 50% after progression on abiraterone and enzalutamide. Furthermore, drastic temporal differences in AR locus structure within individual patients suggested that stronger AR amplification after AR-targeted therapy often arose through selection of a different tumor population relative to baseline mCRPC. This treatment-associated pressure toward increased AR copy number may help explain why most mCRPC with AR amplification initially respond to potent AR inhibitors, and is consistent with observations that AR amplification associates with poor prognosis (41, 42).

Our work promotes LBD truncating AR-GSRs as a major driver of clinical resistance to AR signaling inhibitors. Robust in vitro work has underscored the potential of these genomic variants to drive resistant phenotypes via expression of constitutively active AR protein isoforms (3, 35), but the alternative hypothesis that AR-GSRs represent a neutral side effect of AR amplification has never been ruled out. Here using serial sampling in a live patient population, we show that the frequency of LBD truncating AR-GSRs increases 2.5-fold after treatment with abiraterone/enzalutamide, and that other AR rearrangements are under weaker selection. Combined with mutation and copy-number observations, the consistent change toward more active AR paints a picture of powerful selective pressure toward hyperaugmented AR and emergence of lethal cross-resistance.

The emergence of somatic populations with heightened AR augmentation supports the concept that a large proportion of resistant disease retains some reliance upon the AR, even after potent pathway suppression. As such, renewed targeting of the AR signaling pathway may still have an antitumor effect, if therapeutic strategies are independent of the AR LBD. Nevertheless, it may be more beneficial to temporarily reduce the selective pressure on the AR. This concept is supported by a recent phase III clinical trial demonstrating that second-line taxane-based chemotherapy outperforms treatment with sequential AR signaling inhibitors (6). Correlative studies have also suggested that patients with AR copy gain in plasma cfDNA may respond better to chemotherapy than AR signaling inhibitors (43), supporting the general hypothesis that chemotherapy is less dependent on AR genotype.

The overall abundance of ctDNA is increasingly recognized as an impactful prognostic factor in metastatic cancers (44). We and others have reported on the potential for ctDNA fraction to add to clinical prognostic models in mCRPC (11, 45). Although the current study was not focused on biomarker discovery, we note that ctDNA fraction may help identify patients likely to be responsive to second-line enzalutamide. This concept will be tested in a biomarker-driven phase II trial (NCT04015622).

Our analysis of treatment-driven genomic changes was limited to patients with intermediate outcomes, as sequential collections were not available from patients whose disease responded exceptionally badly (i.e., rapid death) or exceptionally well (i.e., no progression). Among patients with sequential collections, we adjusted for temporal changes in ctDNA fraction, but were forced to exclude samples with low ctDNA fraction, introducing a potential source of bias. Furthermore, this study focused on known driver genes and coding regions, limiting our ability to identify novel changes in the noncoding genome (46), or comprehensively resolve and characterize individual subclonal populations. Future studies should utilize whole-genome sequencing on samples with high ctDNA fractions to inform on clonal structure.

M. Annala reports grants from Jane and Aatos Erkko Foundation and Academy of Finland during the conduct of the study. S. Taavitsainen reports grants from Jane and Aatos Erkko Foundation during the conduct of the study. D.J. Khalaf reports personal fees from Janssen outside the submitted work. D.L. Finch reports other support from Janssen, Astellas, and Bayer outside the submitted work. J. Vergidis reports personal fees from BMS and non-financial support from Bayer outside the submitted work. B.J. Eigl reports personal fees from AstraZeneca, Pfizer, Seagen, Johnson and Johnson, EMD Serono, and Merck outside the submitted work. C.K. Kollmansberger reports personal fees from Astellas, Janssen, Pfizer, BMS, Eisai, Merck, Ipsen, and Bayer outside the submitted work. M. Nykter reports grants from Jane and Aatos Erkko Foundation and Academy of Finland during the conduct of the study. K.N. Chi reports grants from Astellas and Janssen during the conduct of the study. K.N. Chi also reports grants and personal fees from AstraZeneca, Bayer, Novartis, Pfizer, Point Biopharma, Roche, and Sanofi, as well as personal fees from Daiichi Sankyo, Merck, and Bristol Myers Squibb outside the submitted work. A.W. Wyatt reports personal fees from AstraZeneca, Merck, and Astellas; grants and personal fees from Janssen; and grants from ESSA Pharma outside the submitted work. No disclosures were reported by the other authors.

M. Annala: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. S. Taavitsainen: Data curation, formal analysis, investigation, methodology. D.J. Khalaf: Data curation, formal analysis, investigation. G. Vandekerkhove: Investigation, methodology, writing–review and editing. K. Beja: Investigation, methodology. J. Sipola: Data curation, software, methodology. E.W. Warner: Investigation, methodology. C. Herberts: Investigation, methodology, writing–review and editing. A. Wong: Investigation. S. Fu: Investigation. D.L. Finch: Investigation. C.D. Oja: Investigation. J. Vergidis: Investigation. M. Zulfiqar: Investigation. B.J. Eigl: Investigation. C.K. Kollmansberger: Investigation. M. Nykter: Supervision, project administration. M.E. Gleave: Investigation, project administration. K.N. Chi: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. A.W. Wyatt: Conceptualization, resources, data curation, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing.

This work was supported by the Canadian Institutes of Health Research, Canadian Cancer Society Research Institute, Prostate Cancer Canada, Movember Foundation, Prostate Cancer Foundation, Jane and Aatos Erkko Foundation, Academy of Finland, Terry Fox New Frontiers Program, BC Cancer Foundation, Janssen, and Astellas. No funding sources were involved in the design or execution of the study. The authors are thankful to all contributing patients and their families.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Attard
G
,
Parker
C
,
Eeles
RA
,
Schröder
F
,
Tomlins
SA
,
Tannock
I
, et al
Prostate cancer
.
Lancet
2016
;
387
:
70
82
.
2.
Abida
W
,
Cyrta
J
,
Heller
G
,
Prandi
D
,
Armenia
J
,
Coleman
I
, et al
Genomic correlates of clinical outcome in advanced prostate cancer
.
Proc Natl Acad Sci U S A
2019
;
116
:
11428
36
.
3.
Henzler
C
,
Li
Y
,
Yang
R
,
McBride
T
,
Ho
Y
,
Sprenger
C
, et al
Truncation and constitutive activation of the androgen receptor by diverse genomic rearrangements in prostate cancer
.
Nat Commun
2016
;
7
:
13668
.
4.
Quigley
DA
,
Dang
HX
,
Zhao
SG
,
Lloyd
P
,
Aggarwal
R
,
Alumkal
JJ
, et al
Genomic hallmarks and structural variation in metastatic prostate cancer
.
Cell
2018
;
175
:
889
.
5.
Sartor
O
,
de Bono
JS
. 
Metastatic prostate cancer
.
N Engl J Med
2018
;
378
:
645
57
.
6.
de Wit
R
,
de Bono
J
,
Sternberg
CN
,
Fizazi
K
,
Tombal
B
,
Wülfing
C
, et al
Cabazitaxel versus abiraterone or enzalutamide in metastatic prostate cancer
.
N Engl J Med
2019
;
381
:
2506
18
.
7.
Khalaf
DJ
,
Annala
M
,
Taavitsainen
S
,
Finch
DL
,
Oja
C
,
Vergidis
J
, et al
Optimal sequencing of enzalutamide and abiraterone acetate plus prednisone in metastatic castration-resistant prostate cancer: A multicentre, randomised, open-label, phase 2, crossover trial
.
Lancet Oncol
2019
;
20
:
1730
39
.
8.
Beltran
H
,
Hruszkewycz
A
,
Scher
HI
,
Hildesheim
J
,
Isaacs
J
,
Yu
EY
, et al
The role of lineage plasticity in prostate cancer therapy resistance
.
Clin Cancer Res
2019
;
25
:
6916
24
.
9.
Davies
AH
,
Beltran
H
,
Zoubeidi
A
. 
Cellular plasticity and the neuroendocrine phenotype in prostate cancer
.
Nat Rev Urol
2018
;
15
:
271
86
.
10.
Bluemn
EG
,
Coleman
IM
,
Lucas
JM
,
Coleman
RT
,
Hernandez-Lopez
S
,
Tharakan
R
, et al
Androgen receptor pathway-independent prostate cancer is sustained through FGF signaling
.
Cancer Cell
2017
;
32
:
474
89
.
11.
Annala
M
,
Vandekerkhove
G
,
Khalaf
D
,
Taavitsainen
S
,
Beja
K
,
Warner
EW
, et al
Circulating tumor DNA genomics correlate with resistance to abiraterone and enzalutamide in prostate cancer
.
Cancer Discov
2018
;
8
:
444
57
.
12.
Scher
HI
,
Halabi
S
,
Tannock
I
,
Morris
M
,
Sternberg
CN
,
Carducci
MA
, et al
Design and end points of clinical trials for patients with progressive prostate cancer and castrate levels of testosterone: Recommendations of the Prostate Cancer Clinical Trials Working Group
.
J Clin Oncol
2008
;
26
:
1148
59
.
13.
Langmead
B
,
Salzberg
SL
. 
Fast gapped-read alignment with Bowtie 2
.
Nat Methods
2012
;
9
:
357
9
.
14.
Faust
GG
,
Hall
IM
. 
SAMBLASTER: Fast duplicate marking and structural variant read extraction
.
Bioinformatics
2014
;
30
:
2503
5
.
15.
Martin
M
. 
Cutadapt removes adapter sequences from high-throughput sequencing reads
.
EMBnet J
2011
;
17
:
10
.
16.
Quinlan
AR
,
Hall
IM
. 
BEDTools: A flexible suite of utilities for comparing genomic features
.
Bioinformatics
2010
;
26
:
841
2
.
17.
Vandekerkhove
G
,
Struss
WJ
,
Annala
M
,
Kallio
HML
,
Khalaf
D
,
Warner
EW
, et al
Circulating tumor DNA abundance and potential utility in de novo metastatic prostate cancer
.
Eur Urol
2019
;
75
:
667
75
.
18.
Landrum
MJ
,
Lee
JM
,
Benson
M
,
Brown
GR
,
Chao
C
,
Chitipiralla
S
, et al
ClinVar: Improving access to variant interpretations and supporting evidence
.
Nucleic Acids Res
2018
;
46
:
D1062
7
.
19.
Warner
EW
,
Herberts
C
,
Fu
S
,
Yip
SM
,
Wong
A
,
Wang
G
, et al
BRCA2, 
ATM, and CDK12 defects differentially shape prostate tumor driver genomics and clinical aggression
.
Clin Cancer Res
2021
;
27
:
1650
62
.
20.
Vandekerkhove
G
,
Lavoie
J-M
,
Annala
M
,
Murtha
AJ
,
Sundahl
N
,
Walz
S
, et al
Plasma ctDNA is a tumor tissue surrogate and enables clinical-genomic stratification of metastatic bladder cancer
.
Nat Commun
2021
;
12
:
184
.
21.
Van Loo
P
,
Nordgard
SH
,
Lingjærde
OC
,
Russnes
HG
,
Rye
IH
,
Sun
W
, et al
Allele-specific copy number analysis of tumors
.
Proc Natl Acad Sci U S A
2010
;
107
:
16910
5
.
22.
Scher
HI
,
Morris
MJ
,
Stadler
WM
,
Higano
CS
,
Halabi
S
,
Smith
MR
, et al
The Prostate Cancer Working Group 3 (PCWG3) consensus for trials in castration-resistant prostate cancer (CRPC)
.
J Clin Oncol
2015
;
33
:
5000
.
23.
Gundem
G
,
Van Loo
P
,
Kremeyer
B
,
Alexandrov
LB
,
Tubio
JMC
,
Papaemmanuil
E
, et al
The evolutionary history of lethal metastatic prostate cancer
.
Nature
2015
;
520
:
353
7
.
24.
Romanel
A
,
Gasi Tandefelt
D
,
Conteduca
V
,
Jayaram
A
,
Casiraghi
N
,
Wetterskog
D
, et al
Plasma AR and abiraterone-resistant prostate cancer
.
Sci Transl Med
2015
;
7
:
312re10
.
25.
Wyatt
AW
,
Azad
AA
,
Volik
SV
,
Annala
M
,
Beja
K
,
McConeghy
B
, et al
Genomic alterations in cell-free DNA and enzalutamide resistance in castration-resistant prostate cancer
.
JAMA Oncol
2016
;
2
:
1598
606
.
26.
Chen
EJ
,
Sowalsky
AG
,
Gao
S
,
Cai
C
,
Voznesensky
O
,
Schaefer
R
, et al
Abiraterone treatment in castration-resistant prostate cancer selects for progesterone responsive mutant androgen receptors
.
Clin Cancer Res
2015
;
21
:
1273
80
.
27.
Ritch
E
,
Fu
SYF
,
Herberts
C
,
Wang
G
,
Warner
EW
,
Schönlau
E
, et al
Identification of hypermutation and defective mismatch repair in ctDNA from metastatic prostate cancer
.
Clin Cancer Res
2020
;
26
:
1114
25
.
28.
Korpal
M
,
Korn
JM
,
Gao
X
,
Rakiec
DP
,
Ruddy
DA
,
Doshi
S
, et al
An F876L mutation in androgen receptor confers genetic and phenotypic resistance to MDV3100 (enzalutamide)
.
Cancer Discov
2013
;
3
:
1030
43
.
29.
Joseph
JD
,
Lu
N
,
Qian
J
,
Sensintaffar
J
,
Shao
G
,
Brigham
D
, et al
A clinically relevant androgen receptor mutation confers resistance to second-generation antiandrogens enzalutamide and ARN-509
.
Cancer Discov
2013
;
3
:
1020
9
.
30.
Lallous
N
,
Volik
SV
,
Awrey
S
,
Leblanc
E
,
Tse
R
,
Murillo
J
, et al
Functional analysis of androgen receptor mutations that confer anti-androgen resistance identified in circulating cell-free DNA from prostate cancer patients
.
Genome Biol
2016
;
17
:
10
.
31.
Chen
CD
,
Welsbie
DS
,
Tran
C
,
Baek
SH
,
Chen
R
,
Vessella
R
, et al
Molecular determinants of resistance to antiandrogen therapy
.
Nat Med
2004
;
10
:
33
9
.
32.
Viswanathan
SR
,
Ha
G
,
Hoff
AM
,
Wala
JA
,
Carrot-Zhang
J
,
Whelan
CW
, et al
Structural alterations driving castration-resistant prostate cancer revealed by linked-read genome sequencing
.
Cell
2018
;
174
:
433
47
.
33.
Takeda
DY
,
Spisák
S
,
Seo
J-H
,
Bell
C
,
O'Connor
E
,
Korthauer
K
, et al
A somatically acquired enhancer of the androgen receptor is a noncoding driver in advanced prostate cancer
.
Cell
2018
;
174
:
422
32
.
34.
Nyquist
MD
,
Li
Y
,
Hwang
TH
,
Manlove
LS
,
Vessella
RL
,
Silverstein
KAT
, et al
TALEN-engineered AR gene rearrangements reveal endocrine uncoupling of androgen receptor in prostate cancer
.
Proc Natl Acad Sci U S A
2013
;
110
:
17492
7
.
35.
Li
Y
,
Yang
R
,
Henzler
CM
,
Ho
Y
,
Passow
C
,
Auch
B
, et al
Diverse AR gene rearrangements mediate resistance to androgen receptor inhibitors in metastatic prostate cancer
.
Clin Cancer Res
2020
;
26
:
1965
76
.
36.
De Laere
B
,
van Dam
P-J
,
Whitington
T
,
Mayrhofer
M
,
Diaz
EH
,
Van den Eynden
G
, et al
Comprehensive profiling of the androgen receptor in liquid biopsies from castration-resistant prostate cancer reveals novel intra-AR structural variation and splice variant expression patterns
.
Eur Urol
2017
;
72
:
192
200
.
37.
Mateo
J
,
Seed
G
,
Bertan
C
,
Rescigno
P
,
Dolling
D
,
Figueiredo
I
, et al
Genomics of lethal prostate cancer at diagnosis and castration resistance
.
J Clin Invest
2020
;
130
:
1743
51
.
38.
Kumar
A
,
Coleman
I
,
Morrissey
C
,
Zhang
X
,
True
LD
,
Gulati
R
, et al
Substantial interindividual and limited intraindividual genomic diversity among tumors from men with metastatic prostate cancer
.
Nat Med
2016
;
22
:
369
78
.
39.
Lorente
D
,
Mateo
J
,
Zafeiriou
Z
,
Smith
AD
,
Sandhu
S
,
Ferraldeschi
R
, et al
Switching and withdrawing hormonal agents for castration-resistant prostate cancer
.
Nat Rev Urol
2015
;
12
:
37
47
.
40.
Sumiyoshi
T
,
Mizuno
K
,
Yamasaki
T
,
Miyazaki
Y
,
Makino
Y
,
Okasho
K
, et al
Clinical utility of androgen receptor gene aberrations in circulating cell-free DNA as a biomarker for treatment of castration-resistant prostate cancer
.
Sci Rep
2019
;
9
:
4030
.
41.
Conteduca
V
,
Wetterskog
D
,
Sharabiani
MTA
,
Grande
E
,
Fernandez-Perez
MP
,
Jayaram
A
, et al
Androgen receptor gene status in plasma DNA associates with worse outcome on enzalutamide or abiraterone for castration-resistant prostate cancer: A multi-institution correlative biomarker study
.
Ann Oncol
2017
;
28
:
1508
16
.
42.
Jayaram
A
,
Wingate
A
,
Wetterskog
D
,
Conteduca
V
,
Khalaf
D
,
Sharabiani
MTA
, et al
Plasma androgen receptor copy number status at emergence of metastatic castration-resistant prostate cancer: A pooled multicohort analysis
.
JCO Precis Oncol
2019
;
3
:
PO.19.00123
.
43.
Conteduca
V
,
Jayaram
A
,
Romero-Laorden
N
,
Wetterskog
D
,
Salvi
S
,
Gurioli
G
, et al
Plasma androgen receptor and docetaxel for metastatic castration-resistant prostate cancer
.
Eur Urol
2019
;
75
:
368
73
.
44.
Cescon
DW
,
Bratman
SV
,
Chan
SM
,
Siu
LL
. 
Circulating tumor DNA and liquid biopsy in oncology
[Internet].
Nature Cancer
2020
;
1
:
276
90
.
45.
Choudhury
AD
,
Werner
L
,
Francini
E
,
Wei
XX
,
Ha
G
,
Freeman
SS
, et al
Tumor fraction in cell-free DNA as a biomarker in prostate cancer
.
JCI Insight
2018
;
3
:
e122109
.
46.
Mazrooei
P
,
Kron
KJ
,
Zhu
Y
,
Zhou
S
,
Grillo
G
,
Mehdi
T
, et al
Cistrome partitioning reveals convergence of somatic mutations and risk variants on master transcription regulators in primary prostate tumors
.
Cancer Cell
2019
;
36
:
674
89
.e6.

Supplementary data