Resolving Multicopy Duplications de novo Using Polyploid Phasing

Abstract

While the rise of single-molecule sequencing systems has enabled an unprecedented rise in the ability to assemble complex regions of the genome, long segmental duplications in the genome still remain a challenging frontier in assembly. Segmental duplications are at the same time both gene rich and prone to large structural rearrangements, making the resolution of their sequences important in medical and evolutionary studies. Duplicated sequences that are collapsed in mammalian de novo assemblies are rarely identical; after a sequence is duplicated, it begins to acquire paralog-specific variants. In this paper, we study the problem of resolving the variations in multicopy, long segmental duplications by developing and utilizing algorithms for polyploid phasing. We develop two algorithms: the first one is targeted at maximizing the likelihood of observing the reads given the underlying haplotypes using discrete matrix completion. The second algorithm is based on correlation clustering and exploits an assumption, which is often satisfied in these duplications, that each paralog has a sizable number of paralog-specific variants. We develop a detailed simulation methodology and demonstrate the superior performance of the proposed algorithms on an array of simulated datasets. We measure the likelihood score as well as reconstruction accuracy, i.e., what fraction of the reads are clustered correctly. In both the performance metrics, we find that our algorithms dominate existing algorithms on more than 93% of the datasets. While the discrete matrix completion performs better on likelihood score, the correlation-clustering algorithm performs better on reconstruction accuracy due to the stronger regularization inherent in the algorithm. We also show that our correlation-clustering algorithm can reconstruct on average 7.0 haplotypes in 10-copy duplication datasets whereas existing algorithms reconstruct less than one copy on average.

M.J. Chaisson and S. Mukherjee—Joint first authorship.

S. Kannan and E.E. Eichler—Joint last authorship.

A Appendix

After each gradient step, the resultant matrix is projected onto the box. The updates for A and B are as follows:

$$ \tilde{A}^{(t+1)} \leftarrow A^{(t)} - \alpha _{A} \nabla _{A} f(A) $$

Then $A_{ij}^{(t+1)} = {\left\{ \begin{array}{ll} 0, &{} \text {if } \tilde{A}^{(t+1)}_{ij} < 0 \\ \tilde{A}^{(t+1)}_{ij}, &{} \text {if } 0 \le \tilde{A}^{(t+1)}_{ij} \le 1 \\ 1, &{} \text {if } \tilde{A}^{(t+1)}_{ij} > 1 \end{array}\right. }$

$$ \tilde{B}^{(t+1)} \leftarrow B^{(t)} - \alpha _{B} \nabla _{A} f(B) $$

Then $B_{ij}^{(t+1)} = {\left\{ \begin{array}{ll} -1, &{} \text {if } \tilde{B}^{(t+1)}_{ij} < -1 \\ \tilde{B}^{(t+1)}_{ij}, &{} \text {if } -1 \le \tilde{A}^{(t+1)}_{ij} \le 1 \\ 1, &{} \text {if } \tilde{A}^{(t+1)}_{ij} > 1 \end{array}\right. }$

where $f(\cdot )$ is the objective function. The projected gradient descent allows us to incorporate additional constraints on the problem as well. If we further enforce that the sum of each row of A equals 1, then we would have the projection as $A_{ij}^{(t+1)} = \max \lbrace 0, \tilde{A}_{ij}^{(t+1)} - \nu _i \rbrace$ where $\nu _i$ can be computed for each row i using the equality

$$ \sum _{j=1}^S \max \lbrace 0, \tilde{A}_{ij}^{(t+1)} - \nu _i \rbrace =1 $$

We allow a maximum of 50 iteration steps for minimizing each of A and B, and 100 iteration steps for alternating minimization. We exit the iterations if the change in norm is insignificant ($1e-02$) or if the objective value change is below a tolerance ($1e-04$). The learning rate values have to be computed in order to ensure that gradient steps do not diverge. Our choices of learning rates have been

$$ \alpha _A = C \frac{\Vert \nabla f(A^{(t)})\Vert _F^2}{\Vert \mathcal {P}_\varOmega (\nabla f (A^{(t)}) \cdot B^{(t)}) \Vert _F^2} $$

and

$$ \alpha _B = C \frac{\Vert \nabla f(B^{(t)})\Vert _F^2}{\Vert \mathcal {P}_\varOmega ( A^{(t)} \cdot \nabla f (B^{(t)}) ) \Vert _F^2} $$

where $C \in (0,1)$.