Development and validation of a 36-gene sequencing assay for hereditary cancer risk assessment

Institutional review board approval

The protocol for this study was approved by Western Institutional Review Board (IRB number 1145639) and complied with the Health Insurance Portability and Accountability Act (HIPAA). The information associated with patient samples was de-identified in accordance with the HIPAA Privacy Rule. A waiver of informed consent was requested and approved by the IRB.

Multigene panel design

Thirty-six genes associated with hereditary forms of cancer, including breast, ovarian, colorectal, gastric, endometrial, pancreatic, thyroid, prostate, melanoma, and neuroendocrine, were selected for development of the Counsyl Inherited Cancer Screen panel. The genes are: APC, ATM, BARD1, BMPR1A, BRCA1, BRCA2, BRIP1, CDH1, CDK4, CDKN2A, CHEK2, EPCAM, GREM1, MEN1, MLH1, MRE11A, MSH2, MSH6, MUTYH, NBN, PALB2, PMS2, POLD1, POLE, PTEN, RAD50, RAD51C, RAD51D, RET, SDHA, SDHB, SDHC, SMAD4, STK11, TP53, and VHL (Table 1). Twenty nine of the 36 genes were specifically included due to the availability of patient management guidelines by NCCN or other professional societies. Further details regarding the panel are available in Table S1.

The selected genes are tested for SNVs, indels, and CNVs throughout coding exons and 20 bp of flanking intronic sequences. Additionally, known deleterious variants outside the coding regions are sequenced. In EPCAM, only large deletions that include exon 9 are reported as these mutations are known to silence the MSH2 gene (Tutlewska, Lubinski & Kurzawski, 2013). In GREM1, specific pathogenic duplications in the promoter, which are commonly associated with individuals of Ashkenazi Jewish descent, are covered. Specifically, the screen targets the three most common promoter duplications in GREM1 (coordinates with respect to GRCh37/hg19 reference assembly):

• chr15:32,964,939-33,004,759 (40kb)

• chr15:32,986,220-33,002,449 (16kb)

• chr15:32,975,886-33,033,276 (57kb).

For PMS2, exons 11–15 are excluded from the reportable region of interest (ROI) because of high similarity between this portion of PMS2 and its highly homologous pseudogene PMS2CL. In RET, exon 1 is not sequenced due to high guanine-cytosine (GC) content.

Next Generation DNA Sequencing

Our application of next-generation DNA sequencing is performed as described previously (Kang et al., 2016). Briefly, DNA from a patient’s blood or saliva sample is isolated, quantified by a dye-based fluorescence assay, and then fragmented to 200–1,000 bp by sonication. The fragmented DNA is converted to a sequencing library by end repair, A-tailing, and adapter ligation. Samples are then amplified by PCR with barcoded primers, multiplexed, and subjected to hybrid capture-based enrichment with 40-mer oligonucleotides (Integrated DNA Technologies, Coral, IL, USA) complementary to targeted regions. Next generation sequencing of the selected targets is performed with sequencing-by-synthesis on the Illumina HiSeq 2500 instrument to a mean sequencing depth of ∼650×. All target nucleotides are required to be covered with a minimum depth of 20 reads.

Bioinformatics processing

Sequencing reads are aligned to the hg19 human reference genome using the BWA-MEM algorithm (Li, 2013). Single-nucleotide variants and short indels are identified and genotyped using GATK 1.6 and FreeBayes (McKenna et al., 2010; Garrison & Marth, 2012). The calling algorithm for copy number variants is described below. All SNVs, indels, and large deletions/duplications within the reportable range are analyzed and classified by the method described in the section “Variant Classification”. All reportable calls are reviewed by licensed clinical laboratory personnel.

CNV calling algorithm

Copy number variants for samples are determined by inspecting the number of mapped reads observed at targeted positions in the genome across samples in a flowcell lane. Our method is based upon previous successful approaches applying hidden Markov models (HMMs) to exome sequencing data (Plagnol et al., 2012) with modifications presented below that have been optimized for accurate resolution of CNVs based on the particulars of the sequencing technology. As sequencing depth is linearly proportional to the number of copies of the genome at that position, we construct a statistical model for the likelihood of observing a given number of mapped reads d_i,j at a given genomic position i for sample j with copy number c_i,j.

