Sparse Tensor Decomposition for Haplotype Assembly of Diploids and Polyploids

Background Haplotype assembly is the task of reconstructing haplotypes of an individual from a mixture of sequenced chromosome fragments. Haplotype information enables studies of the effects of genetic variations on an organism’s phenotype. Most of the mathematical formulations of haplotype assembly are known to be NP-hard and haplotype assembly becomes even more challenging as the sequencing technology advances and the length of the paired-end reads and inserts increases. Assembly of haplotypes polyploid organisms is considerably more difficult than in the case of diploids. Hence, scalable and accurate schemes with provable performance are desired for haplotype assembly of both diploid and polyploid organisms. Results We propose a framework that formulates haplotype assembly from sequencing data as a sparse tensor decomposition. We cast the problem as that of decomposing a tensor having special structural constraints and missing a large fraction of its entries into a product of two factors, U and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\underline {\mathbf {V}}$\end{document}V_; tensor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\underline {\mathbf {V}}$\end{document}V_ reveals haplotype information while U is a sparse matrix encoding the origin of erroneous sequencing reads. An algorithm, AltHap, which reconstructs haplotypes of either diploid or polyploid organisms by iteratively solving this decomposition problem is proposed. The performance and convergence properties of AltHap are theoretically analyzed and, in doing so, guarantees on the achievable minimum error correction scores and correct phasing rate are established. The developed framework is applicable to diploid, biallelic and polyallelic polyploid species. The code for AltHap is freely available from https://github.com/realabolfazl/AltHap. Conclusion AltHap was tested in a number of different scenarios and was shown to compare favorably to state-of-the-art methods in applications to haplotype assembly of diploids, and significantly outperforms existing techniques when applied to haplotype assembly of polyploids. Electronic supplementary material The online version of this article (10.1186/s12864-018-4551-y) contains supplementary material, which is available to authorized users.

