BagGMM: Calling copy number variation by bagging multiple Gaussian mixture models from tumor and matched normal next-generation sequencing data

Abstract

Copy number variations (CNVs) contribute significantly to human genomic variability, some of which lead to diseases. However, effective detection of CNVs from whole genome next generation sequencing data (NGS) remains challenging. Here, we present BagGMM, a new method to call CNVs using tumor-normal matched samples from NGS data. BagGMM extracts read depth ratios of tumor samples to normal samples, divides the genomic sequences into segments by sliding windows to count the average coverage ratio of each segment, filters candidate deletions and duplications based on a coarse criterion of coverage ratio, and then builds Gaussian Mixture Model (GMM) for remaining ratios to identify the remaining ambiguous copy number states after filtration. Bagging multiple GMMs makes false positive calls descent instead of using one GMM, thus enhancing the detection power of BagGMM. Considering the computation speed of GMMs and false positive calls, we employ a segmentation procedure “large window and then small windows”, which is also helpful to determine boundary of CNV regions. We apply BagGMM to three simulation datasets and two groups of human whole genome sequencing (WGS) data for breast cancer patients and ovarian cancer patients to identify CNVs, respectively. All performed experiments demonstrate that BagGMM has the capability of robustly identification of CNVs with different sizes and states. The performance of this tool is compared to four peer existing CNV detection methods. BagGMM shows a significant improvement in both sensitivity and specificity for detecting both copy number gains and losses.

Introduction

Copy number variation (CNV) is defined as a gain or a loss of genome sequences where the DNA copy number deviates from the normal copy number. CNV, as a type of structural variation that plays an important role in human disease, is of widespread concern [1]. Biomedical research has estimated that a substantial proportion of the human genome consists of CNVs and more than 1000 CNV regions with a frequency above 1% have been identified in the genome [2]. Most of the CNVs have been reported to substantially contribute to phenotypic variations and common diseases by disrupting genes, altering gene dosage, and perturbing their expression levels [3], [4]. Some CNV encompassing genes also have been demonstrated to foster activation of oncogenes or inactivation of tumor suppressors in cancer cells [5], [6]. Therefore, specific focus is drawn on accurate inference of the interest and importance of CNVs in identifying disease-causing regulatory variants and genes.

With the breakthrough technology of high-throughput platforms, next generation sequencing (NGS) technology has been widely applied to measure the relative copy number between patients and normal counterparts, over a set of genomic regions, which achieves both high resolution and accuracy [7]. As a result, whole-genome sequencing (WGS) [8], whole-exome sequencing (WES) [9], and targeted capture sequencing from NGS technologies have become primary strategies for CNV detection and for studying the relationship between these variations and human diseases [10]. For each strategy, many CNV detection tools have been developed correspondingly. CNV detection tools based on targeted capture sequencing are designed by the capture, which is relatively efficient for CNV detection in novel regions. However the performance of this kind of tools is currently limited by the capture design. In contrast, CNV detection tools based on WGS are more effective for detecting CNVs in novel regions than those tools based on target capture sequencing in spite of their higher cost. In this work, all analyses of CNVs are based on WGS data.

Methods existed for the detection of CNVs from NGS data have been reported for whole genome sequences [11], [12]. In general, these can be grouped into four categories: (1) read depth (RD), (2) paired-end, (3) split-read, and (4) assembly. The underlying assumption of RD-based methods is that the number of sequencing reads aligns to a location of the genome is proportional to the number of copies at that location. Furthermore, the copy number of this region is estimated by counting the number of reads aligned to this region. Under this assumption, RD-based methods identify somatic CNVs by comparing the differences of RD in particular genomic regions between case and matched control samples. Methods in this category include mHMM [13], SeqCNV [2], Control-FREEC [14], and XCAVATOR [3]. For paired-end mapping methods, insertions or deletions are detected by comparing the distance between mapped end pairs on a reference genome and the known insert size [15]. This type of approach is mostly used for identifying other types of structural variations (beyond CNVs) such as inversion and translocation. Split-read based methods detect CNVs by finding and evaluating gaps in sequence alignment. In the assembly approach, short reads are used to assemble genomics regions by connecting overlapped short reads or contigs. CNVs are detected by comparing the assembled contigs to the reference genome. Compared to other above-mentioned CNV detection methods, RD-based methods achieve a higher precision solution at a lower cost [16], [17] for CNV detection. For this reason, RD-based methods have recently become a major approach for the identification CNVs.

