Haploid to diploid alignment for variation calling assessment

Abstract

Motivation

Variation calling is the process of detecting differences between donor and consensus DNA via high-throughput sequencing read mapping. When evaluating the performance of different variation calling methods, a typical scenario is to simulate artificial (diploid) genomes and sample reads from those. After variation calling, one can then compute precision and recall statistics. This works reliably on SNPs but on larger indels there is the problem of invariance: a predicted deletion/insertion can differ slightly from the true one, yet both make the same change to the genome. Also exactly correct predictions are rare, especially on larger insertions, so one should consider some notion of approximate predictions for fair comparison.

Results

We propose a full genome alignment-based strategy that allows for fair comparison of variation calling predictions: First, we apply the predicted variations to the consensus genome to create as many haploid genomes as are necessary to explain the variations. Second, we align the haploid genomes to the (aligned) artificial diploid genomes allowing arbitrary recombinations. The resulting haploid to diploid alignments tells how much the predictions differ from the true ones, solving the invariance issues in direct variation comparison. In an effort to make the approach scalable to real genomes, we develop a simple variant of the classical edit distance dynamic programming algorithm and apply the diagonal doubling technique to optimise the computation. We experiment with the approach on simulated predictions and also on real prediction data from a variation calling challenge.

Background

In the study of human genetics, variation calling from high-throughput sequencing reads [1] is a revolutionary technique. Conceptually, the process is remarkably simple. Sequence random short fragments from donor DNA and align them to the reference (consensus) genome. Outside repetitive regions a good alignment is unique, hence if the resulting multiple alignment (read pileup) has columns where reads vote for something differing from the reference genome, the donor is very likely to actually contain this variant in his/her genome. Numerous methods have been proposed to fine-tune this simple scheme. One could argue that this problem is actually an enormous local multiple alignment problem, and so it is not surprising that methods trying to locally improve the alignment are able to improve the accuracy. This standard approach only captures small scale variations, and the methods for discovering large insertions, deletions, translocations, reversals, etc., are more involved [2].

Surprisingly there has been no systematic approach for studying the accuracy issue of variation calling predictions. With real data sets one can always resort to Sanger sequencing to validate the findings, but this is expensive. With simulated data sets one can compare to the ground truth and compute precision/recall statistics. A typical simulated data set is produced by first applying random mutations to a reference genome or selecting a random subset of a predefined set of known frequent variations in the population (or a mixture of these), to create an artificial donor genome. We consider here the setting of a diploid genome, where this simulation is conducted twice, taking into account heterozygousity/homozygousity information on the variants to create a realistic diploid genome. What we assume is that the alignment of the diploid genome to the reference is preserved. Then one can simulate a random pool of short reads from the artificial diploid genome with high enough coverage such that variation calling is plausible. Once a list of predicted variants is produced, one can match it to the list of variants used for creating the artificial diploid.

However, this direct comparison has the shortcoming that invariances among predictions are not properly considered. To see this, consider reference ACGGAAGGT, donor ACGGT, and let the simulated real variant be deletion of GAAG at position 3 (starting indexing from 0). Predicted deletion AAGG at position 4 results in the same donor, but pair-wise comparison of predicted and real variants would not reveal this. While these kind of simple invariances are easy to take into account, things get much more involved with split/merged predictions of various types. Clearly some kind of dynamic programming approach would be required to find best editing of the predicted variants to make them match the real ones. Here we assume that the allowed variant types are only insertions, deletions and substitutions, since rearrangements can always be modelled as series of the former type of operations.

Instead of correcting for the invariances in the manner sketched above, we propose to resort to a natural variant of the familiar sequence level full-genome alignment. The approach is as follows. First apply the predicted variations to the consensus genome to create as many haploid genomes that are necessary to explain the variations. Heterozygous variations are only applied on one haploid and homozygous variations on all. Second, align the haploid genomes to the (aligned) artificial diploid genomes allowing arbitrary recombinations. The intuition is that with short read mapping one cannot usually obtain much haplotype level information, but only determine whether the variations are heterozygous or homozygous. The haploid genome created by applying all non-overlapping predicted variants is better the closer it is to an arbitrary recombination of the diploid alleles. This is what the approach measures by edit distance. If there are overlapping predicted heterozygous variants, multiple haploid genomes are created and compared with a (different) recombination of the diploid alleles. Hence, the resulting haploid to diploid alignments tell how much the predictions differ overall from the true ones, solving the invariance issues in direct variation comparison. To make the approach scalable to real genomes, we develop a simple variant of the classical edit distance dynamic programming algorithm and apply the diagonal doubling technique to optimise the computation. We experiment with the approach on simulated predictions and on real prediction data from a variation calling challenge [3].

