Gene Copy Number Estimation from Targeted Next-Generation Sequencing of Prostate Cancer Biopsies: Analytic Validation and Clinical Qualification

Prostate cancer is a highly heterogeneous disease with distinct genomic underpinnings (1) and generally presents a modest number of single point mutations and small insertions and deletions (indels) at approximately four per megabase (2, 3), but frequently has large-scale copy number and structural alterations (4, 5). We recently showed that 20% to 30% of prostate tumors bear defects in DNA repair genes that render them sensitive to treatment with PARP inhibition, and that these are commonly mediated through deletions in these tumor suppressors (6).

High-throughput, inexpensive, multiplex biomarker assays for molecular stratification are needed to guide therapeutic choices and support clinical trial design. Although high-coverage whole-exome or genome sequencing can reliably assess tumor genomics, cost and bioinformatic demands currently restrict their routine clinical implementation to stratify patients for treatment. Targeted next-generation sequencing (NGS) represents an opportunity for implementing genomics in clinical practice with the advantages of a rapid turnaround time, lower cost, and the ability to concurrently analyze multiple genes with a limited bioinformatic analysis burden. Existing approaches to call germline and somatic single-nucleotide variants (SNV) and short indels from targeted sequencing data are well established and have led to a number of FDA biomarker-dependent drug approvals, including the PARP inhibitor olaparib in the past decade (7, 8). However, in the case of PARP inhibition for homologous recombination repair defective cancers, screening for tumors with these DNA repair defects relies on detecting gene copy number changes.

In this study, we identify gene copy number aberrations (CNA), based on targeted amplicon sequencing of a focused biomarker gene panel dedicated to identifying DNA repair defects. We used a patient cohort of metastatic castration-resistant prostate cancer (CRPC) to explore the utility of this approach. This targeted gene panel focused on the study of 113 genes important to DNA repair machinery, including BRCA2, BRCA1, ATM, PALB2, and CDK12 as well as other important potentially actionable prostate cancer genes, including AR, SPOP, PIK3CA, PTEN, AKT1, AKT2, and MYC (1, 6).

Patient tumor (biopsies) and germline (buccal swab and saliva) samples were collected as part of the Royal Marsden ethics committee–approved CCR2472 protocol. Samples were collected, annotated, stored, and reviewed as described previously (6). All samples were collected between July 25, 2015, and October 17, 2016.

Briefly, biopsies were paraffin embedded, DNA from biopsies was extracted using the QIAamp DNA FFPE Tissue Kit (Qiagen), Quant-iT Picogreen High-Sensitivity double-stranded DNA (dsDNA) Assay Kits (Invitrogen, Thermo Fisher Scientific) and QC check by FFPE QC Kit (Illumina).

A customized Generead v2 DNAseq Panel (Qiagen) panel was used for library construction; the exonic regions of 113 genes were included in the targeted panel, these being selected for being potentially actionable and/or involved in DNA damage repair processes and/or in prostate carcinogenesis (Supplementary Table S1; ref. 6). The panel covered 564 kb of the human genome, and each of the 113 genes was covered by an average of 18.69 probes (SD, 13.64).

Targeted amplicon sequencing was performed using the Illumina MiSeq platform following the manufacturers' protocol. FASTQ files were generated using the Illumina MiSeq Reporter v2.5.1.3 (9, 10). Reads were aligned to the human reference genome (hg19) using BWA. Whole-exome sequencing (WES) was performed at the University of Michigan (Ann Arbor, MI) as described previously (1).

FastQC (v0.11.2) and Samtools (0.1.19) were used to assess sequencing quality. Targeted sequencing samples were excluded from copy number analysis on the basis of: insufficient (<0.5 million) total read counts, low (<95%) percentage of properly paired reads, and low (<99.99%) percentage of on-target reads. Exome sequencing samples were rejected if the FastQC per-base quality score for 75% of the reads was less than Q20 over the first 80 bases, and alignment quality was monitored with Picard, as described previously (1).