The expected number of reads is dependent upon 3 factors: the average depth for that targeted location across samples μ_i, the average depth for that particular sample across targeted positions μ_j, and the local copy number of the sample’s genome at that targeted position. These are first determined by finding the median depth at targeted region across all N_s samples in an analyzed flowcell lane ${\displaystyle {\mu }_i=\frac{\sum _j{d}_{i,j}}{N_s}}$ then the sample dependent factor μ_j is found by taking the median across all N_p positions in genome after normalizing for the expected number of reads at each position ${\displaystyle {\mu }_j=\frac{\sum _i{d}_{i,j}∕{\mu }_i}{N_p}.}$ Combining these factors the observed data are modeled by the negative binomial distribution ${\displaystyle p\left(d_{i,j}|c_{i,j}\right)=NegBinom\left(d_{i,j}|\mu =c_{i,j}{\mu }_i{\mu }_j,r=r_i\right).}$ This characterization has been found to accurately model the observed number of reads from previous targeted sequencing experiments (Anders & Huber, 2010).

In the negative binomial model, the variance parameter r_i accounts for regions of the genome where sequencing depth is observed to follow idealized Poisson statistics in the limit that r → ∞ and regions that are excessively noisy with respect to observed number of reads when r → 0.

r_i may be estimated as ${\displaystyle r_i=\frac{{\mu }_i^2}{Va{r}_j\left[d_{i,j}\right]-{\mu }_i}}$ which is found to closely model the empirical distribution over several orders of magnitude in read depth.

Because duplications and deletions will simultaneously impact the expected depth of all genomic positions encompassing the variant, depth data from spatially adjacent positions are correlated. We leverage the HMM to account for this correlation. The HMM’s state transition probabilities between wild-type and copy-number-variant are parameterized by matching the average length of such variations observed in human population (Sudmant et al., 2015) through setting p_CNV→WT = 1∕6,200 between each subsequent base-pair and a prior on the frequency of such variations ${\displaystyle \frac{p_{WT\to CNV}}{p_{CNV\to WT}}=p_{CNV}.}$ The prior p_CNV = 0.001 was determined by balancing the thresholds for confident calling and retesting of calls to achieve the desired sensitivity and specificity, and the prior was set independently of this validation.

Detecting CNVs using this probabilistic framework invokes the Viterbi algorithm (Korn et al., 2008) to determine the most likely number of copies at every targeted region within a sample. Any contiguous regions of duplication or deletion produce a reported variant, and the confidence of that call is determined by aggregating the posterior probability of the call ∑_iϵCNVp(c_i,j ≠ 2) not being wildtype over the called region.

All called copy number variants are inspected for quality by human review. To avoid false positives, all patient samples with a called CNV are tested a second time, starting with a new DNA extraction and including library preparation, sequencing and bioinformatic analysis. Samples that emit low confidence called variants are additionally rerun to resolve a confident genotype.

Detection of Alu insertions

Alu positives were detected by looking for Alu sequences in reads overlapping with Alu insertion positions. All insertions were only tested at positions where the sequence had been previously confirmed by Sanger sequencing. At the site of an Alu insertion, the Alu sequence is soft-clipped by BWA alignment. These soft-clipped reads were compiled; duplicate reads were discarded; and the remaining reads with sequences matching the known Alu sequence at this site were tallied. Sites with at least three unique reads matching the Alu sequence were called as Alu positive.

Pre- and post-sequencing quality metrics

To ensure the quality of the results obtained from the assay, 27 different review checkpoints (Table S2) were developed. Ancillary quality-control metrics are computed on the sequencing output and used to exclude and re-run failed samples, and include the fraction of sample contamination (<5%), extent of GC bias, read quality (percent Q30 bases per Illumina specifications), depth of coverage (per base minimum coverage ≥20 × and mean coverage of >250×), and region of interest (ROI) coverage (100%). Calls that do not meet criteria listed in Table S2 are set to “no-call.” To ensure clinical calling accuracy, all calls and no-calls for potentially deleterious variants, and all calls for variants of unknown significance are manually reviewed by laboratory personnel and are subject to override if warranted, based on a pre-established protocol.

Variant classification

Variants are classified using multiple lines of evidence according to the ACMG Standards and Guidelines for the Interpretation of Sequence Variants (American College of Medical Genetics and Genomics, 2015; Richards et al., 2015). Variants that are known or likely to be pathogenic are reported; patients and providers have an option to have variants of uncertain significance reported as well. Final variant classifications are regularly uploaded to ClinVar (Landrum et al., 2014), a peer-reviewed database created with a goal of improving variant interpretation consistency between laboratories.

Statistical analysis

