Improving detection of copy-number variation by simultaneous bias correction and read-depth segmentation

Structural variation is an important class of genetic variation in mammals. High-throughput sequencing (HTS) technologies promise to revolutionize copy-number variation (CNV) detection but present substantial analytic challenges. Converging evidence suggests that multiple types of CNV-informative data (e.g. read-depth, read-pair, split-read) need be considered, and that sophisticated methods are needed for more accurate CNV detection. We observed that various sources of experimental biases in HTS confound read-depth estimation, and note that bias correction has not been adequately addressed by existing methods. We present a novel read-depth–based method, GENSENG, which uses a hidden Markov model and negative binomial regression framework to identify regions of discrete copy-number changes while simultaneously accounting for the effects of multiple confounders. Based on extensive calibration using multiple HTS data sets, we conclude that our method outperforms existing read-depth–based CNV detection algorithms. The concept of simultaneous bias correction and CNV detection can serve as a basis for combining read-depth with other types of information such as read-pair or split-read in a single analysis. A user-friendly and computationally efficient implementation of our method is freely available.


Improving detection of copy number variation by simultaneous bias correction and read-depth segmentation
Szatkiewicz et al.    (55)).

Model Deletions Duplications
Step We incorporated various techniques for bias correction, i.e. to account for GC content and mappability, to model the auto-regression in the read-depth data computed from overlapping windows, to account for additional noises by fitting a mixture of negative binomial and uniform distributions, and to apply quality control to prioritize the CNV calls. As these techniques incorporated incrementally, there are five partial versions of our algorithm as shown in Table S2.
In the first step the partial algorithm only corrects for GC content. In the next step, the partial algorithm corrects for GC content and mappability together. Similarly, the auto-regression and the mixture modeling are added incrementally. The algorithm becomes our full algorithm after the CNV quality control step is applied where we (1) removed predicted CNVs that are shorter than 800bps (i.e. those that appear in only one window), and (2)  Sensitivity is computed as the total number of high-confidence CNVs that were overlapped by the GENSENG predicted calls (>50% overlapping) divided by the total number of highconfidence CNVs. We then used the total number and total base pairs of the GENSENG predicted calls as a surrogate measure of specificity. As shown in Table S2, correcting for GC content alone resulted in low sensitivity and specificity. After various techniques were added Page 4 of 19 step-wisely, both the sensitivity and the specificity were improved. Correcting for mappability resulted in the most substantial improvement in both sensitivity and specificity. Quality control of the predicted CNVs substantially improved specificity with minimal loss in sensitivity. The full algorithm demonstrated the best performance and thus was used as the release of our GENSENG software. 1. The high-confidence dataset used in this comparison was generated by Handsaker et al (25) and was downloaded from ftp://ftp.broadinstitute.org/pub/svtoolkit/misc/1kg/NGPaper/. The coordinates of reported deletions were translated from NCBI36 to NCBI37 using liftOver.
2. Default filters were applied. CNVnator filter: the default q0 filters removes any predicted calls that have >50% reads with zero-valued MAPQ (i.e. reads with multiple mapping locations).
GENSENG filter removes any predicted calls that have RDA values lower than 0.75.
3. Sensitivity was calculated by dividing the number of the overlapping events (>1bp or >50% overlap with the high-confidence CNVs) by the total number of high-confidence CNVs.
4. The total numbers of deletions reported by the authors of CNVnator are 4100, 2223, 2352 for NA12878, NA12891, NA12892 respectively based on the same detection (i.e.100bp bin size) and quality control (i.e. default q0 filter) parameters as used in Table S3. These replicable results suggest that our usage of CNVnator is accurate. Note: 1. The released structural variation (SV) calls (4) were downloaded from ftp://ftpmouse.sanger.ac.uk/current_svs/. These SVs have been classified into several categories based on specific paired-end mapping patterns (4). We extracted 2 categories, including deletions and copy number gains (GAINS and TANDEMDUP), to compare to the GENSENG predicted calls.
2. This comparison focuses on the 19 autosomes.
3. GENSENG filter removes any predicted calls that have RDA values ranging between 0.5-2 Page 6 of 19 • Correct for GC content, mappability, and additional noises in the data.
• Use cubic spline smoothing to describe the GC-content vs read-depth relationship for each specific dataset (fitting the data most closely).
• Correct for GC-content only.
• Use a simple linear regression to describe the GC-content vs read-depth relationship (in most cases not appropriate).
• Use all reads  To identify the lower bound that GENSENG can handle, we applied GENSENG to data with varying sequencing coverage generated by simulation and compared the performance to that based on the native coverage (40x) using the same evaluation metrics (as was done in Table 1 of the main texts). In Table S6, we repeated the simulation process as described earlier with the targeted coverage been set as 5x, 10x, 20x, 30x, and 40x. To test the consistency of our simulation, we simulated data at 40x coverage for 10 times and observed replicable results (data not shown). As in Table 1, stringent GENSENG filter removes any predicted calls that span only one window, have RDA values ranging between 0.5-1.25 or mappablility < 0.3. The results demonstrate that: (1) higher sequencing coverage improves CNV detection power; (2) the lower bound of sequencing coverage that yield reasonably good performance of GENSENG is 10x (70% sensitivity and 6% FDR for deletions; 57% sensitivity and 0 FDR for duplications); (3) GENSENG could potentially work on data sequenced to as low as 5x but with much reduced sensitivity (fell by 14% for deletions and fell by 10% for duplications). To identify the lower bound that GENSENG can handle, we applied GENSENG to datasets down-sampled from 1000GP and compared the performance to that based on the native coverage (40x) using the same evaluation metrics (as was done in Table 2 of the main texts). In       Based on the input data computed for chromosome 22 of NA12878, we extracted compatible windows that have the same GC content class (0.4-0.5), the same mappability score class (0.9-1), and have high likelihood to be copy number normal regions (e.g. having not overlapped any high-confidence CNVs or other candidate CNVs). We then compared the observed read-depth distribution from these compatible windows to its theoretical expectations. Blue: observed distribution; red: Poisson distribution; green: negative binomial distribution. As shown in the figure, the Poisson distribution still does not fit the data; whereas the negative binomial distribution fits the data well. From the 716 predicted GENSENG calls (690 deletions and 26 duplications), a CNV-qualitycontrol filter was applied to remove any GENSENG predicted CNVs with mappability less than 0.3, yielding in 77 predicted deletions and 16 predicted duplications shown in this figure. In Figure S4a, the red dots indicate the true discoveries (29) and black the false discoveries (48).