under positive selection [4]. In addition to the SNPs and minor structural variations found in a healthy individual's genome, complex chromosomal aberrations such as translocations and nonreciprocal structural changesincluding aneuploidy -are present in cancer cells. Cancer haplotype assembly enables identification of "driver" mutations and thus helps to understanding the mechanisms behind the disease and discovery of its genetic signatures.
Haplotype assembly from short reads obtained by high-through-put DNA sequencing requires partitioning (either directly or indirectly) the reads into K clusters (K = 2 for diploids, K = 3 for triploids, etc.), each collecting the reads corresponding to one of the chromosomes. If the reads were free of sequencing errors, this task would be straightforward. However, sequencing is erroneous -state-of-the-art platforms have error rates on the order of 10 −3 − 10 −2 . This leads to ambiguities regarding the origin of a read and therefore renders haplotype assembly challenging. For this reason, the vast majority of haplotype assembly techniques attempts to remove the aforementioned ambiguities by either discarding or altering sequencing data; this has led to the minimum fragment removal, minimum SNP removal [5], maximum fragments cut [6], and minimum error correction formulations of the assembly problem [7]. Most of the recent haplotype assembly methods (see, e.g., [8][9][10][11][12]) focus on the minimum error correction (MEC) formulation where the goal is to find the smallest number of nucleotides in reads that need to be changed so that any read partitioning ambiguities would be resolved. It has been shown that finding optimal solution to the MEC formulation of the haplotype assembly problem is NPhard [5,12,13]. In [14], the authors used a branch-andbound scheme to minimize the MEC objective over the space of reads; to reduce the search space, they relied on a bound on the objective obtained by a random partition of the reads. Unfortunately, exponential growth of the complexity of this scheme makes it computationally infeasible even for moderate haplotype lengths. Integer linear programming techniques have been applied to haplotype assembly in [15], but the approach there fails at computationally difficult instances of the problem. More recently, fixed parameter tractable (FPT) algorithms with runtimes exponential in the number of variants per read [16,17] were proposed; these methods are well-suited for short reads but become infeasible for the long ones. A dynamic programming scheme for haplotype assembly of diploids proposed in [18] is also exponential in the length of the longest read. A probabilistic dynamic programming algorithm that optimizes a likelihood function generalizing the MEC objective is developed in [10]; this method is characterized by high accuracy but is significantly slower than the previous heuristics. Authors in [9,11] aim to process long reads by developing algorithms for the exact optimization of weighted variants of the MEC score that scale well with read length but are exponential in the sequencing coverage. These methods, along with ProbHap [10], struggle to remain accurate and practically feasible at high coverages (e.g., higher than 12 [10]).
The computational challenges of optimizing MEC score has motivated several polynomial time heuristics. In a pioneering work [19], a greedy algorithm seeking the most likely haplotypes was used to assemble haplotypes of the first complete diploid individual genome obtained via high-throughput sequencing. To compute posterior joint probabilities of consecutive SNPs, Bayesian methods relying on MCMC and Gibbs sampling schemes were proposed in [20] and [21], respectively; unfortunately, slow convergence of Markov chains that these schemes rely on limits their practical feasibility. Following an observation that haplotype assembly can be interpreted as a clustering problem, a max-cut formulation was proposed in [22]; an efficient algorithm (HapCUT, recently upgraded to Hap-CUT2 [23]) that solves it and significantly outperforms the method in [19] was developed and has widely been used. A flow-graph based approach in [24], HapCompass, re-examined fragment removal strategy and demonstrated superior performance over HapCUT. Other recent diploid haplotype assembly methods include a greedy max-cut approach in [25], convex optimization program for minimizing the MEC score in [26], a communicationtheoretic interpretation of the problem solved via belief propagation (BP) in [27], and methods that use external reference panels such as 1000 Genomes to improve accuracy of haplotype assembly in [28,29]. Note that deep sequencing coverage provided by state-of-the-art highthroughput sequencing platforms and the emergence of very long insert sizes in recent technologies (e.g., fosmid [25]) may enable assembly of extremely long haplotype blocks but also impose significant computational burden on the methods above.
Increased affordability, capability to provide deep coverage, and longer sequencing read lengths also enabled studies of genetic variations of polyploid organisms. However, haplotype assembly for polyploid genomes is considerably more challenging than that for diploids; to illustrate this, note that for a polyploid genome with k haplotype sequences of length m, under the all-heterozygous assumption there are (k − 1) m different genotypes and at least 2 (m−1) (k − 1) m different haplotype phasings. In part for this reason relatively fewer methods for solving the haplotype assembly problems in polyploids have been developed. In fact, with the exception of HapCompass [24], SDhaP [26] and BP [27], the above listed methods are restricted to diploid genomes. Other techniques capable of reconstructing haplotypes for both diploid and polyploid genomes include HapTree [30], a Bayesian method to find the maximum likelihood haplotype shown to be superior to HapCompass and SDhaP (see, e.g., [31] for a detailed comparison), H-PoP [8], the state-ofthe-art dynamic programming method that significantly outperforms the schemes developed in [24,26,30] in terms of accuracy, memory consumption, and speed, and the recently proposed matrix factorization schemes in [32,33].
In this paper, we propose a unified framework for haplotype assembly of diploid and polyploid genomes based on sparse tensor decomposition; the framework essentially solves a relaxed version of the NP-hard MEC formulation of the haplotype assembly problem. In particular, read fragments are organized in a sparse binary tensor which can be thought of as being obtained by multiplying a matrix that contains information about the origin of erroneous sequencing reads and a tensor that contains haplotype information of an organism. The problem then is recast as that of decomposing a tensor having special structural constraints and missing a large fraction of its entries. Based on a modified gradient descent method and after unfolding the observed and haplotype information bearing tensors, an iterative procedure for finding the decomposition is proposed. The algorithm exploits underlying structural properties of the factors to perform decomposition at a low computational cost. In addition, we analyze the performance and convergence properties of the proposed algorithm and determine bounds on the minimum error correction (MEC) scores and correct phasing rate (CPR) -also referred to as reconstruction rate -that the algorithm achieves for a given sequencing coverage and data error rate. To the best of our knowledge, this is the first polynomial time approximation algorithm for haplotype assembly of diploids and polyploids having explicit theoretical guarantees for its achievable MEC score and CPR. The proposed algorithm, referred to as AltHap, is tested in applications to haplotype assembly for both diploid and polyploid genomes (synthetic and real data) and compared with several stateof-the-art methods. Our extensive experiments reveal that AltHap outperforms the competing techniques in terms of accuracy, running time, or both. It should be noted that while state-of-the-art haplotype assembly methods for polyploids assume haplotypes may only have biallelic sites, AltHap is capable of reconstructing polyallelic haplotypes which are common in many plants and some animals, are of particular importance for applications such as crop cultivation [34], and may help in reconstruction of viral quasispecies [35]. Moreover, while many state-of-the-art haplotype assembly methods are computationally intensive (e.g., [10,15]), our extensive numerical experiments demonstrate efficacy of AltHap in a variety of practical settings.