Methods

Problem

Let B¹ and B² represent diploid genome sequences generated by applying mutations to a reference sequence P. Let tables m¹ and m² store for each position in B¹ and B² the corresponding position in P defined by the multiple alignment of P, B¹, and B², which in turn is defined by how P is mutated to B¹ and B². For example, the multiple alignment

   0123456 7 8

P  AGCTGAT-A-C

B1 ACCTGATCACG

B2 ACCTGCT-ACC

is defined by homozygous substitution G->C at 1st position of P , heterozygous substitution A->C at 5th position of P , heterozygous insertion of C after 6th position of P , homozygous insertion of C after 7th position of P, and heterozygous substitution C->G at 8th position of P. The corresponding mappings are

m1 01234566778

m2 0123456778

After simulating reads from B¹ and B², let some variation calling method end up with the above list of variations. Then one can construct a haploid genome A by applying greedily all predicted mutations to P from left to right. In the general case, one could be left with overlapping mutations. Then one should construct another haploid genome applying all homozygous mutations as well; it is possible that there are more overlaps because of prediction errors [4], but since our simulated ground-truth is diploid, generating more predicted genomes is not beneficial in this case.

In our running example, all mutations can be applied simultaneously, as there are no overlapping heterozygous mutations. The result is A = ACCTGCTCACC, which can be aligned to the above multiple alignment with no errors by allowing the alignment to recombine B¹ and B² at the positions they share in common with P. We call such a recombination a valid reference guided recombination. The result is shown below (lower case showing the alignment of A).

A  ACCTGCTCACC

B1 ACCTGATcaCG

B2 acctgct-Acc

Invariance

To highlight the issue of invariance, let us consider a different reference P = GATCAATGAG, a single diploid B and a haploid A given by some variation calling method. Let us assume that B and A differ from P by the following variations

$$\begin{array}{ccc}\text{position}\hfill & \text{P}\hfill & \text{B}\hfill \\ 0\hfill & \text{GA}\hfill & \text{G}\hfill \\ 3\hfill & \text{C}\hfill & \text{CC}\hfill \\ 4\hfill & \text{A}\hfill & \text{C}\hfill \\ 7\hfill & \text{G}\hfill & \text{A}\hfill \end{array}\begin{array}{ccc}\text{position}\hfill & \text{P}\hfill & \text{A}\hfill \\ 1\hfill & \text{A}\hfill & \text{T}\hfill \\ 2\hfill & \text{T}\hfill & \text{C}\hfill \\ 3\hfill & \text{C}\hfill & \text{CC}\hfill \\ 3\hfill & \text{CA}\hfill & \text{C}\hfill \\ 7\hfill & \text{G}\hfill & \text{A}\hfill \end{array}$$

The evaluation method used in Variathon 2013 [3] finds only two common variations: the one in position 3 and the one in position 7. Even if one was to use a method to match two invariant variants that make the same effect to the reference, this would not help in this example; only when the variants are applied altogether, the result is the same genome GTCCCATAAG:

P  GATC-AATGAG          P  GATC-AATGAG

B  G-TCCCATAAG          A  GTCCC-ATAAG

To allow approximate predictions and account for invariance, we propose to compute the unit cost (Levenshtein) edit distance between A and any recombination B¹ and B² according to their alignment.

The computation of such edit distance is an easy extension of the standard dynamic programming, and has been studied earlier under name jumping alignments in the context of amino acid sequences [5]. However, global alignment of full genomes is infeasible using standard quadratic time dynamic programming. Therefore we extend the diagonal doubling optimisation [6] to compute the edit distance in O (dn) time, where d is the least edit distance between the haploid and any valid reference guided recombination of the diploid genomes, and n is the maximum of the input sequence lengths.

Algorithm