We used the software CNVkit (v0.3.5; ref. 11) to analyze sequencing coverage and copy number in the aligned sequencing reads from targeted amplicon sequencing of tumor and germline samples. Sequencing coverage of targeted regions in germline samples was assessed and used to create pooled reference data that included the technical variability at each covered region. Regions that were poorly amplified or mapped were masked from further analysis. The analyses of germline samples were also performed with CNVkit to validate sample quality. The read depths of tumor samples were accessed, normalized (corrected for GC content, target footprint size and spacing, and repetitive sequences), and individually compared with the reference, and the circular binary segmentation (CBS) algorithm (12) was used to infer copy number segments. Very small copy changes (≤3 bins) were treated as artifacts. Poor samples were also identified through high (>0.8) values of the segment interquartile range (IQR), a metric for the spread of bin-level copy ratios within each segment following application of the CBS algorithm. Copy number segments were annotated to genes, and regions bearing a log₂ ratio of at least ±0.4 were identified as suggestive of shallow deletions or gains (11). Segments with log₂ <−1.2 were classified as deep deletions, and those with log₂ >2 were classified as amplifications. Experimental noise was identified as a log₂ ratio SD of approximately 0.2.

WES CNA data were extracted from a previously reported Su2C-PCF dataset (1). Briefly, log₂ ratios were derived on a per-gene basis from circular-binary segmented, Lowess-normalized log₂ transformed coverage ratio between each tumor and matched normal sample.

Tumor DNA from prostate cancer patients and male reference DNA from Agilent were amplified using Sigma WGA2 Kit (Sigma-Aldrich) according to the manufacturer's recommendations, and the amplified DNA was quantified by the Qubit fluorometric quantitation method (Thermo Fisher Scientific). Amplified tumor DNA (500 ng) was then fluorescently labeled with Cy5, and male reference DNA labeled with Cy3, using the SureTag Complete DNA Labeling Kit (Agilent Technologies). Labeled DNA was subsequently hybridized utilizing the Agilent SurePrint G3 Human CGH Microarray Kit, 4 × 180 K according to the manufacturer's instructions. The slides were then scanned and analyzed using the CytoGenomics Software v 4.0.3.12 (Agilent Technologies). To compare array comparative genomic hybridization (aCGH) and NGS results, we used the log₂ ratio of the aCGH segments that overlap with genes present in the NGS panel. Genes not part of a copy-altered segment were included in the analysis by using the mean log₂ ratio for probes, if the value was no greater than 0.25. Genes on the X-chromosome were excluded from comparison with results from the male reference–based panel data as a female reference was used for aCGH.