Problem formulation
We briefly summarize notation used in the paper. Bold capital letters refer to matrices and bold lowercase letters represent vectors. Tensors are denoted by underlined bold capital letters, e.g., M. M ::1 and M denote the frontal slice and the mode-1 unfolding of a third-order tensor M, respectively. For a positive integer n, [ n] denotes the set {1 . . . , n}. The condition number of rank-k matrix M is defined as κ = σ 1 /σ k where σ 1 ≥ · · · ≥ σ k > 0 are singular values of M. SVD k (M) denotes the rank k approximation (compact SVD) of M computed by power iteration method [36,37].
Let H = {h 1 , . . . , h k } denote the set of haplotype sequences of a k-ploid organism, and let R be an n×m SNP fragment matrix where n denotes the number of sequencing reads and m is the length of haplotype sequences. R is an incomplete matrix that can be thought of as being obtained by sampling, with errors, matrix M that consists of n rows; each row of M is a sequence randomly selected from among k haplotype sequences. Since each SNP is one of four possible nucleotides, we use the alphabet A = {(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)} to describe the information in the haplotype sequences; the mapping between nucleotides and alphabet components follows arbitrary convention. The reads can now be organized into an n × m × 4 SNP fragment tensor which we denote by R. The (i, j, :) fiber of R, i.e., a one-dimensional slice obtained by fixing the first and second indices of the tensor, represents the value of the j th SNP in the i th read. Let denote the set of informative fibers of R, i.e., the set of (i, j, :) such that the i th read covers the j th SNP. Define an operator P (.) as P (R) is a tensor obtained by sampling, with errors, tensor M ∈ A n×m having n copies of k encoded haplotype sequences as its horizontal slices. More specifically, we can write M = UV , where V ∈ A m×k contains haplotype information, i.e., the j th vertical slice of V, V :j: , is the encoded sequence of the j th haplotype, and U ∈ {0, 1} n×k is a matrix that assigns each of n horizontal slices of M to one of k haplotype sequences, i.e., the i th row of U, u i , is an indicator of the origin of the i th read. Let = {e 1 , . . . , e k }, where e l ∈ R k is the l th standard basis vector having 1 in the l th position and 0 elsewhere. The rows of U are standard unit basis vectors in R k , i.e., u i ∈ , ∀i ∈[ n]. This representation is illustrated in Fig. 1 where the (1, 1, :) fiber of V specified with dashed lines is mapped to the (1, 1, :) fiber of M which in turn implies that in the example described in Fig. 1 we have u 1 = e 1 . Fig. 1 Representing haplotype sequences and sequencing reads using tensors. Tensor V ∈ A m×k contains haplotype information while matrix U ∈ {0, 1} n×k assigns each of the n horizontal slices of M to one of the k haplotype sequences, i.e., the i th row of U is an indicator of the origin of the i th read DNA sequencing is erroneous and hence we assume a model where the informative fibers in R are perturbed versions of the corresponding fibers in M with data error rate p e , i.e., if the (i, j, :) ∈ fiber in M takes value e l ∈ A, R ij: with probability 1 − p e equals e l and with probability p e takes one of the other three possibilities. Thus, the observed SNP fragment tensor can be modeled as R = P (M + N) where N is an additive noise tensor defined as where the notation U (A\{M ij: }) denotes uniform selection of a vector from A\{M ij: }. The goal of haplotype assembly can now be formulated as follows: Given the SNP fragment tensor R, find the tensor of haplotype sequences V that minimizes the MEC score. Next, we formalize the MEC score as well as the correct phasing rate, also known as reconstruction rate, the two metrics that are used to characterize performance of haplotype assembly schemes (see, e.g., [15,18,38,39]). For two alleles a 1 , a 2 ∈ A ∪ {0}, we define a dissimilarity function d(a 1 , a 2 ) as The MEC score is the smallest number of fibers in R that need to be altered so that the resulting modified data is consistent with the reconstructed haplotype V, i.e., The correct phasing rate (CPR), also referred to as the reconstruction rate, can conveniently be written using the dissimilarity function d(., .). Let V t denote the tensor of true haplotype sequences. Then where M is a one-to-one mapping from lateral slices of V to those of V t , i.e., a one-to-one mapping from the set of reconstructed haplotypes to the set of true haplotypes. We now describe our proposed relaxation of the MEC formulation of the haplotype assembly problem. It is straightforward to observe that the last equality in (6) can equivalently be written as where V is the mode-1 unfolding of the tensor V. Hence, Let U ∈ {0, 1} n×k be the matrix such that for its i th row it holds that u i = e p i . In addition, notice that vec(R i:: ) is the i th row of R. Therefore, from the definition of the Frobenius norm and the fact that P (R) = R we obtain The optimization problem in (7) is NP-hard since the entries of V are binary and the objective function is nonconvex. Relaxing the binary constraint to 1], results in the following relaxation of the MEC formulation, The new formulation can be summarized as follows. We start by finding the so-called mode-1 unfolding of tensors M and V and denote the decomposition M = UV , as illustrated in Fig. 2. As implied by the figure, after unfolding, the entries of the (1, 1, :) fiber are mapped to four blocks of M and V that correspond to the frontal slices of tensors M and V, respectively. Then, to determine the haplotype sequence that minimizes the MEC score, one needs to solve (8) and find the optimal tensor decomposition.