Variant calls were defined as true positive for variants identified by the Counsyl Inherited Cancer Screen and by independent testing (the 1000 Genomes Project or MLPA/Sanger data), false positive for variants identified by the Counsyl test but not by the independent data, and false negative for variants identified by the independent data but not by the Counsyl test. To estimate true negatives, we counted polymorphic sites (positions at which we observed non-reference bases in any sample) with concordant negative results across all considered samples. No-calls were censored from the analysis. As no-calls have the potential to introduce clinically relevant false negatives, we separately examined the no-calls containing potentially deleterious alleles by treating no-calls as homozygous reference and comparing to the 1000 Genomes calls. We found all no-calls when treated as homozygous reference were concordant with the exception that one comparison was inconclusive due to low allele balance in both our data and the exome data from the 1000 Genomes Project (Table S8).

Validation metrics were defined as: Accuracy = (TP +TN) / (TP +FP +TN +FN); Sensitivity = TP / (TP +FN); Specificity = TN / (TN +FP); FDR = FP / (TP +FP), where TP = true positives, TN = true negatives, FP = false positives, FN =false negatives, and FDR = false discovery rate. The confidence intervals (CIs) were calculated by the method of Clopper and Pearson (Clopper & Pearson, 1934). To estimate reproducibility within and between runs, the ratio of concordant calls to total calls was calculated.

Study samples

The validation sample set comprised (a) 111 genomic DNA reference materials purchased from the Coriell Cell Repositories (Camden, NJ) (Table S3), (b) MLH1/MSH2 exon copy number reference panel from the National Institute for Biological Standards and Control (N = 7) (Table S4), and (c) 223 deidentified patient samples used for MLPA- and Sanger-based confirmation (Tables 2 and S4).

The validation set included samples with reference data for SNVs and indels (the 1000 Genomes Project), a broad range of indels (both short ≤10 bp and long >10 bp) characterized by Sanger sequencing, homopolymer-associated variants, Alu element insertions, and both single- and multi-exon copy-number variants characterized by MLPA (Table 3). Validation material was derived from cell lines, blood, and saliva samples. Collectively, the validation set provides broad coverage of known relevant types of genomic variation across the reportable region of the test (Tables 3 and S5). A list of the validation samples from Sanger and MLPA confirmation is provided in Table S4.

Test description

We developed an NGS-based test that interrogates 36 genes associated with hereditary cancer risk (Table 1). The majority of the 36 genes were selected based on the availability of patient management guidelines developed by NCCN or other professional societies. The reportable region of interest (ROI) of the test is 124,245 bp representing coding exons, intron boundaries and non-exonic mutation-containing regions (Table 1). The wet lab protocols and reagents are carefully optimized to ensure 100% coverage of targeted base pairs at an average depth of 650 reads and a minimal depth of 20 reads sufficient for robust detection of multiple classes of genomic alterations: single-nucleotide variations, indels, and copy number variations.

Validation approach

Several regulations, including the Clinical Laboratory Improvement Act of 1988 (CLIA), the ACMG guidelines for analytical validation of NGS methods (Rehm et al., 2013), as well as various quality standards for diagnostic laboratories require rigorous analytical validation of panel tests for clinical use. In contrast to diagnostic assays for a single gene or a limited panel of genes (syndrome-based testing), analytical validation of a NGS-based test assaying 36 genes for multiple types of genomic alterations is a complex task. To address this challenge, we developed a representative validation approach with reference samples selected to cover variant and specimen variability that may affect test accuracy and reproducibility for clinical use.

To measure the accuracy of SNV and indel detection, we tested samples from the 1000 Genomes Projects with reference data for SNVs and indels in all 36 genes. Testing on the 1000 Genomes Project samples allows us to assess the ability to call commonly observed variant types and the ability to test calling in regions that may be difficult for NGS due to considerable sequence homology (e.g., CHEK2, SDHA, and PMS2) or low complexity (homopolymer runs). However, the 1000 Genomes reference samples provide limited validation for technically challenging variants like CNVs, larger indels, and Alu insertions. To build a collection of reference material to test such challenging variants, we identified relevant patient samples tested with a previous version of the Counsyl test (a 24-gene panel) and orthogonally confirmed each of the positive samples by either Sanger or MLPA. Using these cohorts of reference samples (e.g., samples with CNVs), we could then assess call accuracy for each type of technically challenging variant on this newly designed 36-gene panel. Finally, to validate test reproducibility, we examined SNV, indel, and CNV calls in cell line and patient (blood and saliva) samples processed independently in several batches (inter-run reproducibility) or tested repeatedly in the same batch (intra-run reproducibility).

