BayesHammer: Bayesian clustering for error correction in single-cell sequencing

Step (1): computing k-mer statistics

To collect k-mer statistics, we use a straightforward hash map approach [12] that does not require storing instances of all k-mers in memory (as excessive amount of RAM might be needed otherwise). For a certain positive integer N (the number of auxiliary files), we use a hash function h: ∑ ^k →ℤ _N that maps k-mers over the alphabet Σ to integers from 0 to N - 1.

Algorithm 1 Count k-mers

for each k-mer x from the reads R: do

compute h(x) and write x to File_h(x).

for i ∈ [0, N - 1]: do

sort File _i with respect to the lexicographic order;

reading File _i sequentially, compute statistics(s) for each k-mer s from File _i .

Step (2): constructing connected components of Hamming graph

For a set of k-mers X, the Hamming graph HG _τ (X) is an undirected graph with the set of vertices X and edges corresponding to pairs of k-mers from X with Hamming distance at most τ, i.e., x, y ∈ X are connected by an edge in HG _τ (X) iff d(x, y) ≤ τ (Figure 3). To construct HG _τ (X) efficiently, we notice that if two k-mers are at Hamming distance at most τ, and we partition the set of indices [0,k - 1] into τ + 1 parts, then at least one part corresponds to the same subsequence in both k-mers. Below we assume with little loss of generality that τ + 1 divides k, i.e., k = σ (τ + 1) for some integer σ.

For a subset of indices I ⊆ [0, k - 1], we define a partial lexicographic ordering ≺ _I as follows: x ≺ _I y iff x[I] ≺ y[I], where ≺ is the lexicographic ordering on Σ*. Similarly, we define a partial equality = _I such that x = _I y iff x[I] = y[I]. We partition the set of indices [0, k - 1] into τ + 1 parts of size σ and for each part I, sort a separate copy of X with respect to ≺ _I . As noticed above, for every two k-mers x, y ∈ X with d(x, y) ≤ τ, there exists a part I such that x = _I y. It therefore suffices to separately consider blocks of equivalent k-mers with respect to = _I for each part I. If a block is small (i.e., of size smaller than a certain threshold), we go over the pairs of k-mers in this block to find those with Hamming distance at most τ. If a block is large, we recursively apply to it the same procedure with a different partition of the indices. In practice, we use two different partitions of [0, k - 1]: the first corresponds to contigious subsets of indices (recall that $\sigma =\frac{k}{\tau +1}$):

Algorithm 2 Hamming graph processing

procedure HGPROCESS(X, max_quadratic)

Init components with singletons $\mathcal{X}=\left\{\left\{x\right\}:x\in X\right\}$.

for all ϒ ∈ FindBlocks$\left(X,{\left\{I_s^{\mathsf{\text{cnt}}}\right\}}_{s=0}^{\tau }\right)$do

if |ϒ| > max_quadratic then

   for all Z ∈ FindBlocks $\left(\Upsilon ,{\left\{I_s^{\mathsf{\text{str}}}\right\}}_{s=0}^{\tau }\right)$do

      ProcessExhaustively$\left(Z,\mathcal{X}\right)$

else

   ProcessExhaustively$\left(\Upsilon ,\mathcal{X}\right)$.

function FindBlocks $\left(X,{\left\{I_s\right\}}_{s=0}^{\tau }\right)$

for s = 0,...,τ do

sort a copy of X with respect to ${\prec }_{I_s}$, getting X _s .

for s = 0,...,τ do

output the set of equiv. blocks $\left\{\Upsilon \right\}\mathsf{\text{w}}.\mathsf{\text{r}}.\mathsf{\text{t}}.{=}_{I_s}$.

procedure PROCESS EXHAUSTIVELY$\left(\Upsilon ,\mathcal{X}\right)$

for each pair x, y ∈ ϒ do

if d(x, y) ≤ τ then join their sets in $\mathcal{X}$:

   for all $x\in Z_x\in \mathcal{X},y\in Z_y\in \mathcal{X}$ do