The AltHap algorithm
Although the objective function in (8), i.e., is convex in each of the factors when the other factor is fixed, f (U, V) is generally nonconvex. To facilitate computationally efficient search for the solution of (8), we rely on a modified gradient search algorithm which exploits the special structures of U and V and iteratively updates the More specifically, given the current estimates (U t , V t ), the update rules are where α is a judiciously chosen step size, and C denotes the projection operator onto C. Notice that the optimization in (9) is done by exhaustively searching over k vectors in . Since the number of haplotypes k is relatively small, the complexity of the exhaustive search (9) is low. The proposed scheme is formalized as Algorithm 1.

Algorithm 1 Structured Tensor Decomposition
Input: SNP fragment matrix R, step size α, maximum number of iterations T Output: V, an estimate of the true haplotype tensor V t Preprocessing: Encode R to tensor R and find the mode-1 unfolding, R Initialization: Using power method, Compute Note that AltHap differs from a previously proposed SCGD algorithm in [32] as follows: (i) AltHap's novel representation of haplotypes and sequencing reads using binary tensors provides a unified framework for haplotype assembly of diploids as well as biallelic and polyallelic polyploids. The method in [32] is not capable of performing haplotype assembly of polyallelic polyploid Fig. 2 Representing haplotype sequences and sequencing reads using unfolded tensors. Matrix V ∈ {0, 1} 4m×k contains haplotype information while matrix U ∈ {0, 1} n×k assigns each of the n rows of M to one of the k haplotype sequences, i.e., the i th row of U is an indicator of the origin of the i th read genomes. (ii) Unlike [32], AltHap exploits the fact that V is composed of binary entries by imposing the constraint V i,j ∈ C in the MEC relaxation in (8). As our results in Section 5 demonstrate, this leads to significant performance improvements of AltHap over SCGD in a variety of settings. (iii) Lastly, in Section 4 we provide analysis of the global convergence of AltHap and derive explicit analytical bounds on its achievable performance. Such performance guarantees do not exist for the method in [32].