Analytical validation for SNVs and indels

The analytical validation of the Inherited Cancer Screen was performed according to ACMG guidelines (Rehm et al., 2013) and in accordance with the requirements of CLIA for medical laboratories. SNV and indel detection was examined on a 101-sample validation set consisting of reference samples from the 1000 Genomes Project with known SNV and indel sites across the targeted regions (Tables 2 and S5). Counsyl sequence data for 36 genes were compared to reference data obtained from the 1000 Genomes Projects. Out of 42,925 total calls validated, 18 calls were discordant between Counsyl and the 1000 Genomes Project (Table S6). One of the 18 discordances was a potential false positive variant call, identified as a variant by the Counsyl test, but identified as reference by the 1000 Genomes Project. The remaining 17 calls were potential false negative variants identified by the 1000 Genomes Project, but not by the Counsyl test. Manual review of the 1000 Genomes reference data for each of the discordant sites using the Integrated Genomics Viewer (IGV) (Robinson et al., 2011; Thorvaldsdóttir, Robinson & Mesirov, 2013) found that a large portion of the discordant calls came from hard-to-sequence (e.g., highly homologous SDHA gene) or low-coverage regions, which is a reported limitation in the 1000 Genomes Project (1000 Genomes Project Consortium et al., 2012). With that in mind, each of the discordant sites was subjected to Sanger sequencing as an independent testing method and the data from Sanger sequencing supported all 18 of Counsyl’s calls as true positives or true negatives (Table S6).

Analytical validation results of Counsyl’s test for SNV and indel detection is presented in Table 4. Counsyl’s test identified 5,182 true positive calls, 37,743 true negative calls, and no false positive nor false negative calls, resulting in 100% sensitivity (95% CI [100%–99.93%]), 100% specificity (95% CI [100%–99.99%]), and 0% FDR (95% CI [0–0.0007%]) of the test for detecting SNVs and indels.

Validation of challenging variants

Challenging indels

To measure accuracy for detecting indels, we built a cohort (N = 82) of patient samples with variants of a range of sizes, including both short (≤10 bp) and the more technically challenging long (>10 bp) deletions or insertions (Tables 3 and S4a). These samples were identified using a previous version of the Counsyl test (a 24-gene panel) and orthogonally confirmed by Sanger. We then tested these samples with the newly developed 36-gene panel and confirmed all of the expected indel calls; no false-positives nor false-negatives were observed in the 36-gene panel results (Table 5).

Alu insertions

Alu elements represent a special class of insertions and are known to be clinically important (Belancio, Roy-Engel & Deininger, 2010). Alu insertions have been reported in ATM, BRCA1, BRCA2, and BRIP1 (Belancio, Roy-Engel & Deininger, 2010; Kennemer et al., 2016), including known examples of Alu insertion founder mutations (e.g., c.156_157insAlu in BRCA2 exon 3 in Portuguese populations) (Peixoto et al., 2014). Accurate detection of Alu insertions is challenging, especially for traditional Sanger sequencing where longer Alu-containing alleles are usually out-competed during PCR (De Brakeleer et al., 2013). To test the sensitivity of our assay and bioinformatics pipeline for Alu insertion detection, we included 7 positive cases (Portuguese founder mutation in exon 3 of BRCA2, Alu insertion in BRCA2 exon 25 and intronic Alu insertions in ATM and MSH6) in our validation study (Table 6). We confirmed that the Alu insertions identified by the Counsyl Inherited Cancer Screen were also detected by Sanger sequencing.

Reproducibility

In addition to establishing the test’s analytical sensitivity and specificity, Counsyl’s Inherited Cancer Screen was validated for intra- and inter-run call reproducibility. Intra-run reproducibility of SNV and indel calls was established by testing eight cell lines and 13 blood or saliva samples in 2–3 replicates in the same batch, split across sequencer lanes. Inter-run reproducibility was validated by testing eight cell lines and 84 patient blood or saliva samples in 2–3 different batches run by two operators, on different instruments and on different dates (Table S7a). Concordance between replicates was >99.99%, with just one discordant call at a known benign homopolymer site in an intron of ATM (Table S7a).

For CNVs, intra-run and inter-run reproducibility was established using the Coriell sample NA14626 with a duplication of BRCA1 exon 12 (Table S7b). Concordance between eight replicates was 100%, with no differences between inter- and intra-run replicates observed.