BAYES HAMMER uses a two-step procedure, first splitting with respect to ${\left\{I_s^{\mathsf{\text{cnt}}}\right\}}_{s=0}^{\tau }$ (Figure 4) and then, if an equivalence block is large, with respect to ${\left\{I_s^{\mathsf{\text{str}}}\right\}}_{s=0}^{\tau }$. On the block processing step, we use the disjoint set data structure [12] to maintain the set of connected components. Step (2) is summarized in Algorithm 2.

Step (3): Bayesian subclustering

Under this model, the maximum likelihood center of a cluster is simply its consensus string [8].

In BAYES HAMMER, we further elaborate upon HAMMER's model. Instead of a fixed ε, we use reads quality values that approximate probabilities q _x [i] of a nucleotide at position i in the k-mer x being erroneous. We combine quality values from identical k-mers in the reads: for a multiset of k-mers X that agree on the j^th nucleotide, it is erroneous with probability Π_x∈Xq _x [j].

where c _i is the center (consensus string) of the subcluster C _i .

In the subclustering procedure (see Algorithm 3), we sequentially subcluster each connected component of the Hamming graph into more and more clusters with the classical k-means clustering algorithm (denoted m-means since k has different meaning). For the objective function, we use the likelihood as above penalized for overfitting with the Bayesian information criterion (BIC) [13]. In this case, there are |C| observations in the dataset, and the total number of parameters is 3 km + m - 1:

Algorithm 3 Bayesian subclustering

for all connected components C of the Hamming graph do

m := 1

ℓ₁ := 2 log L₁(C) (likelihood of the cluster generated by the consensus)

repeat

m := m + 1

do m-means clustering of C = C₁ ∪...∪ C _m w.r.t. the Hamming distance; the initial approximation to the centers is given by k-mers that have the least error probability

ℓ _m := 2 · log L _m (C₁,...,C _m ) (3 km + m - 1) · log |C|

until ℓ _m ≤ ℓ _m-1

output the best found clustering C = C₁ ∪...∪ C_m-1

for subclustering into m clusters; we stop as soon as ℓ _m ceases to increase.

Steps (4) and (5): selecting solid k-mers and expanding the set of solid k-mers

In contrast to HAMMER, we do not distinguish whether the cluster is a singleton (i.e., |C| = 1); there may be plenty of superfluous clusters with several k-mers obtained by chance (actually, it is more likely to obtain a cluster of several k-mers by chance than a singleton of the same total multiplicity).

Initially we mark as solid the centers of the clusters whose total quality exceeds a predefined threshold (a global parameter for BAYES HAMMER, set to be rather strict). Then we expand the set of solid k-mers iteratively: if a read is completely covered by solid k-mers we conclude that it actually comes from the genome and mark all other k-mers in this read as solid, too (Algorithm 4).

Step (6): reads correction

After Steps (1)-(5), we have constructed the set of solid k-mers that are presumably error-free. To construct corrected reads from the set of solid k-mers, for each base of every read, we compute the consensus of all solid k-mers and solid centers of clusters of all non-solid k-mers covering this base (Figure 5). This step is formally described as Algorithm 5.

Algorithm 4 Solid k-mers expansion

procedure ITERATIVE EXPANSION(R, X)

while ExpansionStep(R, X) do

function EXPANSION STEP(R, X)

for all reads r ∈ R do

if r is completely covered by solid k-mers then

   mark all k-mers in r as solid

Return TRUE if X has increased and FALSE otherwise.

Algorithm 5 Reads correction

Input: reads R, solid k-mers X, clusters $\mathcal{C}$.

for all reads r ∈ R do

init consensus array υ: [0, |r| - 1] × {A, C, G, T} → ℕ with zeros: υ(j, x[i]):= 0 for all i = 0,...,|r| - 1 and j = 0,...,k - 1

for i = 0,...,|r| - k do

if r[i, i + k - 1] ∈ X (it is solid) then

   for j ∈ [i, i + k - 1] do

      υ(j, r[i]):= υ(j, r[i]) + 1

if r[i, i + k - 1] ∈ C for some C ∈ $\mathcal{C}$then

   let x be the center of C

   if x ∈ X (r belongs to a cluster with solid center) then

      for j ∈ [i, i + k - 1] do

         υ(j, x[i]):= υ(j, x[i]) + 1

for i ∈ [0, |r| - 1] do

r[i]:= arg max_a∈Συ(i, a).