Accurate and exact CNV identification from targeted high-throughput sequence data
© Nord et al; licensee BioMed Central Ltd. 2011
Received: 19 October 2010
Accepted: 12 April 2011
Published: 12 April 2011
Skip to main content
© Nord et al; licensee BioMed Central Ltd. 2011
Received: 19 October 2010
Accepted: 12 April 2011
Published: 12 April 2011
Massively parallel sequencing of barcoded DNA samples significantly increases screening efficiency for clinically important genes. Short read aligners are well suited to single nucleotide and indel detection. However, methods for CNV detection from targeted enrichment are lacking. We present a method combining coverage with map information for the identification of deletions and duplications in targeted sequence data.
Sequencing data is first scanned for gains and losses using a comparison of normalized coverage data between samples. CNV calls are confirmed by testing for a signature of sequences that span the CNV breakpoint. With our method, CNVs can be identified regardless of whether breakpoints are within regions targeted for sequencing. For CNVs where at least one breakpoint is within targeted sequence, exact CNV breakpoints can be identified. In a test data set of 96 subjects sequenced across ~1 Mb genomic sequence using multiplexing technology, our method detected mutations as small as 31 bp, predicted quantitative copy count, and had a low false-positive rate.
Application of this method allows for identification of gains and losses in targeted sequence data, providing comprehensive mutation screening when combined with a short read aligner.
Massively parallel sequencing technology can be used to efficiently interrogate multiple targeted genomic regions for clinically relevant mutations . The use of short unique sequence indexes, or barcodes, increases sequencing throughput and allows analysis of multiple samples in a single sequencing run. Application of this technology will allow cost-efficient screening of genes known or suspected to harbor clinically relevant pathogenic mutations. While small mutations are readily detectable with current analysis platforms [2, 3], methods for the identification of copy-number variants (CNVs) are not well established for targeted data. The ability to detect the full spectrum of mutations is critical to the success of targeted sequencing projects.
Current methods for detecting structural variation from massively parallel data use either paired-end mapping or depth of coverage methods (see  for a recent review). Paired-end mapping methods perform well for mate-pair sequence data and have the advantage of identifying both balanced (e.g. translocations or inversions) and unbalanced (deletions and duplications) structural variation. Comprehensive CNV scanning using paired-end mapping requires mulitple insert sizes and the presence of paired tags that bridge the breakpoint of the event. This approach is widely used in the assembly of whole genomes but cannot be applied to targeted enrichment protocols, as hybridization-based capture methods for massively parallel sequencing typically require small insert fragments (200 bp) and breakpoints of structural variants may lie outside of targeted regions.
Depth of coverage methods do not require breakpoint capture or have specific size restrictions, yet such methods are hampered by low signal to noise ratio due to sparse sampling (a product of low overall coverage) and sampling bias based on sequence GC-content, bait coverage, and other factors affecting capture efficiency. Further, depth of coverage methods require cross-sample normalization and comparison in order to deal with variation in coverage across genomic regions and samples. We previously showed that a basic depth of coverage approach performs well to identify CNVs from high-coverage (~1000×) targeted sequence data. Depth of coverage methods have achieved higher signal to noise by averaging data across a genomic region[6–8], however, doing so results in a loss of the ability to detect smaller mutations. In addition, depth of coverage methods are poor at localizing breakpoints .
Here we describe a method that combines CNV detection through high-resolution depth of coverage analysis with call confirmation through partially-mapped reads. This allows for the identification of DNA gain or loss even where breakpoints are outside the targeted regions, and exact breakpoint characterization where breakpoints are within the targeted regions. By combining these orthologous approaches, CNVs of any size can be identified from targeted high-throughput sequence data with a low false-positive rate.
Our method uses two independent approaches, depth of coverage and a scan for partially-mapped reads at CNV edges. We implemented a straightforward depth of coverage algorithm and internally validated CNV calls by looking for a signature of partially-mapped reads that confirm the CNV breakpoints where possible. Partially-mapped reads are defined as high-quality full length (76 bp) reads where the best alignment to targeted regions is a perfect match less than 76 bases in length. A minimum match length of 30 bp was used for this analysis. Where partially-mapped reads match to a CNV edge, it is assumed that the unaligned portion of the read maps to sequence flanking the other breakpoint and the partially-mapping reads can be aligned to each other to form a consensus sequence that flanks the breakpoint. For example, in the case of a deletion, the partially mapped reads will correspond to sequence flanking the deleted region, giving the exact breakpoints of the CNV. Partially-mapped reads map to the reference sequence differently for gains and losses, as described for paired-end mapping methods of CNV detection . Where CNV breakpoints lie within non-targeted sequence, partially-mapped reads are unavailable and depth of coverage is the only evidence used for CNV calling.
We obtained median raw coverage of 350× across the 909 kbp of targeted sequence (Additional file 1, Figure S1). Low coverage regions were primarily GC-rich, consistent with other targeted sequence data, and account for the majority of bases excluded from screening for structural variation. 99.8% of targeted bases had median coverage of at least 10×, permitting high-confidence detection of single base and short indels for nearly all sequence assayed. In order to estimate the effect of random noise on the false-negative rate, we simulated data where one CNV of 50 bp, 100 bp, 200 bp, 500 bp, or 1000 bp was present within 1 mb of sequence and tested for detection of the CNV across varying signal-to-noise ratios (Additional file 2, Figure S2). This analysis indicated that our methods can detect a 200 bp CNV with 87% sensitivity at a signal-to-noise ratio of 6, a criteria met by >92.3% of the bases targeted in this experiment, and a 100 bp CNV with 80% sensitivity at signal-to-noise ratio of 8, a criteria met by 75.8% of the bases.
CNVs identified in 94 subjects across ~1 Mb targeted sequence
Median Sample Coverage
Median Lane Coverage
DoC + Seq
DoC + Seq
DoC + Seq
DoC + Seq
The first deletion, of 510 bp, removes an exon of BRCA1 (chr17:41,200,740-41,201,249). We confirmed a 31 bp homozygous deletion in BRCA2 (chr13:32,960,705-32,960,735) that matches a known indel (rs56213495), indicating that our method is sensitive to small CNVs. The final two deletions mapped to the same PTEN coordinates and overlap with a known CNV previously identified by array comparative genomic hybridization (aCGH) (DGV ID: 657536). In comparison to aCGH methods, we were able to resolve exact breakpoints for this 899 bp polymorphic deletion (chr10:89,652,824-89,653,723). Importantly, mapping sequence tags to the breakpoints allowed exclusion of a large number of false positive calls. All seven CNVs with breakpoints predicted to be within repetitive non-targeted sequence were validated by Sanger sequencing or MLPA. In the 909 kbp in 94 samples there was one CNV called that could not be confirmed. This 109 bp duplication could not be verified through sequence mapping or Sanger sequencing and is a likely false positive. Breakpoint resolution via partial read mapping reduced the number of false positives to a single duplication across 94 samples and 909 kbp.
Targeted high-thoughput sequencing data is used to produce high-coverage (often >50×) sequence data on specific regions of the genome. Targeted regions are often non-contiguous, as in the case of coding sequence capture or interruptions due to repeat sequences. We created a method that is specific to high-coverage, non-contiguous sequence data for robust CNV detection from targeted sequence data. Our method does not require sequencing across CNV breakpoints (as paired-end methods do), but leverages sequence level data for exact CNV characterization where breakpoints are within targeted regions. Using a test dataset of regions that contribute to inherited breast and ovarian cancer susceptibility, we identified 10 mutations (7 known to be pathogenic and 3 benign), localizing 4 mutations to exact genomic breakpoints. The minimum size for a CNV that we detected was 31 bp, and there was only one suspected false positive among all samples that passed quality control.
Our method has several advantages over previous methods used to detect CNVs from high-throughput sequencing. While sequence-based CNV signal (e.g. partially-mapped reads or split paired-end reads) is preferable over relative coverage, sequence data is unavailable in targeted data whenever the CNV breakpoints are not in the targeted region. This introduces a major issue for methods that rely on sequence-based data for CNV detection (e.g. [11, 12]). Additionally, as described earlier, paired-end mapping strategies are not currently applicable when targeted enrichment is used. In comparison, relative depth of coverage is problematic due to local variation in coverage within samples and variation in coverage across samples. This leads to problems with false-positive signal when scanning for regions that depart from expected diploid coverage. This issue can be partially addressed by averaging read depth across a genomic region (e.g. [6–8]), however, doing so reduces the ability to detect small mutations. By combining ortholgous algorithms, we eliminated from consideration a large number of false-positive CNVs called by depth of coverage methods where no signal of partially-mapped reads was present. The ability to computationally resolve the majority of false-positive calls permitted us to use high-resolution analysis even where signal-to-noise ratio was intermediate. In addition, depth of coverage analysis has the benefit of producing a quantitative estimate of copy-count, as demonstrated by the clear distinction between duplication and triplication ratios.
The major limitation of depth of coverage methods arises in duplicated genomic regions, such as segmental duplications or smaller regions of homology between genomic regions. Within such regions, the expected ratio for deletions and duplications is relative to the copy count in the genome. For depth of coverage, a copy number change from two to one (hemizygous deletion) or two to three (duplication) produces a strong signal, however, copy number changes in duplicated regions of the genome produce weaker signal that can be easily overwhelmed by noise. For example, variation between samples will overwhelm the signal of a change from six copies to five, where the expected ratio would be 0.83. As such, we limit analysis to unique regions of the genome. In addition, depth of coverage methods do not detect balanced structural variation.
The methods we describe here were optimized for multiplexed high-coverage data, but we expect the fundamental approach will be applicable to experiments with a larger volume of targeted sequence, such as exome capture. While we used multiplexed data from a relatively small number of targeted bases, the method we describe is scalable and could be used on any dataset where signal-to-noise ratio allows reliable depth of coverage CNV calling. Where local variation is too high or coverage too low, average coverage across a genomic window can be used instead of single base coverage, although this lowers resolution to the size of the window used. However, scanning across the window for partially-mapped reads as we do here would still permit exact CNV breakpoint characterization when breakpoints are in targeted sequence.
Targeted massively parallel sequencing permits screening of genomic regions for multiple samples simultaneously, and thus is a powerful and cost-effective tool for characterizing mutations in contexts where multiple genes or pathways are involved[1, 13–15]. We demonstrate that a combination of depth of coverage and analysis of partially-mapped reads allows detection with high-confidence of CNVs within unique sequence, to a resolution of 31 bp. Furthermore, we demonstrate that ratio can be used to estimate absolute copy number and that depth of coverage alone is sufficient for CNVs where breakpoints are not represented within targeted regions. The method we describe here is robust, yet simple to implement; while we employ specific mapping and segmentation algorithms, the overall framework is amenable to improvement and variation in the algorithms as massively parallel sequencing technology progresses.
Sample selection and DNA extraction were as described previously. Library preparation, sequencing, and mapping of reads to the genome performed as described per Illumina and Agilent protocols http://www.illumina.com and http://www.agilent.com. Three micrograms of high quality genomic DNA was sonicated on a Covaris S2 to a peak size of 150-200 bp. DNA was end-repaired, A-tailed, ligated to standard Illumina adapters and amplified with flanking primers. 500 ng of paired-end library was hybridized individually to a custom design of cRNA oligonucleotides (Agilent, ELID:0279281, BED file available upon request) and amplified with index specific primers from the Illumina Mutiplex Oligonucleotide kit. Equal concentrations (1.3 pM) of 12 paired end indexed libraries were pooled per lane on a v4 flow cell and subjected to cluster amplification on a cBot instrument. Raw sequence coverage was determined by tabulating the number of tags that overlapped each targeted base as mapped by the standard Illumina GAIIx bioinformatics pipeline with the demultiplexing script http://www.illumina.com. Coverage data for the targeted regions is shown in Additional file 3, Table S1.
The computational methods described here can be separated into four steps: raw coverage normalization, correction for sample-specific coverage biases, CNV calling, and partially-mapped read analysis. All analysis was done using the R platform http://www.r-project.org with custom scripts, which are available at request. We emphasize that the methods described below are one specific implementation within a general framework. The approach we describe is amenable to alternative normalization methods, application of any segmentation algorithm that can employ ratio data, and use of other sequence mapping algorithms suitable for short read data.
Single sample coverage was normalized to the other eleven samples run in the same lane. Normalization was performed using invariant set methods. Briefly, bases that had a similar coverage rank in the sample and median data were selected and a line was fit using a smooth spline function across the range of coverage counts. Normalization reduced the median standard deviation from 85.7 for raw coverage to 27.1 after normalization (Additional file 1, Figure S1).
where d is the distance from the first targeted base and the constant 10 is used to account for off-target capture where the first captured base is outside of the region of hybridization.
Normalized coverage was corrected using simple linear regression with the transformed GC-content and bait capture likelihood terms as the independent variables and log(normalized coverage) as the dependent variable. Bases with coverage < 50 were excluded from contributing during curve fitting, as these bases show high levels of variation due to sparse sampling effects. Signal to noise ratio (mean/standard deviation) went from 12.3 for uncorrected normalized coverage, to 13.2 after correction (Additional file 1, Figure S1). 99% of bases had corrected coverage signal to noise ratio > 5 and 83% had signal to noise ratio > 10. Bases with signal to noise < 2 were excluded from CNV analysis.
Ratio data was calculated by comparing sample corrected coverage to the median corrected coverage of the other 11 samples in the lane. We used a sliding window to identify regions where the majority of bases had a ratio < = 0.6, indicating copy number loss, or > = 1.4, indicating copy number gain. Window size was 20 bp, with a minimum of 18 bp required to meet the criteria for either gain or loss. Copy-number variant regions were generated by extending variant windows and merging neighboring variant regions. The median z-score across the CNV region was calculated by comparing sample ratio to mean and standard deviation of the remaining 11 samples in the lane. Z-scores closer to zero could reflect either diploid regions with high variation in coverage across samples or regions with an underlying common copy number polymorphism. As we were interested in rare mutations, CNV regions were required to meet the criteria of median absolute z-score > = 5 and signal to noise ratio > = 9 for acceptance. These values were set based on the minimum z-score and signal to noise ratio for CNVs that could be confirmed with Sanger sequencing.
Across diploid regions (where no CNV was both called and confirmed), 98.9% of the ratio values were between .6 and 1.4 (Additional file 1, Figure S1), compared to only 9.3% for confirmed CNV regions. The majority of bases that had ratio values < = 0.6 or > = 1.4 represented false positive signal, as no CNV could be confirmed. However, these false positive bases showed low signal to noise and z-scores were closer to zero versus true positive bases (Additional file 4, Figure S3). For ratio values for regions where a deletion was present and ratio values were < 0.6, 89.2% had signal to noise values > 10 and 91.1% had z-score values < -5. For regions with a confirmed duplication and ratio values > = 1.4, 78.6% had signal to noise values > 10 and 91% had z-scores > 5. The differences in signal to noise and z-score between confirmed CNV bases and false positive CNV bases indicates that ratio is a stronger indicator of true copy number when signal variation is taken into account.
CNVs identified using depth of coverage methods were then examined for the presence of partially-mapped reads within the predicted breakpoint region. This is only possible for CNVs where at least one breakpoint is within targeted sequence. CNVs where both breakpoints appear to be in non-targeted regions (e.g. repeat sequence or sequence outside the targeted region boundaries) were excluded from this step. Repeat regions often mediate CNV formation[18–20] and minimum contiguous targeted region size was limited by bait probe size to 80 bp in this experiment, the majority of calls where both breakpoints are in repeats were confirmed by Sanger sequencing (7/8) and all calls were >100 bp in length, indicating that such calls are likely to be true positives. For calls where at least one breakpoint is unique we used BLAST to test for partial reads that aligned to the specific breakpoint region. Partial reads were 76 bp in length and mapped at least 30 bp to the region. These parameters were selected as 76 bp is a full length read, and 30 bp was set as minimum cut-off to avoid mapping repeat sequence fragments. These lengths are dependent on platform read length and should change accordingly. Mean and standard deviation were calculated for the number of partial reads across the targeted breakpoint region, and z-score for each base was calculated and used to detect breakpoint based on over-representation of tags with the same start/stop position. Where no over-representation of start/end sites for partial reads was present, CNV calls were considered to be false positives. If a signature of partial reads was present, we mapped the consensus sequence of aligned partial reads to the genome to identify the exact CNV breakpoints.
We thank Anne M. Thornton for technical advice and assistance. This work was supported by the Breast Cancer Research Foundation, the Susan G. Komen Foundation for the Cure, and the Department of Defense Congressionally Directed Medical Research Programs for Breast Cancer and Ovarian Cancer. M.-C.K. is an American Cancer Society Research Professor. A.S.N. is supported by NIH grant T32 ES015459.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.