We simultaneously fill two dynamic programming matrices $\left(d_{i,j_1}^1\right)$ and $\left(d_{i,j_2}^2\right)$ such that $\left(d_{i,j_1}^1\right)$ (resp. $\left(d_{i,j_2}^2\right)$) gives the minimum edit distance between A[0 ... i] and R[0 ... j'] such that R[0 ... j'] is a valid reference guided recombination of B¹[0 ... j₁] and B²[0 ... j₂] ending at B¹ (resp. ending at B²). Filling the matrices is an easy extension of standard dynamic programming for edit distance: take minimum across the matrices at columns j₁ and j₂ such that B¹ and B² can recombine. Such j₁ and j₂ values have the property that

$$\begin{array}{c}{m}^1\left[j_1-1\right]=m^2\left[j_2-1\right]=j\text{or}\\ m^1\left[j_1-1\right]=m^2\left[j_2\right]-1=j\end{array}$$

and

$$m^1\left[j_1-2\right]\ne j\ne m^2\left[j_2-2\right]$$

This results into algorithm 1, where call d = diploid_align(A, B 1, B 2, m¹, m², reference length) returns the minimum edit distance d between A and any valid reference guided recombination of B¹ and B².

Proof of correctness for Algorithm 1

Definition 1. Let A and B be strings. Now ed(A, B) is the Levenshtein edit distance between the strings A and B.

Definition 2. Let A, B¹ and B² be strings. R _j is the reference guided recombination of B¹[0 ... k] and B²[0 ... v] that minimises ed(A, R _j ), where k = max{i | m¹[i] ≤ j} (likewise for v). $R_j^1$ is as R _j but the last character must be from the string B¹. $R_j^2$ is defined similarly.

Lemma 1. Let A, B¹ and B² be strings and R _j , $R_j^1$ and $R_j^2$ as above. Now