RB1 FISH was performed using a standard formalin-fixed, paraffin-embedded (FFPE) hybridization method (13) on 3-μm FFPE tissue slices adjacent to hematoxylin and eosin sections that were confirmed to contain a minimum of 50 intact cells. Briefly, RB1 status was determined using Vysis LSI 13 RB1 (13q14) probe (catalog # 08L65-020; Abbott Laboratories) and a reference probe Vysis 13q34 (catalog # 05N34-020; Abbott Laboratories). Nineteen z-stacks were used to assess cell status.

Freshly collected metastatic CRPC (mCRPC) biopsies (n = 110) were collected from patients with progressing disease (Supplementary Table S2) and assessed for copy number changes by sequencing with our targeted panel. Germline DNA samples (buccal swab) were taken from 34 consenting patients to use as a baseline reference and for additional analyses. Tumor biopsies were from: bone (49.1%), lymph nodes (35.5%), liver (9.1%), soft tissue disease (3.6%), and transurethral resection of the prostate (3.6%).

Accurately assessing the depth of coverage of genes can be biased by high GC content and repetitive regions (11, 14), so we first evaluated CNA reproducibility on a per-gene basis between technical replicates. Thirteen samples were sequenced in duplicate to assess technical variation; 12 of these duplicate datasets included repeated library preparation from the same DNA extraction, while one involved resequencing of the same library twice.

We identified 15 genes that produced less concordant results between replicates, including five genes with insufficient number of targeted regions (four or less probes) per gene (NRAS, NFKBIA, FANCF, FAM46C, CDKN1B), and 10 (MAP2K2, CDKN2A, HRAS, RECQL4, NOTCH1, FGFR3, STK11, TSC2, SMARCA4, AXIN1) that contained repetitive or polymorphic exonic regions that impeded accurate read alignment and possessed high coverage variability (Supplementary Fig. S1). These genes were excluded from downstream CNA analyses.

Following removal of the aforementioned genes, CNA estimates between duplicates showed high correlation (Pearson r correlation coefficient = 0.92; Fig. 1A). All genes with deep deletions (log₂ ratio, <−1.2) or amplifications (log₂ ratio > 2) had a similar (at least ±0.4) CNA estimation replicated in the duplicate sample. A total of 95% of samples with a log₂ ratio change of at least ±0.4 had a similar result in the duplicate.

A matched germline sample is commonly used for identifying deviations in coverage ratios, so we sought to compare this with our pooled reference approach. For the pooled reference, we sequenced 34 unmatched germline samples and assessed the depth of coverage as described previously (11). Regions with low variability were identified in the pooled reference and used to weight downstream CNA estimates (11). This approach produced a highly similar result (overall Pearson r = 0.93) when compared with the CNA estimation approach using a matched normal as reference, when matched normal samples were available (Fig. 1B); however, the CNA estimation confidence score (IQR value) was improved when using pooled germline samples (paired t test P = 0.05). This was in line with other reports (14) that indicated that using pooled germline sequencing data as reference offered a robust method to identify gene CNAs.

We cross-compared copy number log₂ ratios for all evaluable genes on the targeted panel with log₂ ratios generated by exome sequencing (n = 13; Fig. 1C), resulting in a Pearson r correlation coefficient of 0.86. Moreover, all the genes with a log₂ ratio <−1.2 by targeted sequencing, in keeping with putative “deep” deletions, exhibited a similar (at least ±0.4) result in the exome sequencing. For example, one sample bore AR amplification (log₂ ratio 5.7 in targeted NGS and 4.91 in WES) alongside WRN deletion (−1.53 vs. −1.17) and ATM deletion (−0.68 vs. −0.56). We also compared log₂ ratios from aCGH segments (n = 9) with gene copy number estimations by panel NGS among the genes covered by both platforms and found log₂ ratios were also concordant (Fig. 1D), with a Pearson r value of 0.7.

We then sought to further validate copy number calls for specific genes of clinical relevance, such as BRCA2, MYC, PIK3CA, PTEN, and RB1 (Fig. 1E and F). For these genes, we found that our panel produced results highly concordant with exome sequencing and aCGH data, with Pearson r values of 0.91 and 0.84, respectively. All clinically relevant deletions of BRCA2 in samples with exome or aCGH data were also identified in our panel results (Table 1).

In the analyses of this cohort of advanced prostate cancer samples, the prevalence of the CNAs detected with this assay was consistent with that of previously reported studies (Fig. 2). CNAs of AR (amplification; log₂ > 2) were found in 45.5% (50/110) of tumor samples. Broad deletions of chromosome 13 (1) involving the shared loss of BRCA2 and RB1 were detected as previously described with 23 of 26 (88%) samples with any BRCA2 loss also having some loss of RB1, while overall 23 of 66 (34.8%) samples with RB1 loss had loss of BRCA2 (Fig. 2A).

Overall, the evaluation of the deep deletions of commonly aberrant CRPC genes revealed results in keeping with previous reports from two recent studies of CRPC genomics (1, 15); PTEN genomic deep deletions were detected in 17 of 110 (15.5%), BRCA2 deep deletion in 5 of 110 (4.5%), and RB1 deep deletion in 9 of 110 (8.2%) CRPC samples. In addition, gain of the MYC locus was found in 36 of 110 samples (32.7%). We found that the copy number alteration frequencies in our cohort were in line with expected CRPC frequencies.

To further explore how tumor purity altered the capacity of the assay to detect CNAs, we pursued the serial dilution of tumor DNA acquired from a sample estimated by pathology review as having 80% tumor content, with same-patient germline DNA. This sample harbored a somatic, BRCA2 homozygous deletion (independently confirmed WES). We then sequenced the resulting dilutions (Fig. 3A) and found that homozygous deletions were easily distinguishable with tumor content purities >60% and that clonal CNAs were detectable with purity as low as 30%.

We also performed microdissection of a tumor biopsy (pathologist tumor purity estimate, ∼40%) bearing borderline BRCA2 homozygous deletion (log₂ ratio of −1.15) and resequenced the resulting tumor-enriched sample (purity estimate, ∼80%). This step further shifted the CNA estimates away from normal, with the log₂ ratio for BRCA2 dropping to −2.82, an unambiguous deletion (Fig. 3B).

The log₂ copy ratios estimated from sequencing a bulk of cells, as in our protocol, can be biased by tumor sample heterogeneity. To assess this relationship, we performed FISH assays of the RB1 gene, counting the number of cells where each copy number (0 through 3+) was observed (Fig. 3C). Our results showed an association between the proportion of cells with 1 or 0 copies of RB1 (hemizygous or homozygous deletion) and putative copy loss inferred from our targeted sequencing panel (unpaired t test P = 0.02, n = 18, df = 15.884). However, the full complexity of aneuploidy in individual tumor cells is masked by bulk tissue sequencing, which effectively averages the estimated DNA content across all cells in a sample. FISH analyses revealed that many tumor samples contained a mixture of cells with different copy number changes at the RB1 locus, in highly heterogeneous cell populations (Supplementary Fig. S2).

The ability to rapidly, inexpensively, and accurately identify cancers with actionable genomic defects in tumor biopsies is now becoming critically important for advanced prostate cancer care. We have previously shown that a proportion of metastatic CRPCs have DNA repair defects and are sensitive to treatment with PARP inhibitors (6), with emerging studies of platinum-based chemotherapies also indicating important antitumor activity in this population (16). These and other data support further studies of genomic stratification for treatment selection of patients suffering from mCRPC. Here, we demonstrate the analytic validation and clinical use of a targeted NGS assay that allows for copy number screening for predictive biomarkers of response, enabling the integration of more precise patient stratification into trial design.

Our method reliably estimates gene copy number changes from amplicon-based targeted NGS with the acquired results being comparable with widely used assays, including exome sequencing, aCGH, digital droplet PCR, and FISH (17). Critically, the data processing time for this panel-based CNA estimation is substantially shorter of that for exome sequencing and can be easily completed even with a personal computer. In addition, we have confirmed that using pooled germline samples as a reference can produce equivalent and robust results to utilizing matched normal samples. This feature provides the possibility of diagnosing patients with only tumor samples for whom germline DNA is not available, and can dramatically reduce the overall sequencing costs of large studies.

This approach is currently being applied in the TOPARP-B (NCT01682772) study as part of an ongoing attempt to clinically qualify DNA repair pathway aberrations as a predictive biomarker for treatment with PARP inhibition with olaparib (16). However, this assay has utility beyond this study. The gene panel we utilized covers several key genomic aberrations that are of emerging interest as possible clinical predictive biomarkers that may be actionable in prostate cancer, including loss of PTEN, gain of PIK3CA (both may sensitize tumors to PI3K/AKT pathway inhibition; refs. 18, 19), and amplification of MYC (MYC-driven tumors have been reported to sensitize to AURA and BET inhibition and can result in replication stress; ref. 20). Prospectively identifying tumors bearing these CNAs and other deleterious mutations, utilizing such validated methods, is likely to become critically important in defining future treatment strategies for this commonest of male cancers.

We explored the functional limits of assessing CNAs and found that our pipeline can confidently estimate deletion when tumor content is reduced to 60% (homozygous and heterozygous). As tumor content decreases, the pipeline struggles to differentiate between homozygous and heterozygous deletions, and we estimate that interpreting results from samples with <30% purity may be extremely challenging for targeted NGS panels. We have shown, however, that the use of microdissection to enrich tumor content offers a possible solution to resolve these issues.

This approach also suffers from the nature of amplicon-based NGS, including difficulty identifying precise breakpoints, very focal gains or deletions, unequal amplification among genes, and coverage variation across polymorphic repeat regions. Furthermore, because the number of SNVs identified in targeted panels will be much lower when compared with exome or whole-genome sequencing, SNV allele frequencies may not be sufficient to assist with copy number estimation or to detect allelic imbalance.

It is important to note that we did not systematically compare software tools in this study and that different clinical requirements will require individual optimizations. As the field develops, other software tools will emerge to improve the precision of CNA identification in targeted sequencing, and independent validations will be critical in guiding further work.

In conclusion, we have described the validation and clinical application of an assay evaluating 100 genes utilized to assess both deleterious tumor gene copy number, and mutations, from targeted NGS of prostate cancer samples from the most common metastatic sites, including bone biopsies. We have validated these results by orthogonal methods, including not only exome sequencing, but also aCGH and FISH. Such targeted panel NGS assays are likely to become central to future patient care in cancer therapy.

S. Sandhu is a consultant/advisory board member for Amgen, Bristol-Myers Squibb, and Janssen. J.S. de Bono is a consultant/advisory board member for AstraZeneca. No potential conflicts of interest were disclosed by the other authors.

Conception and design: G. Seed, W. Yuan, J. Mateo, S. Carreira, G. Boysen, S. Miranda, S. Sandhu, A.M. Chinnaiyan, J.S. de Bono

Development of methodology: G. Seed, S. Carreira, S. Miranda, E. Talevich, S. Sandhu, A.M. Chinnaiyan, J.S. de Bono

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): G. Seed, J. Mateo, S. Carreira, C. Bertan, R. Ferraldeschi, M. Crespo, D.N. Rodrigues, D.R. Robinson, Y.-M. Wu, S. Sandhu, J.S. de Bono

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): G. Seed, W. Yuan, J. Mateo, S. Carreira, M. Lambros, D.N. Rodrigues, E. Talevich, D.R. Robinson, L.P. Kunju, Y.-M. Wu, R. Lonigro, J.S. de Bono

