{'fraction_non_alphanumeric': 0.03407604580187705, 'fraction_numerical': 0.00855916840680144, 'mean_word_length': 4.618151347170298, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Iterative Learning for Reference-Guided DNA Sequence Assembly from Short Reads: Algorithms and Limits of Performance

Recent emergence of next-generation DNA sequencing technology has enabled acquisition of genetic information at unprecedented scales. In order to determine the genetic blueprint of an organism, sequencing platforms typically employ so-called shotgun sequencing strategy to oversample the target genome with a library of relatively short overlapping reads. The order of nucleotides in the reads is determined by processing the acquired noisy signals generated by the sequencing instrument. Assembly of a genome from potentially erroneous short reads is a computationally daunting task even in the scenario where a reference genome exists. Errors and gaps in the reference, and perfect repeat regions in the target, further render the assembly challenging and cause inaccuracies. In this paper, we formulate the reference-guided sequence assembly problem as the inference of the genome sequence on a bipartite graph and solve it using a message-passing algorithm. The proposed algorithm can be interpreted as the well-known classical belief propagation scheme under a certain prior. Unlike existing state-of-the-art methods, the proposed algorithm combines the information provided by the reads without needing to know reliability of the short reads (so-called quality scores). Relation of the message-passing algorithm to a provably convergent power iteration scheme is discussed. To evaluate and benchmark the performance of the proposed technique, we find an analytical expression for the probability of error of a genie-aided maximum a posteriori (MAP) decision scheme. Results on both simulated and experimental data demonstrate that the proposed message-passing algorithm outperforms commonly used state-of-the-art tools, and it nearly achieves the performance of the aforementioned MAP decision scheme.


I. INTRODUCTION
In the last decade, rapid development of next-generation DNA sequencing technologies has enabled cheap and fast generation of massive amounts of sequencing data [1], [2], [3]. Determining the order of nucleotides in a long target DNA molecule typically involves use of shotgun sequencing strategy where multiple copies of the target are fragmented into short templates. Each template is then analyzed by a sequencing instrument which provides reads (i.e., information about nucleotide content of the templates) that are used to assemble the desired long sequence. Shotgun sequencing strategy is illustrated in Fig. 1. The majority of next-generation sequencing methods detect the order of nucleotides in a template by facilitating enzymatic synthesis of a complementary strand on the template, and acquiring a signal that indicates the type of a nucleotide successfully incorporated into the complementary strand. Quality of the acquired signal is adversely affected by various sources of uncertainty. Base calling algorithms attempt to infer the order of nucleotides in short templates from the acquired noisy signals [4]. Confidence in the accuracy of base calls is expressed by their quality scores, which provide a measure of the probability of base calling error. Conventional base calling algorithms rely on various heuristics to estimate quality scores, while more recent methods employ Bayesian inference schemes to evaluate posteriori probabilities of the bases in the reads [5], [6], [7]. Due to their computational efficiency, the heuristic methods are the preferred choice for assessing confidence of the base calls -however, the resulting quality scores are not necessarily an accurate reflection of the actual probability of base calling errors [8]. In turn, inaccurate quality scores may adversely impact reliability of the assembly process.
The short reads generated by a sequencing instrument are used to assemble the target genome. The assembly may be performed with or without referring to a previously determined sequence related to the target (genome, transcriptome, proteins). De novo assembly refers to a scenario where the reconstruction is performed without a reference sequence. This is a computationally challenging task, difficult due to the presence of perfect repeat regions in the target sequence and short lengths of the reads [11], [12]. In re-sequencing projects where the goal, for instance, may be to study genetic variations among individuals or to discover new strains of bacteria [9], [10], a reference is available and used to order the reads. Such reference-guided assembly is still challenging due to the errors in the reads and because the reference often contains errors and gaps [13], [14]. Many of the assembly challenges are ameliorated if the target sequence is significantly oversampled and thus the information provided by short reads is highly redundant. This redundancy is quantified by means of a sequencing coverage -the average number of times a base in the target sequence is present in the overlapping reads. However, the demands for higher throughput and lower sequencing costs often limit the coverage to medium (5-20X) or low (≤ 5X). As an example, the ongoing 1000 Genomes Project has opted for trading-off sequencing depth for the number of individuals being sequenced [18]. In its preliminary phase, the project has focused on sequencing a large number of individuals at a very low 3X coverage.
In the reference-guided assembly, the short reads are first mapped to a reference sequence using an alignment algorithm (e.g., [16], [17]). Then each position along the target is determined by combining information provided by all the reads that cover that particular position. Due to the errors in base calls, short length of the reads, and repetitiveness in the target, both the mapping and the sequence assembly steps are potentially erroneous. The widely used tools to analyze and assemble genome sequence from high-throughput sequencing data include SAMtools [13] and Genome Analysis Toolkit (GATK) [14]. Note that both of these packages rely on the quality scores provided by the sequencing platform to infer the assembled sequence.
In this paper, we formulate the reference-guided assembly problem as the inference of the genome sequence on a bipartite graph and solve it using a message-passing algorithm. Unlike existing state-ofthe-art methods, the proposed algorithm seeks the target sequence without needing to know reliability of the short reads (i.e., their quality scores). Instead, it infers reliability of a base in the assembled sequence by combining the information of all the reads covering that particular position. The proposed algorithm can be interpreted as the classical belief propagation under a certain prior. Binary reformulation of the problem leads to an alternative solution in the form of another message passing algorithm that is closely related to the so-called power iteration method. The power iteration algorithm approximates the solution to the sequence assembly problem by the leading singular vector of a matrix comprising read data. The power iteration method has guaranteed convergence, and its careful examination provides relation between the algorithm accuracy and the number of iterations. To evaluate and benchmark performance of the proposed techniques, we find an analytical expression for the probability of error of a genieaided maximum a posteriori (MAP) sequence assembly scheme which is an idealized assembler with perfect quality score information and error-free mapping of the reads to their locations. Results on both simulated and experimental data obtained by sequencing Escherichia Coli and Neisseria Meningitidis at UT Austin's Center for Genomic Sequencing and Analysis demonstrate that our proposed messagepassing algorithm performs close to the aforementioned genie-aided MAP assembly scheme and is superior compared to state-of-the-art methods (in particular, it outperforms the aforementioned SAMtools and GATK software packages). Note that the developed algorithms as well as simulation and experimental studies are focused on haploid genomes -while modifications that enable application to diploid/polyploid genomes are relatively straightforward, they are beyond the scope of the current manuscript.
The paper is organized as follows. In Section II, we introduce the bipartite graphical model and a message passing based sequence assembly algorithm. In Section III, we show that this message passing algorithm can be interpreted as the classical belief propagation under a certain prior. In Section IV, we derive an alternative message passing scheme based on a binary reformulation of the sequence assembly problem. In Section V, we derive an expression for the probability of error of a genie-aided MAP assembly scheme. In Section VI, we present simulations as well as experimental results obtained by applying the proposed methods on the Escherichia Coli and Neisseria Meningitidis data sets. Section VII concludes the paper and outlines potential future work.
Preliminary work on the basic message passing scheme discussed in Section III was presented in [19], [20]. Implementation code of the algorithm in C++ is available at https://sourceforge.net/projects/mpsequencing/.

II. GRAPHICAL MODEL AND THE MESSAGE-PASSING ASSEMBLY ALGORITHM
To facilitate processing of the short reads generated by next-generation sequencing instruments, we introduce a bipartite graph representing the reads and bases in the target sequence that needs to be assembled. The fundamental building blocks of a sequence -the nucleotides A, C, G, and T -are numerically represented using 4-dimensional unit vectors containing a single non-zero component whose position indicates type of a nucleotide. In particular, the 4-dimensional unit vectors that we use are For convenience, we will assume that all the short reads at our disposal are generated by the same sequencing instrument and thus have identical read length l; note, however, that there is no loss of generality and that our scheme can combine reads generated by sequencing the same target on different instruments and of different read lengths. Let us denote the set of short reads by R = {r j }, j = 1, 2, . . . , n. In general, the base calls in these reads are erroneous due to various uncertainties in the underlying sequencing-by-synthesis process. Average base-calling error rates of most current next-generation sequencing systems are on the order of 10 −2 .
In applications where a reference sequence is available (i.e., in the so-called reference guided sequence assembly scenario), the short reads are mapped onto the reference using one among many recently developed short-sequence alignment algorithms [16], [17]. Note that the reads comprising bases with low quality scores are often discarded by the alignment algorithms. Ideally, the remaining reads (the ones of high fidelity) are accurately mapped to their corresponding locations on the reference sequence. However, for some reads there may exist several candidate positions which leads to possible mis-alignments and, consequently, might provide erroneous information about the regions of the target sequence where the read is mis-aligned. As pointed out in [14], the misalignment rate for reads from genome regions which contains homozygous indels (where the chromosomes in a homologous pair have the same sequence but contain insertion or deletion as compared to the reference genome) can be as high as 15%.
We represent the reference-guided assembly problem by a graph shown in Fig. 2. The bipartite graph G(b 1:L ∪ r 1:n , E) illustrated in the figure has L base nodes (representing the target genome sequence) and n read nodes. Since we assume that all reads are of the same length l, each read node is connected to exactly l base nodes. The edge (i, j) in the edge set E connecting b i and r j is associated with a unit Fig. 2. Illustration of the reference-guided DNA sequence assembly problem using short reads. Nodes bi represent bases in the target DNA sequence and rj represent reads. Each read node is connected to l base nodes, where l denotes the read length.
vector e ij indicating information about the type of base b i provided by the read r j . Note that the bipartite graph described here is reminiscent of the graphical representation of the crowdsourcing problem in [21]. Motivated by the iterative learning scheme proposed there, we employ a message-passing algorithm to infer the target genome sequence using overlapping reads. Note that the previously mentioned problem of having multiple candidate locations for mapping reads can be incorporated in the proposed graphical representation and resolved using the algorithm that we describe next.
The message passing algorithms rely on the exchange of messages between neighboring nodes in the graph [24]. Our algorithm operates on real-valued base messages {x i→j } (i,j)∈E and read messages {y j→i } (i,j)∈E . A base message x i→j is a 4 × 1 vector representing the likelihood of the base b i being A, C, G, or T, while a read message y j→i represents the reliability of read j. Read messages are initialized from a random distribution, and the message update rules at iteration k are given by where ∂i and ∂j denote collection of the neighboring nodes of nodes i and j, respectively, and 1 is a 4 × 1 vector containing all 1's. Note that 2e ij − 1 has element 1 in the position corresponding to the nucleotide base b ij represented by e ij , and −1's elsewhere. Hence a read with positive reliability value y j →i will increase the likelihood of b ij and decrease the likelihood of other bases. Finally, the likelihood of a base being A, C, G, or T is calculated as the sum of the information provided by the reads weighted by each read's reliability. The symbol with the highest likelihood is chosen as the estimate of the base in the corresponding position. The estimate rule for the i th base iŝ where the decision vector j→i . Here, k m denotes the number of iterations performed and x {t} i denotes the likelihood corresponding to symbol t ∈ {A, C, G, T } in the vector The procedure is formalized as Algorithm 1. Note that Algorithm 1 needs to be appropriately initialized. In our experimental studies presented in Section VI, we initialize y j→i by drawing from both Gaussian distribution N (1, 1) and uniform distribution U [0, 1]. For the data sets under consideration, it turns out that different initializations lead to identical solutions. The algorithm is terminated when the reliability increment between subsequent iterations is small, i.e., |y j→i | < . As pointed out earlier, the algorithm does not require exact knowledge of quality scores, and iteratively infers reliability of individual reads.
Since the reads originating from a single sequencing instrument have identical lengths, the degree of the read nodes in the graph is uniform. On the other hand, degree of a base node is the number of reads that Algorithm 1 Message passing for sequence assembly cover the corresponding base, usually referred to as the sequencing coverage. Typically, coverage varies from one position to another and, consequently, degree of the base nodes varies. Note that fragmentation of multiple copies of the target sequence -a fundamental step in shotgun sequencing procedure -can be viewed as a uniform sampling from the original DNA strand. The resulting coverage is a random variable that can be described well by a Poisson distribution [26], [28].
Letc denote the average sequencing coverage. The computational complexity of the base message updating step (1), which needs to be performed in each iteration of Algorithm 1, is O(nlc) on average, while the complexity of the read message updating step (2) is O(nl 2 ). Since Lc = nl, the complexity of the algorithm is O(k m nl(l+c)) = O(k m Lc(l+c)), where k m denotes the number of iterations (i.e., the number of message updates). On the other hand, simple plurality voting scheme has complexity O(Lc logc). Our experimental studies show that k m ≤ 30 is sufficient for the convergence of the algorithm. We tested the algorithm on a broad range of parameters (in particular, for read lengths l ≤ 100, coveragec ≤ 60), and found that the runtimes are comparable to those of the state-of-the-art techniques (SAMtools and GATK) -a specific comparison of runtimes is reported in Section VI.

III. RELATION TO STANDARD BELIEF PROPAGATION
As an alternative to the intuitively pleasing but basically heuristic message passing scheme proposed in Section II, we can also derive a standard belief propagation algorithm for the reference-guided sequence assembly. To this end, we seek the sequenceb 1:L that maximizes the joint probability P (b 1:L , p 1:n ), where p 1:n denotes confidences of the aligned read data and p j ∈ [0, 1]. This maximization can be formalized as where D(p i ) denotes the prior distribution on p i andp j = 1 − p j . δ(·) denotes an indicator function taking value 1 if its argument is true and is 0 otherwise. The joint optimization is computationally challenging and thus often practically not feasible. As an alternative, belief propagation provides an approximate solution to (4) by computing the marginal distributions of the optimization variables and selecting their most likely values according to the computed distributions. A thorough review of theoretical and practical aspects of the belief propagation method can be found in [23]. For the graphical model proposed in Section II, we define two messages to facilitate belief propagation:x i→j andỹ j→i . The former is the belief onb i and essentially represents a distribution over the four possible nucleotide bases {e A , e C , e G , e T }. The latter is a probability of p j on [0, 1]. In the k th iteration of the belief propagation algorithm, message update rules are given by (see, e.g., [23] and the references therein) After the completion of the iterative procedure, the bases b i in the target genome are estimated by first whereb i ∈ {e A , e C , e G , e T }, and then choosing the base with the highestx i value. Note that, by exploiting the symmetry of the expression (4), we can writẽ For the brevity of notation, we denotex Assuming that the prior distribution on p j , D(p j ), is Beta(0,0) (which is essentially as same as the Bernoulli(1/2) distribution), the read confidence is a binary variable, Define a log-likelihood ratio (11) After substituting (5) in (8), we obtain . (9) .
It is straightforward to write y A closer examination of the first element of x i→j (1), leads to simplification shown in (11), where we implicitly used the assumption that p j is binary. We can obtain similar expressions to (11) for other components of x (k) i→j . As a result, the updating rule for x where the vector 2e ij −1 has element 1 in the position corresponding to the nucleotide base b ij represented by e ij , and −1's elsewhere. Therefore, the belief propagation update rule (12) is identical to the update rule (1) of our message passing algorithm presented in Section II. Moreover, update rule (10) is identical (up to the scaling factor) to the message update rule (2). Therefore, message passing scheme proposed in Section II can be interpreted as the belief propagation under a specific prior on the confidence of the aligned data p j -in particular, p j should come from a Beta(0,0) distribution, i.e., be treated as a binary variable.
IV. BINARY REPRESENTATION, MESSAGE PASSING, AND POWER ITERATION ALGORITHM So far, we discussed reference-guided assembly schemes that rely on a representation of the nucleotide basis with 4-dimensional vectors {e A , e C , e G , e T }. As an alternative, in this section we rely on a binary representation of nucleotides to formulate a message passing scheme and discuss the provably convergent power iteration algorithm for finding the target genome sequence. The power iteration scheme finds the desired sequence by computing the leading singular vectors of an appropriately defined data matrix.
The four-letter alphabet {A, C, G, T } in DNA sequencing data can be represented using binary symbols, e.g., {+1, −1}. In particular, we encode the nucleotide basis as A = {−1, −1}, C = {−1, +1}, G = {+1, −1}, and T = {+1, +1}, and represent reads as binary sequences comprising {±1}. Similar to how we built a model utilizing 4-dimensional vectors {e A , e C , e G , e T } in Section II, we define a bipartite graph where each base b i is represented by two binary nodesb 2i−1 andb 2i . Using the output of an alignment algorithm, each read node of the bipartite graph is connected to 2l binary base nodes in the node setb 1:2L , where l denotes read length and L is the length of the target sequence. For convenience, let us denote the resulting graph by G(b 1:2L ∪ r 1:n ,Ẽ). The edge (k, j) inẼ connectingb k and r j is assigned a variable e kj ∈ {±1}, the binary representation ofb k provided by read r j . Given such a graphical representation, we can apply a binary message passing algorithm as in [21]. In particular, the read and base messages are scalars and the update equations are given by After the iterative procedure reaches a stopping criterion, the binary string representing unknown target DNA sequence is obtained as the weighted averagẽ The above algorithm is known to converge to the optimal solution when the bi-partite graph is regular [22]. In our application, however, the graph is not regular since the sequencing coverage varies. Nevertheless, we find that the binary message-passing algorithm performs very well in both simulations and on experimental data, as we demonstrate in Section VI. The binary message passing algorithm is also closely related to the so-called power iteration scheme for computing the leading singular vector of an appropriately defined data matrix. We next examine the power iteration algorithm and argue its convergence.
With the adopted binary encoding of nucleotides, we can represent sequencing reads by a sparse n × 2L matrix D. The 2L columns of D correspond to the L positions in the target sequence whereas the j th row of D comprises binary data representing read r j . In each row, only 2l entries are non-zero (representing an l-long read) while the remaining ones are filled with zeros. Therefore, matrix D has entries D ij ∈ {0, +1, −1}. Since the percentage of nonzero entries of D is 2l L and L l, D is a sparse matrix. It is easy to show (see, e.g., [22]) that if each row of D has the same number of nonzero entries, and the same holds for each column, the left singular vector corresponding to the largest singular value of D is a reliable estimate of the target genome sequence when the measurement noise (i.e., read error rate) is low. Here is an illustration. Let s denote the 2L × 1 binary vector with alphabet {−1, +1} representing the true sequence of length L, and let the number of nonzero entries in each columns of D be c. Consider the case where the reads are error-free and s is a 2L × 1 all one vector 1 2L . Since DD T 1 2L = 2Lc1 2L , then s is an eigenvector of DD T . Here D is a non-negative matrix with entries 0s and 1s and thus, by Perron-Frobenius theorem, 1 2L is a left singular vector corresponding to D's largest singular value. In the general case where s consists of both 1 and −1, we can represent s = S1 2L where S is a 2L × 2L diagonal matrix with diag(S) = s. In this case, it is straightforward to generalize the above analysis and show that s remains to be proportional to the leading singular vector of the matrix D.
Performing singular value decomposition is roughly cubic in the dimension of D and, for our problem dimensions, clearly infeasible. Fortunately, we only need to find u, the leading singular vector of D, and then estimate the target sequence s as sign(u). This can be done in a computationally efficient way using the power iteration technique due to sparsity of D. In particular, the power iteration procedure entails computing To demonstrate convergence of the power iteration scheme (16), let us denote the singular values of D as With a random initialization y (0) , power iterations will converge to the singular vector u if the inequality σ 1 (D) > σ 2 (D) holds strictly. The speed of the convergence of power iterations depends on the ratio σ 2 (D)/σ 1 (D). This can be easily shown by an analysis of the consecutive projections of the iteratively updated vectors x (k) onto the singular vector u. In particular, the projection of x (k) onto u is (u T x (k) )u. A closer look into the singular value decomposition shows that Clearly, power iterations will converge with any initialization if σ 1 (D) > σ 2 (D), and the speed of convergence depends on the ratio of σ 1 (D) and σ 2 (D) -the larger the ratio, the faster the convergence. On the other hand, from (16) it directly follows that the update equations for the entries of x (k) and y (k) can be written as x Note that the power iterations (17) differ from the message update rules (13) and (14) in only one term. As our results in Section VI show, accuracy of message passing and power iterations is essentially identical, while the former converges in significantly fewer iterations than the latter. Moreover, both message-passing schemes -the one based on the representation of basis via 4-dimensional vectors {e A , e C , e G , e T } as well as the one relying on the binary representation of nucleotides -converge after approximately the same number of iterations.
V. BENCHMARKING PERFORMANCE OF THE PROPOSED ASSEMBLY SCHEMES To assess and benchmark performance of the proposed iterative learning schemes, in this section we analyze the probability of error of a genie-aided maximum a posteriori (MAP) estimator of the bases in the target genome sequence. In this problem, "genie-aided" is referring to an idealized scenario where short reads are mapped to the reference genome with no errors, i.e., there are no misplacements of the reads along the reference sequence and the MAP estimator knows exact probabilities of mis-calling the bases in the short reads (i.e., has exact quality score information). Recall that neither our message-passing schemes nor the power iteration algorithm make such practically unrealistic assumptions and, in fact, do not require prior knowledge of quality scores.

A. Genie-aided MAP estimator
Let b k denote the k th base in the target sequence, and let y (i) k denote the signal generated by the sequencing platform as it examines b k , i = 1, 2, . . . , c k , where c k stands for the total number of reads covering b k . Assume that the probability of erroneously calling b k in the i th read is p (i) k . Given the base calls of the reads covering b k , y (i) k , the MAP estimateb k is readily found aŝ where we introduced w . Therefore, the MAP estimateb k formed by combining the information provided by c k reads covering b k is given bŷ In the absence of prior information P (b k = x), the MAP estimation of b k in (18) is identical to the socalled weighted plurality voting [25]. Note that if p (i) k = p for all i and k, (18) becomes the well-known plurality voting scheme.
Intuitively, we expect that the performance of the MAP decision scheme improves as we increase the coverage c k . Note that the above expressions are predicated on the assumption of error-free read mapping.

B. Performance of the genie-aided MAP estimator
For notational convenience, let us write the expression for the estimate in (18) where . 1 The probability of error is defined as To characterize the probability of error of the MAP decision scheme, we rely on the so-called universal generating functions often used in reliability analysis of multi-state systems [27]. Consider n independent discrete random variables X 1 , . . . , X n with probability mass functions (pmf) represented by vectors (x i , p i ) (e.g., P (X i = x ij ) = p ij ). In order to evaluate the pmf of an arbitrary function f (X 1 , . . . , X N ), one has to find the vector y of all the possible values of f (·) and the vector q of the corresponding probabilities. The total number of possible combinations (X 1 , . . . , X N ) is K = n i=1 (k i + 1), where (k i + 1) is the number of different realizations of X i . Since the variables are independent, the probability of each unique combination is equal to the product of the probabilities of the realizations of arguments composing this combination. The probability of the j th combination of the realizations of the variables is q j = n i=1 p ij and the corresponding value of the function is f j = f (x 1j 1 , . . . , x njn ). If different combinations produce the same value of the function, then the probability that f (·) takes that value is equal to the sum of probabilities of the combinations resulting in it. As an illustration, let A h denote the set of combinations resulting in the particular function value f h . If the total number of different values that the function of random variables f (X 1 , . . . , X n ) may assume is H, then its probability mass function is completely specified with a pair of vectors (y, q) defined as A compact representation of the probability mass function of a random variable X i , (x i0 , x i1 , . . . , x ik i , p i0 , p i1 , . . . , p ik i ), is given by a z-transform that takes the polynomial form Such a representation is convenient since, to find the probability that X i ∈ Φ, one can use operator δ defined as Moreover, the z-transform representation enables straightforward calculation of the probability mass function of an arbitrary function f of n independent random variables. This can be facilitated via the composition operator ⊗ f applied to z-transform representations of the probability mass functions of the variables, The technique for finding probability mass functions that relies on the z-transform and composition operators ⊗ f is referred to as the universal z-transform or the universal (moment) generating function (UGF) technique. In the context of this technique, the z-transform of a random variable for which the operator ⊗ f is defined is often referred to as its U-function. For additional background on this subject, we refer an interested reader to [27]. Here, we rely on this technique to characterize the probability of error of the genie-aided sequence assembly scheme. Consider the U-function (similar to (20)) defined for each read position where r i x j denotes the probability that the symbol from read i is x j , and w i (x j ) is the weight associated with the information provided by read i (essentially given by the quality scores, which the genie-aided scheme assumes to be perfectly known). To obtain a U-function of the decision for two positions having respective U-functions U 1 (z) and U 2 (z), the following composition operator can be used, Note that some combinations of v 1 m and v 2 m may lead to the same v 1 m + v 2 m and hence there may be multiple terms in (24)  is at most 4 2 = 16, but may be smaller due to aforementioned grouping of the terms that involve identical vectors.
The above procedure leads to a representation of the probability of error of arriving at the decision for a particular sequence position by combining information provided by two reads. Given an arbitrary subset of reads λ (e.g., so far we discussed λ = {1, 2}), it is straightforward to obtain the U-function for an extended subset λ ∪ j with an arbitrary j / ∈ λ as We can further simplify and arrive at more explicit expressions in the following way. Consider the Ufunction of an arbitrary weighted voting classifier (WVC) over a subset of reads λ, Let W Λ = j∈Λ w j be the total weight of all the votes belonging to the WVC, and let the total weight of the subsystem λ be given by W λ . The weight not belonging to λ can be expressed as Note that if W is the largest element of the vector v λ m and v λ m (W ) − v λ m (i) > σ, then any element v λ m (i) = W can be set to zero since this does not affect the probability of reliability even if all of the remaining votes are given to i. Similarly, vectors satisfying v λ m (W ) − v λ m (1) > σ, W = 1 can be removed from further consideration for the same reason. After these simplifications, the probability of correctly identifying the base is given by where Λ = {1, 2, . . . , c k }. The steps for computing the error probability of decision for a given set of read positions are summarized below.
Having computed the probability of error P (b k = b k ) for a fixed coverage c k (where the MAP estimator formsb k by combining information from c k reads that cover the k th base), we can readily evaluate the probability of error for a random coverage. In particular, the average assembly error probability P error can be found by evaluating Note that the fraction of bases not covered by any read is given by exp(−λ) and thus the overall probability of error is given by P total = 3 4 exp(−λ) + P error . For distributions other than Poisson, e.g., empirical distributions inferred from data, one can still apply the above approach to perform a semi-analytical evaluation of P error .

VI. EXPERIMENTAL RESULTS
In this section we present performance studies using both simulations and experimental data sets. First, using realistic synthetic data, we compare the performance of the message passing algorithm from Section II (Algorithm 1), the binary message passing algorithm and the power iteration algorithm. Moreover, we examine the convergence properties of all these schemes and benchmark their accuracy by comparing it with the genie-aided MAP estimation employed in the idealistic scenario where the exact error probabilities of the reads are known. Then we proceed by testing the algorithms on the experimental data we obtained by sequencing E. Coli and N. Meningitides using the Illumina's HiSeq sequencing instrument that provides 100-bp long reads. In particular, we compare the performance of our developed reference-guided sequence assembly algorithms with the commonly used sequencing data analysis tools including GATK and SAMtools.

A. Simulation data
We simulated reference-guided sequence assembly of the genome of a strain of Neisseria Meningitidis. The reference sequence is obtained from GenBank (http://www.ncbi.nlm.nih.gov/nuccore) database and is L = 2, 184, 406 bases long. The reference is used to generate target sequences having 1% variation rate. We then uniformly select starting positions along the sequence and simulate short reads of length l = 76 (mimicking Illumina's Genome Analyzer II platform). Sequencing errors in these reads are simulated according to the position-dependent base calling error profile typical of this particular sequencing platform [6]. The average error rate of the base calling procedure is 0.015 (averaged over all reads and bases in the reads). To construct the bipartite graphical model, we map the reads to the reference sequence using an alignment algorithm based on the Burrows-Wheeler transform [16] and thus establish connections (i.e., edges) between the read nodes and their aligned base nodes. The read nodes with multiple candidate mapping positions are replicated (where each replica may be assigned different confidence score), and each replica is connected to its corresponding set of base nodes. The bipartite graph with binary base nodes introduced in Section IV is constructed in the same way. We apply both the message passing algorithms from Section II and Section IV to infer the target sequence (note that since the algorithms are randomly initialized, the stopping points and hence the resulting assembled sequences may be different). We also form the binary data matrix representing all the short read data and employ the power iteration method to infer the target genome sequence. While the analysis in Section IV gives a guarantee of convergence of the power iteration algorithm, we found that its convergence is usually faster than the theoretical bound. We set the stopping criterion for all these iterative learning methods as |y j→i | < = 0.01L. It turns out that both message passing algorithms need ∼ 30 iterations to converge, while the power iterations converge in ∼ 50 iterations. We initialize all these algorithm by generating y (0) j→i from Gaussian distribution N (1, 1) and uniform distribution U (0, 1) -our extensive simulation studies indicate that different initializations lead to the same error rate of the considered iterative schemes.
For a comparison, we also consider the plurality voting based decision scheme often used in practice (see, e.g., [26]). Here, multiple calls for a base in any given position along the target sequence are consolidated by performing plurality voting. Notice that, in both message passing and plurality voting, we assume the error profiles of the reads (i.e., base calling error rates) are unknown. Plurality voting assumes all reads has equal reliability while message passing scheme iteratively infers the reliability of each read. We also consider probability of error of the MAP decision scheme in Section III which assumes perfect knowledge of the positions of reads along the target sequence and exact information about position-dependent base calling errors (both assumptions are unrealistic in practice). The error rates of these algorithms are shown in Fig. 3 for various sequencing coverages (horizontal axis shows the average coverage). As can be seen from Fig. 3, the message-passing scheme and the power iteration algorithm outperform plurality voting. The binary message passing algorithm has almost identical accuracy as power iterations, while being slightly worse than Algorithm 1. Moreover, we see that the error rates of message passing are close to the genie-aided MAP decision scheme.

B. Experimental data
In addition to the simulation studies, we tested the performance of our proposed iterative learning schemes for reference-guided sequence assembly using two experimental data sets. In particular, we sequenced Escherichia Coli (from strain MG1655, 4.64 × 10 6 bases long) and Neisseria Meningitidis (from strain FAM18, having length 2.2 × 10 6 ) at the Center for Genomic Sequencing and Analysis of the University of Texas at Austin. The data is obtained using Illumina's HiSeq platform that provides 100bp-long paired-end reads, and the performance of our proposed methods are compared with that of the widely used sequencing analysis packages SAMtools and GATK. Both SAMtools and GATK process aligned next-generation sequencing data stored in SAM format, the alignment file format provided by the majority of frequently used alignment tools (e.g., BWA). These files contain the aligned reads, their positions and the quality scores of the bases. SAMtools calculates empirical quality scores from the alignment information and uses them to recalibrate the raw quality scores provided by the sequencing platform. The assembled sequence is formed using the aligned bases weighted by these new quality scores. In addition to the quality score recalibration, GATK also performs a local realignment procedure to correct misaligned reads, especially from the target genome region containing indels compared to the reference genome. After performing sequence assembly using quality score information, these software packages can also perform downstream single nucleotide polymorphism (SNP) detection, while GATK also incorporates a machine learning tool to separate true variation from sequencing platform artifacts. The two genomes are sequenced using 67% of an HiSeq platform lane having approximately 30 × 10 6 reads, resulting in the coverage greater than 200. This enables accurate inference of the true E. Coli and N. Meningitides sequences using any of the techniques discussed in the paper, providing us with the ground truth. To determine the accuracy of our proposed schemes in realistic scenarios where the coverage is limited, we uniformly subsample the data to emulate low coverage situations. The resulting error rates are shown in Table I. As can be seen there, the developed message passing schemes outperform both SAMtools and GATK in terms of the accuracy. The number of iterations for each message passing scheme was set to 30, which at coverage c = 20 resulted in the average CPU runtimes of 65 and 37 minutes for processing E. Coli and N. Meningitidis data sets, respectively (the algorithms were coded in C++, run on a 3.07G Hz single core machine). The corresponding runtimes for SAMtools are 50 and 28 minutes, and for GATK 53 and 30 minutes. As seen from the table, increasing the coverage can dramatically improve accuracy of the assembly -recall the discussion from Section V where we showed that the probability of error of the genie-aided MAP estimator decreases exponentially with the coverage. However, increasing coverage also increases the cost of the sequencing project.
Note that the sequenced genome might contain insertions as compared to the reference or, equivalently, the reference sequence contains gaps. This structural variation can be detected in the alignment stage by using paired-end reads [29], [30]. The paired-end reads have a known range of lengths of inserts between the reads in a pair. The gaps in the reference can be detected by relying on a multi-read alignment of the pairs of reads and comparing the aligned positions with the insert lengths. We used the scheme in [30] to perform the alignment of our E. Coli data set and detected 34 gaps in the reference. We includes the gap positions as additional base nodes in our graphical model and uses our Algorithm 1 to identify the order of nucleotides in the gaps. As a result, 31 out of 34 gaps were reconstructed (i.e., closed).

VII. SUMMARY AND CONCLUSION
We studied reference-guided sequence assembly from short reads generated by next-generation sequencing technologies, specifically focusing on the problem of obtaining the target genome sequence from potentially erroneous and misaligned reads. We cast the problem as the inference of the target sequence on an appropriately defined bipartite graph and proposed iterative learning algorithms for solving it. In particular, we developed message passing algorithms that rely on both binary as well as representation of nucleotide bases by 4-dimensional vectors. It was shown that the derived message passing algorithm (in particular, Algorithm 1 in Section II) can be interpreted as the standard belief propagation under a certain prior. In addition, the problem was rephrased so that the power iteration algorithm, employed to find the leading singular vector of a matrix collecting all short reads, results in a good approximation of the target sequence. Convergence of power iterations is guaranteed, while the convergence of message passing algorithms is studied empirically. Unlike existing methods, the proposed algorithms find the desired sequence without using reliability information (i.e., quality scores) of the short reads -in fact, message passing algorithms infer the aforementioned quality score information.
To assess achievable accuracy of the proposed iterative learning techniques, we analyzed the probability of error of a genie-aided maximum a posteriori decision scheme in the idealized scenario where the base calling error rates and read mapping locations are known perfectly. It was shown empirically that the iterative learning schemes perform close to the genie-aided estimation scheme, and that they outperform state-of-the-art software packages for downstream processing of sequencing data.