Figure S4 RDA ranking and prioritization of CNV prediction from simulation study
It is evident that ranking the RDA values (the X-axis) correctly prioritized the predictions made by GENSENG. In Figure S4c, all 16 predicted duplications had RDA>1.25 and were all true discoveries. The GC content and mappability of the region are plotted in the fifth and the sixth panels

Figure S5 Example CNVs identified from the human HTS dataset
respectively. The last panel shows the locations of segmental duplication (purple) and repeatmask repetitive DNAs (orange). Note that beginning at approximately 3684100 bp, there appears to be a deletion-like artifact in all strains with zero-valued read-depth. Close inspection suggested that it was caused by zero-valued mappability and was not called as deletion by GENSENG because of its ability to correct biases.

Supplemental Methods
Supplemental Methods for "Improving detection of copy number variation by simultaneous bias correction and read-depth segmentation"

Computation of covariate values
The set of covariates include GC content and mappability score.
GC content is computed as in the following steps. (1) Calculate the proportion of G or C bases in each window from a given reference genome.
(2) Apply a cubic spline smoothing and then transform the GC proportion based on the fitted curve so that the transformed GC proportion and logarithm of the read-depth are linearly correlated. (3) The transformed GC proportion is median-centered and is referred to as GC content hereafter.
Mappability score is calculated a priori in four steps: (1) Identify K-mers where each K-mer consists K consecutive bases starting at each base position from the reference genome. (2) Align the K-mers back to the reference genome using a desired aligner, e.g. BWA (37). Ideally, the aligner and the alignment parameters are chosen to match what was used for generating read alignment files from the sample genomes. (3) Identify mappable base positions where the corresponding K-mers map back to themselves unambiguously (i.e. there is a single best hit and it is the true position of the K-mer). For example, the X0 field produced by BWA (37) relates a K-mer from a specific place in the genome to the number of best hits of that K-mer in the entire genome. If a K-mer has a X0 value of 1, the corresponding base can be identified as a mappable base. (4) Compute mappability score as the proportion of mappable bases in a given window, which measures the uniqueness of specific regions of the reference genome.
In summary, the input data is a triplet for each window represent by o t denotes the read-depth, g t denotes the GC content, and l t denotes the mappability score of the t th window.

HMM setup
We use a time-homogeneous discrete hidden Markov model (HMM) to segment the genome to regions of same copy number. In our HMM, time represents the sliding windows tiled along a chromosome, denoted by t.
The state represents the underlying copy number (CN). The state variable q t = CN t is hid-den and discrete with N possible values, (0, 1, ..., N − 1), where N , is derived from the data by K-mean clustering the logarithm of the read-depth. A particular sequence of the states is described by q = (q 1 , ..., q T ), where T is the total number of sliding windows of a chromosome. Let π j be the initial state probability, the probability that the state of the first window is state j. The underlying hidden Markov chain is defined by state transitions P (q t |q t−1 ) and is represented by a time-independent stochastic transition matrix A = {a jz } = P (q t = z|q t−1 = j).
Each copy number state emits an observation, the read-depth. The observation variable, O t , is a discrete count variable. A particular sequence of the observations is described by o = (o 1 , ..., o T ).
The emission probability of a particular observation at a particular time t for state j is described by e(t, j) = P (O t = o t |q t = j). For a detailed description of the emission probability, see Section