Writing, review, and/or revision of the manuscript: G. Seed, W. Yuan, J. Mateo, S. Carreira, M. Lambros, G. Boysen, S. Miranda, I. Figueiredo, R. Riisnaes, D.N. Rodrigues, E. Talevich, R. Lonigro, S. Sandhu, A.M. Chinnaiyan, J.S. de Bono

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): G. Seed, I. Figueiredo, D.N. Rodrigues, J.S. de Bono

Study supervision: W. Yuan, J.S. de Bono

We thank the staff at the Prostate Targeted Therapy Group and all the patients who agreed to participate in this study.

This study was supported by Cancer Research UK, Prostate Cancer UK, and the Stand Up To Cancer – Prostate Cancer Foundation Prostate Dream Team Translational Research Grant (grant number: SU2C-AACR-DT0712). Research grant administered by the American Association for Cancer Research, the Scientific Partner of SU2C. Support was also provided by an Experimental Cancer Medicine Centre grant and a Medical Research Council/Prostate Cancer UK and PCF Young Investigator Award to J. Mateo. G. Seed was supported by a Prostate Cancer UK PhD studentship (TLD-S15-006), and G. Boysen was supported by a Marie Curie International Incoming Fellowship (625792). This study was also supported by research funding from the EU FP7 project CTCTrap #305341.

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.