Convergence analysis of AltHap
In this section, we analyze the convergence properties of AltHap and provide performance guarantees in different scenarios.
In the Additional file 1 we show that, a judicious choice of the step size α according to where C ∈ (0, 2) is a constant, guarantees that the value of the objective function in (8) decreases as one alternates between (9) and (10), which in turn implies that AltHap converges. The key observation that leads to this result is that f (U, V) is a convex function in each of the factor matrices and that C =[ 0, 1] is a convex set; hence the projection C in (10) leads to a reduction of f U t , V t in each iteration t. It is important however to determine the conditions under which the stationary point of AltHap coincides with the global optima of (8). To this end, we first provide the definition of incoherence of matrices [40].
Let each fiber in MT be observed uniformly with probably p. Let C snp = mp denote the expected number of SNPs covered by each read, and C seq = np denote the expected coverage for each of the haplotype sequences. Theorem 1 built upon the results of [41][42][43] states that with an adequate number of covered SNPs, the solution found by AltHap reconstructs M up to an error term that stems from the existence of errors in sequencing reads.
Suppose the condition number of M is κ. Then there exist numerical constants C 0 , C 1 > 0 such that if is uniformly generated at random and C snp > max C 0 3 μ 4 k 14 κ 12 C seq , p e k 2 κ 6 2C 1 (13) with probability at least 1− 1 m 3 , the solution (U * , V * ) found by AltHap satisfies The proof of Theorem 1 relies on a coupled perturbation analysis to establish a certain type of local convexity of the objective function around the global optima. Thus, under (13) there is no other stationary point around the global optima and hence, starting from a good initial point, AltHap converges globally. We employ the initialization procedure suggested by [42] -summarized in the initialization step of Algorithm 1 -which is based on a low cost singular value decomposition of R using power method [36,37] and with high probability lies in the described convexity region of f (U, V).  (14) is the bound on the error of the solution generated by AltHap which increases with the sequencing error rate p e and ploidy k, and decreases with C snp and the number of reads n, as expected.

Remark 2
If M is well-conditioned, i.e., M is characterized by a small incoherence parameter μ and a small condition number κ, the recovery becomes easier; this is reflected in less strict sufficient condition (13) and improved achievable performance (14). In fact, as we verified in our simulation studies, by using the proposed framework for haplotype assembly, the parameters μ and κ associated with M are close to 1 (the ideal case). Theorem ( 15) suggests that the expected MEC increases with the length of the haplotype sequences, sequencing error, number of haplotype sequences, and sequencing coverage. A higher sequencing coverage results in a larger fragment data which in turn leads to higher MEC scores.

Remark 4 As intuitively expected, the bound ( 16) suggests that AltHap's achievable expected CPR improves with the number of sequencing reads and the SNP coverage; on
the other hand, the CPR deteriorates at higher data error rates. Finally, assuming the same sequencing parameters, ( 16) implies that reconstruction of polyploid haplotypes is more challenging than that of diploids.

Results and discussion
We evaluated the performance of the proposed method on both experimental and simulated data, as described next. AltHap was implemented in Python and MATLAB, and the simulations were conducted on a single core Intel Xeon E5-2690 v3 (Haswell) with 2.6 GHz and 64 GB DDR4-2133 RAM. The benchmarking algorithms include Belief Propagation (BP) [27], a communication-inspired method capable of performing haplotype assembly of diploid and biallelic polyploid species, HapTree [30], integer linear programming (ILP) technique [15], SCGD [32], and H-PoP [8], the state-of-the-art dynamic programming algorithm for haplotype assembly of diploid and biallelic polyploid species shown to be superior to HapTree [30], HapCompass [24], and SDhaP [26] in terms of both accuracy and speed [8,31]. Following the prior works on haplotype assembly (see, e.g., [15,18,38,39]) we use MEC score and CPR to assess the quality of the reconstructed haplotypes. For clarity, in the tables we report the CPR percentage, i.e., CPR× 100.