RD-based methods for the detection of CNVs are based on the assumption that the sequencing process is uniform. Consequently, the number of reads mapped to a region (the total coverage of a region) is expected to be proportional to the number of times that the region appears in the DNA genome [18]. Following this idea, the absolute copy number can be inferred from the read counts ratio of diseased samples to matched normal samples.

Although the aforementioned RD-based methods are capable to detect CNVs with high accuracy, many studies have shown that each RD-based approach detects variation events with specific structural characteristics [19], [20]. Each method has its advantages; however, its detection power and sensitivity are still limited by low coverage data. Control-FREEC segments and normalizes read counts, and then provides an estimate of polymerization to detect CNV break points. Nevertheless, Control-FREEC is more sensitive to copy gains than copy loss inferring from the result of following experiments. m-HMM uses a K-means method to form windows by joining adjacent genomic sites and then establishes a Hidden Markov Model to discriminate copy number states, which has capability to detect CNVs with high sensitivity but high false positive rate. SeqCNV builds a statistical method, using maximum penalized likelihood estimation to evaluate CNV boundaries and copy number ratios. However, it just provides the types of copy number (gains or loss) while not providing the copies of CNV events. The detection power of SeqCNV and its sensitivity are also limited by low coverage depth data. XCAVATOR builds a normalization procedure for read counts using the statistical properties and biases of read counts, and combines SLM and FastCall algorithms to detect CNV regions. Being capable to discriminate five discrete copy number states (two- and single-copy deletions, neutral, three- and multiple-copies duplication), XCAVATOR detects duplications with lower sensitivity than deletions calls. When the variance of RD signals is large, XCAVATOR shows a lack of sensitivity for detecting signal shifts.

Currently available tools, except for XCAVATOR, are able to classify three states (deletion, normal, and amplification) without the ability of discriminating homozygous deletions (single-copy deletions) from heterozygous deletions (two-copies deletions). This keeps the use of sequencing data limited for the prediction of exact DNA copies.

To overcome this limitation of discriminating homozygous deletions from heterozygous deletions with high specificity and recalls, we have developed BagGMM, which is an analytic toolset that uses bagging of multiple Gaussian mixture models to detect CNVs from a pair of case–control genome sequencing data. In this method, we start to divide genomic sequences into segments by non-overlapping sliding windows according to genomic coordinates. The underlying assumption is that genomic loci within a segment share the same copy number state. For each segment, we divide its read count in the case sample by that of the control sample as coverage ratio of that segment. Here, we use a Gaussian assumption that all coverage ratio values in a certain copy number state (e.g., homozygous deletion, heterozygous deletion, normal, and duplication) in a genome follow a Gaussian distribution. For example, if there are three kinds of copy number status, the coverage ratios will be clustered into three different distributions. Therefore, there are multiple Gaussian distribution models of coverage ratios when a genome has different copy number status. Specially, we employ a segmentation procedure: “large windows with small windows” to divide genomic sequences into segments with a focus on reducing the computing time and determining CNV boundaries. Then, we designate segments as distinct copy number status via a coarse criterion of coverage ratio and leave segments with ambiguous copy number status. As a further step, based on the above-mentioned assumption, we build a Gaussian Mixture Model (GMM) for coverage ratios of remaining segments with ambiguous copy number states to identify CNVs, which models different Gaussian kernels to represent different copy number distributions. However, there are many potential false positive events in the GMM. To solve such a problem, we employ a bagging idea, implying that multiple GMMs are built by repeatedly re-sampling coverage ratios. In summary, BagGMM is suited for the detection of CNVs with high specificity. The code implementing our method is publicly available and can be accessed at https://github.com/tudui123/BagGMM.

To demonstrate the power of our method, we use BagGMM to analyze simulation datasets and real WGS datasets and compare its performance to that of four other state-of-the-art tools. All performed analyses demonstrate that our computational pipeline is capable to detect CNVs in WGS data outperforming all other compared tools.