CoNVEX: copy number variation estimation in exome sequencing data using HMM
© Amarasinghe et al.; licensee BioMed Central Ltd. 2013
Published: 21 January 2013
Skip to main content
© Amarasinghe et al.; licensee BioMed Central Ltd. 2013
Published: 21 January 2013
One of the main types of genetic variations in cancer is Copy Number Variations (CNV). Whole exome sequenicng (WES) is a popular alternative to whole genome sequencing (WGS) to study disease specific genomic variations. However, finding CNV in Cancer samples using WES data has not been fully explored.
We present a new method, called CoNVEX, to estimate copy number variation in whole exome sequencing data. It uses ratio of tumour and matched normal average read depths at each exonic region, to predict the copy gain or loss. The useful signal produced by WES data will be hindered by the intrinsic noise present in the data itself. This limits its capacity to be used as a highly reliable CNV detection source. Here, we propose a method that consists of discrete wavelet transform (DWT) to reduce noise. The identification of copy number gains/losses of each targeted region is performed by a Hidden Markov Model (HMM).
HMM is frequently used to identify CNV in data produced by various technologies including Array Comparative Genomic Hybridization (aCGH) and WGS. Here, we propose an HMM to detect CNV in cancer exome data. We used modified data from 1000 Genomes project to evaluate the performance of the proposed method. Using these data we have shown that CoNVEX outperforms the existing methods significantly in terms of precision. Overall, CoNVEX achieved a sensitivity of more than 92% and a precision of more than 50%.
Commercial products of Next Generation Sequencing (NGS) Technologies such as Roche/454 FLX, Illumina Genome Analyzer/HiSeq, Applied Biosystems SOLiD™System and Helicos Heliscope™have enabled the sequencing of DNA much faster and cheaper than before . These have shifted the paradigm of biological sequence analysis to a new level. Currently these are not only being used for the sequencing of whole genome, but also for sequencing of known exons and transcriptomes as well. The main motivations behind the technology of targeted resequencing (TR) include the following among others. The actual coding regions or the exons of the human genome account only for ~1% of the total sequences, which consequently gives about 30 Mb data compared to 3 Gb data in WGS . Currently, getting higher coverage of targeted regions using NGS technologies is about six times  cheaper and faster compared to achieving the same coverage of whole genome. On the other hand approximately 85% of disease causing mutations lie in the coding regions . Targeted resequencing has been mainly used in medical sequencing to find disease causing genetic variations (a review can be found in ). Recent studies on TR and WES data have successfully detected cancer specific mutations (somatic mutations) in breast cancer [6–8], ovarian cancer  and prostate cancer . Although, exome sequencing has been successfully used to find small variations in cancer genomes, its potential to find large structural variations such as CNV has not yet been fully explored.
Cancer arises due to the acquisition of many somatic variations by the DNA of cancer cells . Copy Number Alterations (CNA) play a major part in the progression of this deadly disease . Until recently, the most common method to detect CNV in cancer DNA was to use micro array based technologies. However, during last 2 - 3 years many algorithms have been developed to identify CNV in cancerous data generated by whole genome sequencing [11–15], making use of the vast amount of data produced by NGS technology. The higher resolution that can be achieved through NGS data has helped to detect new variations that were undetectable previously and include CNVs which are as small as 50 bp . These methods use the number of reads mapped to a particular region in the genome, to find copy number varying regions in one genome compared to one or more other genomes. Some of these methods have been adapted from the methods used in aCGH. For example Circular Binary Segmentation (CBS)  and Hidden Markov Model (HMM) . However, methods in whole genome sequencing cannot be directly applied to whole exome sequencing data due to the small size and sparseness of these data . On the other hand, the useful signal will be hindered by the intrinsic noise present in the exome sequencing data itself due to various biases introduced in target capturing and sequencing phases. To address these issues and to utilize the advantages provided by targeted resequencing, new algorithms have to be developed. Since the end of 2011, very few bioinformatics methods for detecting copy number variations in targeted resequencing data have been published. The method in  describes the use of TR data to detect CNV in cancer samples. However, the targeted regions in this method are larger in size (~ 40 kb), where as exons are much shorter (200 bp -300 bp). Methods such as [21, 22] are developed to find CNV in non cancerous exome data, such as in population studies. CONTRA  is a recent method proposed to evaluate cancer TR data using a pooled or a matched normal sample. ExomeCNV  and Var Scan 2  are specifically designed for whole exome sequencing of cancer samples. A limitation in these approaches is that they have a higher number of false positives which result in a very low precision (further discussed in Results and discussion section for ExomeCNV).
In this work, we present CoNVEX, a method that evaluates exon level depth of coverage ratios to assess variation in copy number of whole exome capture data produced from cancer samples. We propose to use Discrete Wavelet Transformation denoising to reduce the variability of coverage ratios and then use HMM to detect copy number variations. Our method reduces the number of false positives by efficient pre-processing of the data, which results in a mean precision of more than 50%.
Where and are the mean normalized DOC of tumour and control respectively.
Here, is the true signal of copy number variation with additive noise, ϵ i . This noise can be assumed to be iid with N(0, σ 2)where σ is the standard deviation of the distribution. We have used DWT smoothing  on R i , to detect true signal to increases the ability of actual copy number prediction. The DWT smoothing procedure starts by first taking discrete wavelet transformation of ratios using "HAAR" wavelet. The fundamental assumption behind discrete wavelet transform is that, there is a correlation between the two neighbouring samples or data points. This is very much true in predicting CNVs as they span multiple successive exons. The selection of HAAR wavelet family was based on the fact that it computes the wavelet coefficients as the difference between two near by blocks of data points. This feature helps to retain the information regarding copy number aberration points. The shrinking of the DWT coefficients were done using soft thresholding function and the threshold value was calculated by Stein's unbiased risk estimator (SURE) for each level of DWT. Finally, the modified coefficients were used to reconstruct the de-noised signal at i th location of chromosome j, , by taking the inverse transform.
The total number of hidden states in the model is given by K and those are denoted by S = S 1, S 2,..., S K . If there are L exons in the sample of consideration, the state of l th exon (e l ) equals to S k where 1 ≤ l ≤ L and 1 ≤ k ≤ K.
Here, represents the Gaussian distribution. Mean (μ k ) of that distribution vary with different states and the provided normal cell contamination percentage and ploidy. We used a common standard deviation, σ, to all states.
The above HMM can be represented compactly as λ = (A, B, π) where A, B and π represent transition probability matrix, emission probability distribution and initial state distribution. When fitting the above HMM, the K states must be fixed at first and normal contamination and tumour ploidy must be given as inputs.
The optimal λ is selected by optimizing the negative log-likelihood [27, 28]. The initial state distribution π is chosen such that higher probability is attached to the most abundantly expected state or the normal state (i.e. copy 2 in normal humans). Similarly, the transition probability matrix A, is chosen such that, a higher probability is assigned to remain in the same state and lower probability is assigned to transition to another state. Also the transition to normal state has higher probability than transition to a CNV state. Then we used Viterbi Algorithm to assign the most appropriate copy number state for each exon.
where P T is the most abundant ploidy in the tumour sample.
We randomly selected six samples, NA18536, NA18543, NA18544, NA18548, NA18557, NA18558, from 1000 Genome project, which share some common attributes, to evaluate the performance of the proposed method. These selected individuals have been studied by the HapMap project http://www.hapmap.org. The common features in these individuals are (i) exome sequencing was performed by the Beijing Genome Institute, hence a common exome-capture (NimbleGen V2) has been performed, (ii) male individuals and (iii) from CHB population.
First, we retain only the copy number neutral regions in each sample. The CNV information were downloaded from the HapMap project genotype file.
We selected one sample (NA18536) as the Control and others as the tumour with known CNV.
To do the simulation of gains and losses we randomly selected a region in the Chr 1 and reduce (e.g: multiplied by 0.05) or amplified (e.g: multiplied by 2) the number of reads in that particular region. For each variation type, we perform 100 simulations.
When we evaluated the performance using only one sample (NA18543), we used 100 simulations for each variation type. When we used 5 samples in simulations, 20 variations were simulated in each individual sample.
To incorporate contamination in the simulation, we mix the control sample and simulated sample as per the relationship (α∗ Control + (1 - α) ∗ Tumor) where α is the contamination proportion.
Different methods have been proposed to reduce the experimental biases present in TR data. These include GC content bias reduction using regression methods [20, 21], taking base level ratios between normal and control samples [19, 23] and bait capture bias reduction using log transformation . Those methods, adapted from aCGH or whole genome sequencing based approaches, try to reduce different experimental biases separately. Hsu et. al.  proposed DWT smoothing as an effective method to extract true copy number variations from aCGH data.
In this work, we propose to combine the strengths of both DWT and HMM to robustly predict copy number variations in cancer samples. The main novelty of our approach is the use of DWT smoothing to reduce experimental biases present in whole exome sequencing data prior to applying a Hidden Markov Model. These experimental biases are modelled here as additive noise to the true signal. The wavelet coefficients, which are the differences between two nearby data blocks, can be used to reduce noise. This is achieved through approximating some coefficients that do not by pass a certain threshold to zero. After thresholding step when the inverse transform is performed on these wavelet coefficients, we can generate a smoother version of the input signal. Exon level ratios, before and after DWT smoothing, for data downloaded from 1000 Genome project (http://www.1000genomes.org) are given in Figure 1.
After smoothing, we applied an HMM described in Methods section to detect copy gains and losses. Hidden Markov Models have been previously used to detect CNV in exome data (an R package called ExomeCopy) , but not used in this manner to detect CNV in tumour samples. The differences between ExomeCopy and CoNVEX are,
ExomeCopy uses HMM to identify CNVs in male patients with X-linked Intellectual Disabilities (XLID)
They have used depth of coverage of exons as observations or emissions of hidden states
The robustness in copy identification is achieved by pooling coverage data from all patients
Therefore, it fails to identify relative copy number in cancer samples against a matched normal.
We carried out a comparison between the proposed method and the existing method, ExomeCNV . Using simulated data, we were able to assess the performance of CoNVEX and ExomeCNV for different size ranges.
A true positive (TP) is identified when the gain or loss of an exon is correctly identified by the algorithm and a false positive (FP) identification is defined in the same manner. When using ExomeCNV, we used their primary CNV detection method (here after referred to as ExomeCNV1) and the extension which combines DNACopy  (here after referred to as ExomeCNV2) separately on our simulated data sets. The DNACopy version of ExomeCNV is applied to make sure that we get results for all exons that pass the default cut-off level of the coverage (this is a direction given by the authors of the paper). We used default parameters given in ExomeCNV R package for CNV prediction, except for read length and admixture rate, which we set to 90 and 0.0 in our evaluation.
Performance of proposed method for 100 simulations.
Deletions (1 k -1 M bp)
97.82 ± 12.37%
99.94 ± 0.081%
79.25 ± 23.23%
99.94 ± 0.081%
Duplications (1 k -1 M bp)
95.25 ± 19.64%
99.93 ± 0.082%
77.04 ± 26.43%
99.93 ± 0.085%
Performance of ExomeCNV1 for 100 simulations.
Deletions (1 k - 1 M bp)
97.91 ± 2.81%
86.20 ± 1.57%
8.76 ± 6.54%
86.24 ± 1.56%
Duplications (1 k - 1 M bp)
90.68 ± 9.02%
86.26 ± 1.55%
8.96 ± 8.57%
86.28 ± 1.54%
Performance of ExomeCNV2 for 100 simulations.
Deletions (1 k - 1 M bp)
99.26 ± 2.11%
96.00 ± 1.67%
8.69 ± 6.50%
96.01 ± 1.66%
Duplications (1 k - 1 M bp)
99.98 ± 0.16%
96.06 ± 1.65%
9.62 ± 9.25%
96.08 ± 1.64%
When compared with ExomeCNV2, our method showed superior performance in terms of specificity, precision and accuracy. Slight decrease in sensitivity was observed in CoNVEX, this is mainly due to the detecting short variations involving 1 or 2 exons. This can be attributed to the smoothing step we performed using DWT. Because of this we separately tested the performance of CoNVEX for shorter variations sizes as described below. Both versions of ExomeCNV, showed very poor performance when it comes to precision, as it tries to detect as many as possible variations to maintain a higher sensitivity rate.
Performance of CoNVEX against other methods.
ExomeCopy and CONTRA did not identify any of the variations present in the test sample. This can be attributed to the fact that these are specifically designed for using a background sample  or a robust baseline . VarScan 2 was able to identify the hemizygous duplication in the region with 60% sensitivity, however the number of false positives reported by the method was very high (false positive rate of 32%). CoNVEX performed well with 90% sensitivity and 0.06% false positive rate.
Median sensitivity of CoNVEX for small variation detection is 100%. Every deletion of size, more than 200 bp was detected by our method. Hence, giving a mean sensitivity of 100% for detecting deletions. Mean sensitivity of detecting each duplication size was more than 85%. As seen in the graph, almost every variation of size of more than 800 bases can be detected by the proposed method. Also, a median precision of more than 30% can be achieved.
Exome sequencing data can be used to detect copy number variations as an initial screening procedure. It is a cheap and time efficient method. We have successfully applied the proposed method on exome data to identify CNVs spanning one to thousands of exons. However, actual breakpoint of the CNV would not necessarily lie in the coding region. This limits the use of WES in identifying actual breakpoints of the CNV.
As discussed in the Results and Discussion section, we have achieved a higher precision than existing methods in detecting variations due to the data smoothing step. However, detection of some of the small variations may be missed by this smoothing step, as these can be recognised as noise. Further analysis is needed in order to better detect these variations among higher level of noise.
Although, we have used a matched normal sample to detect CNVs, the CNV identification can be done based on a pooled normal sample as described in . This might give an advantage in finding CNVs in familial studies assuming all members have a median copy number of two.
We thank Dr. Isaam Saeed and Dr. Suhinthan Maheswararajah for initial discussions on HMM. We used resources from both University of Melbourne and Peter MacCallum Cancer Centre for data processing and analysis. This work is partially funded by Australian Research Council (grant DP1096296).
This article has been published as part of BMC Bioinformatics Volume 14 Supplement 2, 2013: Selected articles from the Eleventh Asia Pacific Bioinformatics Conference (APBC 2013): Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S2.
The funding for open access charges were provided by The University of Melbourne.
This article is published under license to BioMed Central Ltd. 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.