Experimental data
We first tested performance of AltHap in an application to haplotype reconstruction of a data set from the 1000 Genomes Project -in particular, the sample NA12878 sequenced at high coverage using the 454 sequencing platform. In this work, we take the trio-phased variant calls from the GATK resource bundle [44] as the true haplotype sequences. We compare the MEC score, CPR, and running time achieved by AltHap to those of H-PoP, BP, HapTree, SCGD and ILP. All the algorithms used in the benchmarking study were executed with their default settings. The results are given in Table 1. As seen there, among the considered algorithms AltHap achieves the highest CPR for majority of the chromosomes (nine), followed by H-PoP and ILP (five each). As expected, ILP achieves the lowest MEC scores among all the methods but this comes at a computational cost much higher than that of AltHap. Notice that lower MEC score does not necessarily imply better CPR. MEC is the error evaluated on observed SNPs positions, i.e., the training data points, while CPR is related to the generalization error that is calculated on unobserved SNPs positions, i.e., the test data points. Since the sequencing reads are erroneous, an algorithm might over-fit while trying to minimize the MEC score.
Fosmid pool-based sequencing provides very long fragments and is characterized by much higher ratio of  the number of SNPs to the number of reads than the standard data sets generated by high-throughput sequencing platforms. We consider the fosmid sequence data for chromosomes of HapMap NA12878 and again take the trio-phased variant calls from the GATK resource bundle [44] as the true haplotype sequences. We compare the performance of AltHap to those of H-PoP, BP, HapTree, SCGD and ILP and report the results in Table 2. As can be seen from Table 2, AltHap achieves the best CPR for most of the chromosomes (thirteen out of 22) followed by  The best results in each simulation setting are in bolface font H-PoP (four). As with the 1000 Genome Project Data, ILP achieves the best MEC scores but is much slower and significantly inferior to AltHap in terms of CPR. Note that since HapTree could not finish assembling haplotype of the 6 th chromosome in 48 hours, that result is missing from the table.

Simulated data: the diploid case
To further benchmark performance of the proposed scheme, we test it on the synthetic data from [39] often used to compare methods for haplotype assembly of diploids. These data sets emulate haplotype assembly under varied coverage, sequencing error rates and haplotype block lengths. We constrain our study to the assembly of haplotype blocks having length m = 700 bp (the longest blocks in the data set). The results, averaged over 100 instances of the problem, are given in Table 3. As evident from this table, AltHap outperforms other algorithms for nearly all the combinations of data error rates and sequencing coverage and is also much faster than SCGD, ILP, BP and HapTree while being slightly slower than H-PoP. Note that ILP could only finish assembling haplotype of two settings with p e = 0.1 and coverages of 5 and 8, in 48 hours. Hence, the results for other settings are missing from the table.

Simulated data: the biallelic polyploid case
The performance of AltHap in applications to haplotype assembly for polyploids was tested using simulations; in particular, we studied how AltHap's accuracy depends on coverage and sequencing error rate. The generated data sets consist of paired-end reads with long inserts that emulate the scenario where long connected haplotype blocks need to be assembled. We simulate sampling of the entire genome using paired-end reads and generate SNPs along the genome with probability 1 in 300. In other words, the distance between pairs of adjacent SNPs follows a geometric random variable with parameter p snp = 1 300 (the SNP rate). To emulate a sequencing process capable of facilitating reconstruction of long The best results in each simulation setting are in bolface font haplotype blocks, we randomly generate paired-end reads of length 2 × 250 with average insert length of 10,000 bp and the standard deviation of 10%; sequencing errors are inserted using realistic error profiles [45] and genotyping is performed using a Bayesian approach [44]. At such read and insert lengths, the generated haplotype blocks are nearly fully connected. Each experiment is repeated 10 times. AltHap is compared with H-PoP, BP and SCGD.
We also tried to run HapTree. However, HapTree could not finish the simulations for the considered block size in 48 hours. Table 4 compares the CPR, MEC score, and running times of AltHap with those of H-PoP, BP and SCGD for biallelic triploid genomes with haplotype block lengths of m = 1000 for several combinations of sequencing coverage and data error rates. As can be seen there, in terms of CPR AltHap outperforms all other meth- ods in all the scenarios while in terms of MEC score it outperforms other methods in the vast majority of the scenarios. Note that unlike H-PoP's, the complexity of AltHap scales gracefully with coverage (i.e., although H-PoP is very fast at low coverages, at the highest coverage AltHap becomes much faster than H-PoP). As can be seen in Table 6, overall CPR score (MEC score) of all algorithms decreases (increases) as the probability of error increases. This is expected -and also reflected in our theoretical results -since with higher data error rate haplotype assembly becomes more challenging. Additionally, MEC scores increases with coverage since higher coverage implies more sequencing reads. Therefore, the total number of observed SNP positions increases which in turn results in higher MEC scores.
The results of tests conducted on simulated biallelic tetraploid genomes are summarized in Table 5, where we observe that AltHap outperforms the competing schemes in terms of both accuracy and running time. To investigate how the performance and complexity of AltHap vary with coverage and read length, in Table 6 we report its CPR, MEC, and runtimes obtained by simulating assembly of biallelic triploid haplotypes using paired end reads of length 2 × 250, 2 × 300, and 2 × 500 and coverage 10, 20 and 30 (block length is set to m = 1000 and data error rate is p e = 0.002). The results imply that AltHap's runtime scales approximately linearly with both read length and coverage over the consider range of these two parameters. Additionally, as we see in Table 5, MEC score slightly increases with read length. The impact of read length in this matter is similar to that of sequencing coverage. longer sequencing reads provide more observed SNP positions and hence the MEC might increase, as also predicted by our theoretical results.