3.2.
We use the Baum-Welch algorithm (Baum et al 1970) to find the maximum likelihood estimates (MLE) of the HMM parameters. Following Bilmes (1998), we define the complete-date likelihood and solve the Q function in order to find the maximum likelihood estimates (MLE) of the HMM parameters.
A standard HMM assumes the Markov property, P (q t |q t−1 , q t−2 , q t−3 , ...q 1 ) = P (q t |q t−1 ). An additional assumption that is often employed is that the observations are independent given the states, O t ⊥ O i (i = t)|q t , which is valid when the windows are non-overlapping. When the windows are overlapping, this assumption is invalid; and instead, the observations are drawn from an autoregression process (Juang and Rabiner 1985). We have implemented an autoregressive HMM to model this feature of the data. For details of the autoregressive HMM, see Section 5.

Emission probability
The emission probability of the read-depth, e(t, j) = P (O t = o t |q t = j), is modeled as a mixture of a uniform distribution and a negative binomial distribution.
where c is the proportion of the random uniform component and is fixed as constant for each state; and R m is the largest read-depth among all windows and thus 1/R m is the uniform density.
To describe the negative binomially distributed component, e N B (t, j), we first explain the re-5 lationship between the Poisson and the negative binomial distributions. The Poisson distribution imposes that the variance equals to the mean. The negative binomial distribution allows overdispersion. Specifically, if O follows a Poisson distribution with mean µ, and µ follows a gamma distribution, the resulting distribution for O is a negative binomial distribution. The variance of negative binomial distribution is µ t + φµ 2 t , where φµ 2 t is the overdispersion part of the variance. As φ → 0, f N B (o t ; µ t , φ) reduces to a Poisson distribution with mean µ t and variance µ t .
Next, the mean value of the negative binomially distributed component is expressed as a function of a set of covariates to account for confounders.
( 2) where t denotes the t th window, j is the index of the copy number state, j emphasizes the dependency of the mean µ t on the copy number CN t , l t is the mappability score, g t is the GC content. For computational convenience, we set CN t = 0.5 when j = 0, and set CN t = j when j > 0.
We then employ a log link function to acknowledge the fact that µ tj > 0 and obtain: log(µ tj ) = β 0 + β 1 * log(CN t ) + β 2 * log(l t ) + β 3 * log(g t ) β 0 , β 1 , β 2 , β 3 are the regression coefficients. Specifically, β 0 = log(α 0 ), is the intercept parameter and is interpreted as the average level of read-depth signal when all covariates are equal to zero. β 1 is the amount of increase of read-depth for every unit increase of copy number, CN. β 2 is the amount of increase of read-depth for every unit increase of the mappability score, l. β 3 is the amount of increase of read-depth for every unit increase of the GC content, g.
Thus, given the above regression model for the mean, the negative binomial probability distribution function is expressed as the following: The complete emission probability is then expressed as the following: Because we fix c and R m as constant, parameter estimation will only concern the negative binomially distributed component.
To estimate the negative binomial parameters, a weighted GLM function is implemented in C++.
The argument "family" means either Poisson or negative binomial. The argument "prior" means the probability that each observation belongs to the negative binomially distributed component.
For the t th window and state j, the "prior" is denoted by p t,j and is computed as the following: Following the implementation function MASS/glm.nb in R (Venables and Ripley, 2002), we use an alternating iterative estimation procedure to obtain the new estimate of log(µ tj ), denoted by log(µ tj ), and the new estimate of φ, denoted by φ.
• First, we fix φ and compute log(µ tj ) by fitting weighted GLM using the iteratively reweighted least squares (IRLS) method. For details, see Section 4.3.4.
• Then, we fix log(µ tj ) and compute φ using the Newton-Raphson method with weight. For details, see Section 4.3.5.
• The above two steps alternated until convergence. If GC content (G) is the only covariate, M = 1 and define the covariate vector as: x t = {g t , ...g t , .., g t } (g t repeats for N times)

Perform a score test
The score test (Dean 1992) is used to test whether the overdispersion parameter, φ, is significantly greater than 0. If the score test is significant, we • Estimate φ using the Newton-Raphson method as described in Section 4.3.5.

Estimation of overdispersion using the Newton-Raphson method
Given log(µ tj ), we use the Newton-Raphson method to estimate the overdispersion parameter, φ.
In this study, we estimate one overdispersion parameter φ jointly for all states, and set φ j = φ for j = 0..N − 1.
The following weighted log-likelihood and its first, second derivatives are used in the Newton-Raphson method to estimate φ. The weighted log-likelihood is the same as Equation (28). where Ψ(x) = ∂ log Γ(x)/∂x, the digamma function. The observed Fisher information is where ψ(x) = ∂ 2 log Γ(x)/∂x 2 , the trigamma function.
We use the Newton-Raphson method given the score function and the fisher information.