$$\text{ed}\left(A_i,R_j\right)=\text{min}\left\{\begin{array}{c}\text{ed}\left(A_i,R_j^1\right)\hfill \\ \text{ed}\left(A_i,R_j^2\right).\hfill \end{array}\right.$$

Algorithm 1 Haploid to diploid alignment algorithm.

  function DIPLOID _ALIGN(A, B¹, B², m¹, m², reference_length)

   d¹[length(B¹) + 1][length(A) + 1]

   d²[length(B²) + 1][length(A) + 1]

    /* Initialise as in Levenshtein distance. */

  function CALCULATE_COLUMN(matrix, col, B_char)

    for i ← 1, length(A) do     

$$matrix\left[col\right]\left[i\right]\leftarrow \text{min}\left\{\begin{array}{c}matrix\left[col-1\left[i\right]+1\right]\hfill \\ matrix[col-1\left[i-1\right]+\left[A\left[i\right]\ne B\text{\_}char\right]\hfill \\ matrix\left[col\right]\left[i-1\right]+1\hfill \end{array}\right.$$

  function MIN_FROM_TO(from, to)

    for i = 1 → length(A) do     

$$to\left[i\right]\leftarrow \text{min}\left\{\begin{array}{c}to\left[i\right]\hfill \\ from\left[i\right]\hfill \end{array}\right.$$

  column1 ← 1

  column2 ← 1

  for j ← 0, reference_length − 1 do

    if m¹[column1 − 1] = j then

    CALCULATE_COLUMN(d¹, column1, B¹[column1 − 1])

     column1 ← column1 + 1

    if m²[column2 − 1] = j then

     CALCULATE_COLUMN(d², column2, B²[column2 − 1])

      column2 ← column2 + 1

    if m¹[column1 − 2] = j and (m²[column2 − 1] = j + 1 or m²[column2 − 2] = j) then

     MIN_FROM_TO(d¹[column1 − 1], d²[column2 − 1])

    if m²[column2 − 2] = j and (m¹[column1 − 1] = j + 1 or m¹[column1 − 2] = j) then

     MIN_FROM_TO(d²[column2 − 1], d¹[column1 − 1])

    while m¹[column1 − 1] ≤ j do

     CALCULATE_COLUMN(d¹, column1, B¹[column1 − 1])

      column1 ← column1 + 1

    while m²[column2 − 1] ≤ j do

     CALCULATE_COLUMN(d², column2, B²[column2 − 1])

      column2 ← column2 + 1 

$$return\text{min}\left\{\begin{array}{c}{d}^1\left[\text{length}\left(B^1\right)\right]\left[\text{length}\left(A\right)\right]\hfill \\ d^2\left[\text{length}\left(B^2\right)\right]\left[\text{length}\left(A\right)\right]\hfill \end{array}\right.$$

Lemma 2. When a recombination

$$B^2\left[0\dots \upsilon -2\right]B^1\left[k-1\dots \right]$$

is possible, that is,

$$m^2\left[\upsilon -2\right]+1=m^1\left[k-1\right]\text{or}{m}^2\left[\upsilon -2\right]=m^1\left[k-2\right],$$

it holds that

$$\begin{array}{c}\text{min}\left\{\begin{array}{c}\text{min}\left\{\begin{array}{c}\text{ed}\left(A_i,R_{m^1\left[k-2\right]}^1\right)\hfill \\ \text{ed}\left(A_i,R_{m^2\left[\upsilon -2\right]}^2\right)\hfill \end{array}\right.+1\hfill \\ \text{min}\left\{\begin{array}{c}\text{ed}\left(A_{i-1},R_{m^1\left[k-2\right]}^1\right)\hfill \\ \text{ed}\left(A_{i-1},R_{m^2\left[\upsilon -2\right]}^2\right)\hfill \end{array}\right.+\delta \hfill \\ \text{ed}\left(A_{i-1},R_{m^1\left[k-1\right]}^1\right)+1\hfill \end{array}\right.\\ =\text{ed}\left(A_i,R_{m^1\left[k-1\right]}^1\right)\hfill \end{array}$$

where δ = 1 if A[i] ≠ B¹[k − 1] and 0 otherwise.

Proof. First we observe that

$$m^{\text{2}}\left[\upsilon -\text{2}\right]+\text{1}=m^{\text{1}}\left[k-\text{1}\right]$$

implies

$$m^1\left[k-\text{2}\right]\le m^{\text{2}}\left[\upsilon -\text{2}\right]=m^{\text{1}}\left[k-\text{1}\right]-\text{1}.$$

Using this we get

$$\begin{array}{c}\text{min}\left\{\begin{array}{c}\text{min}\left\{\begin{array}{c}\text{ed}\left(A_i,R_{m^1\left[k-2\right]}^1\right)\hfill \\ \text{ed}\left(A_i,R_{m^2\left[\upsilon -2\right]}^2\right)\hfill \end{array}\right.+1\hfill \\ \text{min}\left\{\begin{array}{c}\text{ed}\left(A_{i-1},R_{m^1\left[k-2\right]}^1\right)\hfill \\ \text{ed}\left(A_{i-1},R_{m^2\left[\upsilon -2\right]}^2\right)\hfill \end{array}\right.+\delta \hfill \\ \text{ed}\left(A_{i-1},R_{m^1\left[k-1\right]}^1\right)+1\hfill \end{array}\right.\\ =\text{min}\left\{\begin{array}{c}\text{min}\left\{\begin{array}{c}\text{ed}\left(A_i,R_{m^2\left[\upsilon -2\right]}^1\right)\hfill \\ \text{ed}\left(A_i,R_{m^2\left[\upsilon -2\right]}^2\right)\hfill \end{array}\right.+1\hfill \\ \text{min}\left\{\begin{array}{c}\text{ed}\left(A_{i-1},R_{m^2\left[\upsilon -2\right]}^1\right)\hfill \\ \text{ed}\left(A_{i-1},R_{m^2\left[\upsilon -2\right]}^2\right)\hfill \end{array}\right.+\delta \hfill \\ \text{ed}\left(A_{i-1},R_{m^1\left[k-1\right]}^1\right)+1\hfill \end{array}\right.\end{array}$$

and by Lemma 1

$$\begin{array}{c}=\text{min}\left\{\begin{array}{c}ed\left(A_i,{R_{m^2}}_{\left[\upsilon -2\right]}\right)+1\hfill \\ ed\left(A_{i-1},R_{m^2\left[\upsilon -2\right]}\right)+\delta \hfill \\ ed\left(A_{i-1},{R^1}_{m^1\left[k-1\right]}\right)+1\hfill \end{array}\right.\\ =ed\left(A_i,{R^1}_{m^1\left[k-1\right]}\right)\end{array}$$

which is as claimed. Note that if m²[v − 2] = m¹[k − 2], the claim follows directly without the observation.

Lemma 3. If the strings B¹ and B² do not overlap, that is

$$\begin{array}{c}{m}^{\text{1}}\left[\text{length}\left(B^{\text{1}}\right)-\text{1}\right]<m^{\text{2}}\left[0\right]\text{or}\hfill \\ m^{}\left[\text{length}\left(B^{\text{2}}\right)-\text{1}\right]<m^{\text{1}}\left[0\right]\hfill \end{array}$$

then we can compute the diploid_align as the Levenshtein distance between the string A and the appropriate concatenation B¹B² or B²B¹.

Theorem 1. The Algorithm 1 calculates

$$\begin{array}{c}{d}^1\left[k\right]\left[i\right]=\text{ed}(A_i,R_{m^1\left[k-1]\right)}^{{}_1}1)\text{and}\\ d^2\left[v\right]\left[i\right]=\text{ed}\left(A_i,R_{m^2\left[\upsilon -1]\right)}^2\right).\end{array}$$

Proof. The proof is by induction over the index sum i + k for the case

$$d^{\text{1}}\left[k\right]\left[i\right]=\text{ed}\left(A_i,{R^{\text{1}}}_{m^{\text{1}}\left[k-\text{1]}\right)}\right).$$

The other case is symmetric.

First consider the base cases. The algorithm initialises the matrices as in the Levenshtein distance algorithm. By Lemma 3 we can assume that the strings B¹ and B² overlap and thus no recombination needs to be considered before the first iteration of the algorithm. The base cases are thus as required. Let us assume that the claim holds for i + k < p + q. If the condition

$$\begin{array}{c}{m}^{\text{2}}\left[column2-\text{2}\right]==j\text{and}\\ (m^1\left[column1-1\right]==j+1\text{or}\\ m^1\left[column1-1\right]==j)\end{array}$$

(1)

held on a previous iteration, then the minimum of the columns d¹[q − 1] and d²[column2 − 1] was assigned to d¹[q − 1]. It also means that

$$m^{\text{2}}\left[column2-\text{2}\right]+\text{1}=m^{\text{1}}\left[q-\text{1}\right].$$

(2)

Let δ = 1 if A[p] ≠ B¹[q − 1] and 0 otherwise. Now the algorithm calculates

$$\begin{array}{c}{d}^1\left[q\right]\left[p\right]=\text{min}\left\{\begin{array}{c}\text{min}\left\{\begin{array}{c}{d}^1\left[q-1\right]\left[p\right]\hfill \\ d^2\left[column2-1\right]\left[p\right]\hfill \end{array}\right.+1\hfill \\ \text{min}\left\{\begin{array}{c}{d}^1\left[q-1\right]\left[p-1\right]\hfill \\ d^2\left[column2-1\right]\left[p-1\right]\hfill \end{array}\right.+\delta \hfill \\ d^1\left[q\right]\left[p-1\right]+1\hfill \end{array}\right.\\ =\text{min}\left\{\begin{array}{c}\text{min}\left\{\begin{array}{c}\text{ed}\left(A_p,R_{m^1\left[q-2\right]}^1\right)\hfill \\ \text{ed}\left(A_p,R_{m^2\left[column2-2\right]}^2\right)\hfill \end{array}\right.+1\hfill \\ \text{min}\left\{\begin{array}{c}\text{ed}\left(A_{p-1},R_{m^1\left[q-2\right]}^1\right)\hfill \\ \text{ed}\left(A_{p-1},R_{m^2\left[column2-2\right]}^2\right)\hfill \end{array}\right.+\delta \hfill \\ \text{ed}\left(A_{p-1},R_{m^1\left[q-1\right]}^1\right)+1\hfill \end{array}\right.\\ =\text{ed}\left(A_p,R_{m^1\left[q-1\right]}^1\right)\end{array}$$

where the first equality holds by the induction assumption. The second equality holds as by equation (2) we can use Lemma 2 with the assignment v := column2 and k := q.

If condition (1) did not hold on the previous iteration, no minimisation was performed. Thus the algorithm calculates

$$\begin{array}{c}{d}^1\left[q\right]\left[p\right]=\text{min}\left\{\begin{array}{c}{d}^1\left[q-1\right]\left[p\right]+1\hfill \\ d^1\left[q-1\right]\left[p-1\right]+\delta \hfill \\ d^1\left[q\right]\left[p-1\right]+1\hfill \end{array}\right.\\ =\text{min}\left\{\begin{array}{c}\text{ed}\left(A_p,{R^1{}_{m^1}}_{\left[q-2\right]}\right)+1\hfill \\ \text{ed}\left(A_{p-1},{R^1}_{m^2\left[q-2\right]}\right)+\delta \hfill \\ \text{ed}\left(A_{p-1},{R^1}_{m^1\left[q-1\right]}\right)+1\hfill \end{array}\right.\\ =\text{ed}\left(A_p,{R^1}_{m^1\left[q-1\right]}\right)\end{array}$$

which is correct, as no recombination could be made at this point.

Diagonal doubling

The runtime optimisation known as diagonal doubling [6] can be used also with Algorithm 1. The idea of the original algorithm is as follows. Consider checking whether the Levenshtein edit distance between two sequences of length n is at most k, where k is a given cutoff. Then it is enough to do computation on diagonals i − j such that |i − j| ≤ [k/2], since every change of diagonal costs 1; any path visiting a cell outside this diagonal zone results into edit script with cost more than k. With unknown cutoff, it is easy to see that starting with k = 1, doubling the value of k at each iteration, recomputing the diagonal zone, and stopping the doubling when the computed distance value remains unchanged, results in an O (dn) time algorithm, where d is the edit distance. The idea works with technical changes for two sequences of different length. With some care, one can perform the computation in O (d) space.

To modify the optimisation for our haploid to diploid alignment, the following details need to be taken into account. Because the two columns calculated by the algorithm at every iteration may not completely overlap, one needs to make sure that the minimisation step is done correctly. The minimum of the two columns can only be taken for the overlapping part. Because of the two Levenshtein calculations, we also need to choose the starting cutoff values with care. In particular one needs to start with max(| length(A) − length(B¹)|, | length(A) − length(B²)|), so that the diagonal is big enough to cover the bottom right corner in both matrices.

The diagonal doubling does not affect the correctness of the Algorithm 1. The same argument as in the case of Levenshtein distance applies here. Let us assume that we are interested only on edit distances less than a certain cutoff value. As both matrices contain a Levenshtein distance, it must be that the values grow monotonically when we deviate from the diagonal. Thus, when taking the minimum between the two columns, the values outside the overlapping area must be greater than the current cutoff value. This means that the optimum edit distance must derive from the values in the overlapping area.

Results and discussion

To test the algorithm, we implemented it in C with the diagonal doubling optimisation [6]. The diploid_align algorithm was first tested with generated data. Two diploid strings were generated from a DNA string of appropriate length. Single inserts, deletes and substitutions were applied each with probability of 0.003, and with probability 0.5 they were applied to both diploids. This gives approximately n/ 100 variations in a string of length n. The algorithm was then run five times with the original string as the haploid input and an average of the runtimes was taken. The test machine was a laptop with a Intel Core 2 Duo T7300 2GHz processor running Linux. The results are shown in the Table 1. The actual distance is always about half the generated amount of variations, as nearly half the variations are heterozygous and in those cases one allele has the same base as the haploid input (original sequence).

Table 1 Results of running diploid_align with generated data.

The second test used artificially generated diploids from human chromosome 20 created for a pilot variation calling challenge, Variathon 2013 [3]. From each of the two submissions taking part in the challenge we created one predicted haploid by a simple script that applied the predicted variants to the reference genome. These predicted haploids were then aligned to the diploid. Table 2 shows the original evaluations from the challenge (first 6 columns) and our new evaluation with edit distance. As can be seen, our evaluation agrees with the evaluation done for the challenge.

Table 2 Results with predicted haploids from Variathon.

We also ran a small experiment to highlight the issue of invariance. With the example given in Section Methods, Invariance, the evaluation method used in Variathon 2013 [3] finds only two common variations, as claimed. Not surprisingly, the edit distance between the diploid and the haploid is zero, as in both cases the variations when combined produce the same genome GTCCCATAAG.

Conclusions

We proposed an approach for assessing variation calling predictions in the case of artificial diploid genomes, using a modification of global alignment with a diagonal doubling optimisation to compute it on large inputs. The motivation for the approach was to avoid the invariance problems of direct variation comparison.

We tested the approach on a haploid instance to show its robustness and scalability to complete (small) genomes. We also compared the output of our algorithm with applicable evaluations from a variation calling challenge and the results agreed.