Simulated data: the polyallelic polyploid case
We further studied the performance of AltHap on triploid and tetraploid organisms having polyallelic sites and the results are summarized in Tables 7 and 8, respectively. Notice that none of the competing schemes are capable of handling polyallelic genomes. Evidently, AltHap is able to reconstruct underlying haplotype sequences with competitive performance at a low computational cost. The results of these extensive simulations imply that, as expected, haplotype assembly becomes more challenging as the number of haplotype sequences (i.e., the ploidy) increases. Nevertheless, in all the conducted studies, AltHap consistently reconstructs haplotype sequences accurately and with practically feasible computational cost. In addition, the results of Tables 4 and 5 demonstrate that the computational time of AltHap grows significantly slower with coverage than the computational time of the competing schemes. In particular, for high coverages that are characteristic of high-throughput sequencing technologies, AltHap is the most efficient among the considered algorithm.

CPR lower bound
Finally, we use the results obtained by running AltHap on simulated biallelic triploid data (i.e., the results summarized in Table 4) to examine tightness of the theoretical bounds on the CPR stated in Theorem 2. In particular, theoretical bounds on CPR are compared to the CPRs empirically computed for various combinations of coverage and data error rates (averaged over 10 independent problem instances). In Fig. 3, the theoretical bound and experimental CPR results are shown as functions of the data error rate for coverage 15. We observe that the bound is reasonably close to the experimental results over the considered range of data error rates. In Fig. 4, the theoretical bound and experimental CPR results are plotted against sequencing coverage for the data error rate p e = 0.002. This figure, too, implies that the theoretical CPR bound is relatively close to the experimental results. Notice that as Fig. 3 shows, the CPR score as well as the lower bound derived in Theorem 2 decrease when the sequencing error increases. On the other hand, Fig. 4 depicts that higher coverage improves AltHap's CPR score, which again is reflected in our theoretical results.

Conclusion
In this paper, we developed a novel haplotype assembly framework for both diploid and polyploid organisms that relies on sparse tensor decomposition. The proposed algorithm, referred to as AltHap, exploits structural properties of the problem to efficiently find tensor factors and thus assemble haplotypes in an iterative fashion, alternating between two computationally tractable optimization tasks. If the algorithm starts the iterations from an appropriately selected initial point, AltHap converges to a stationary point which is with high probability in close proximity of the solution that is optimal in the MEC sense. In addition, we analyzed the performance and convergence properties of AltHap and found bounds on its achievable MEC score and the correct phasing rate. AltHap, unlike the majority of existing methods for haplotype assembly for polyploids, is capable of reconstructing haplotypes with polyallelic sites, making it useful in a number of applications involving plant genomes. Our extensive tests on real and simulated data demonstrate that AltHap compares favorably to competing methods in applications to haplotype assembly of diploids, and significantly outperforms existing techniques when applied to haplotype assembly of polyploids. As part of the future work, it is of interest to extend the sparse tensor decomposition framework to viral quasispecies reconstruction and recovery of bacterial haplotypes from metagenomic data.