Sequential and parallel algorithms for all-pair $k$-mismatch maximal common substrings

Highlights

• First and only known sub-quadratic algorithms for the k-mismatch maximal common substring problem.

• Sequential algorithm runs in $O(N{log}^{k}N+\mathsf{occ})$ expected time, where occ is the output size.

• Our distributed memory parallel algorithm runs in $O\left(\left(\left(\frac{N}{p}\right)logN+\mathsf{occ}\right){log}^{k}N\right)$ expected time using $O\left({log}^{k+1}N\right)$ expected communication rounds, where p is the number of processors.

• Parallel algorithm is demonstrated with high throughput genomic datasets ranging in size from 18 million to over 270 million reads, on up to 1024 processor cores.

Abstract

Identifying long pairwise maximal common substrings among a large set of sequences is a frequently used construct in computational biology, with applications in DNA sequence clustering and assembly. Due to errors made by sequencers, algorithms that can accommodate a small number of differences are of particular interest. Formally, let $\mathcal{D}$ be a collection of $n$ sequences of total length $N$, $\varphi$ be a length threshold, and $k$ be a mismatch threshold. The goal is to identify and report all $k$-mismatch maximal common substrings of length at least $\varphi$ over all pairs of strings in $\mathcal{D}$. Heuristics based on seed-and-extend style filtering techniques are often employed in such applications. However, such methods cannot provide any provably efficient run time guarantees. To this end, we present a sequential algorithm with an expected run time of $O(N{log}^{k}N+\mathsf{occ})$, where $\mathsf{occ}$ is the output size. We then present a distributed memory parallel algorithm with an expected run time of $O\left(\left(\frac{N}{p}logN+\mathsf{occ}\right){log}^{k}N\right)$ using $O\left({log}^{k+1}N\right)$ expected rounds of global communications, under some realistic assumptions, where $p$ is the number of processors. Finally, we demonstrate the performance and scalability of our algorithms using experiments on large high throughput sequencing data.

Introduction

Sequence matching algorithms are at the core of many applications in computational biology. Next Generation Sequencing (NGS) [15] instruments sequence hundreds of millions of short reads, that are typically randomly sampled from one or multiple genomes. Deciphering pairwise alignments between the reads is often the first step in many applications. For example, one may be interested in finding all pairs of reads that have a sufficiently long overlap, such as suffix/prefix overlap (for genomic or metagenomic assembly [20]), or substring overlap (for read compression [8], finding RNA sequences containing common exons [9], [17], etc.). Much of modern-day high-throughput sequencing is carried out using Illumina sequencers, which have a small error rate ($<$ 1%–2%) and predominantly ($>99\text{\%}$) substitution errors [16]. Thus, algorithms that tolerate a small number of mismatch errors can yield the same solution as the much more expensive alignment computations. Motivated by such applications, we formulate the following all-pair $k$-mismatch maximal common substrings problem:

Problem 1

Given a collection $\mathcal{D}=\{S_1,S_2,\dots ,S_n\}$ of $n$ sequences with ${\sum }_i\left|S_i\right|=N$, a length threshold $\varphi$, and a mismatch threshold $k\ge 0$, report all $k$-mismatch maximal common substrings of length $\ge \varphi$ between all pairs of sequences in $\mathcal{D}$.

A pair of two equal length substrings $S_i[x..(x+t-1)]$ and $S_j[y..(y+t-1)]$ is a $k$-mismatch common substring if the hamming distance between them is $\le k$. Also, they are a $k$-mismatch maximal common substring if neither $S_i[(x-1)..(x+t-1)]$ and $S_j[(y-1)..(y+t-1)]$, nor $S_i[x..(x+t)]$ and $S_j[y..(y+t)]$ are a $k$-mismatch common substring pair.

In this paper, we present efficient solutions for this problem in both sequential as well as parallel settings. Our sequential algorithm runs in $O(N{log}^{k}N+\mathsf{occ})$ expected time, where $\mathsf{occ}$ is the output size. Our distributed memory parallel algorithm runs in $O\left(((N∕p)logN+\mathsf{occ}){log}^{k}N\right)$ expected time using $O\left({log}^{k+1}N\right)$ expected communication rounds, where $p$ is the number of processors. Here we make a reasonable assumption that the number of occurrence of any $\tau$-long substring across all sequences in $\mathcal{D}$ is $O(N∕p)$. Under this assumption, our algorithm enforces an effective partitioning of a series of modified suffix trees to localize processing within each processor. We demonstrate the scalability and performance of our parallel algorithm using genomic datasets ranging in size from 18 million to over 270 million reads, on up to 1024 processor cores.

To solve such problems in practice, seed-and-extend style filtering approaches are often employed. The underlying principle is: if two sequences have a $k$-mismatch common substring of length $\ge \varphi$, then they must have an exact common substring of length at least $\tau =⌈\frac{\varphi -k}{k+1}⌉$. Therefore, using a fast hashing technique, all pairs of sequences that have a $\tau$-length common substring are identified. Then, by exhaustively checking all such candidate pairs, the final output is generated. In the sequential setting, the filtering heuristics can be broadly classified under three categories: suffix filtering [12], [23], spaced seeds filtering [3], and substring filtering [21]. In case of parallel heuristic methods, filtering based methods mainly focused on the corresponding applications. For example, [11] and [19] proposed suffix tree based parallel clustering of EST data. Clearly, a filtering-based algorithm cannot provide any run time guarantees and often times the candidate pairs generated can be overwhelmingly larger than the final output size. Recently, [22] published the first and only known sub-quadratic exact sequential algorithm for this problem that includes insertions and deletions along with mismatches. Work done on accelerating pairwise edit distance estimations among $n$ sequences using sequence alignment can also be applied to this problem [18]. However, for a large number of short sequences (with low error rate, mostly mismatches), all pair sequence alignment is impractical because of its quadratic time complexity.

This paper is organized as follows. In Section 2, we introduce notations and data structures used in our algorithm. Due to the absence of a provably efficient sequential algorithm for this problem, we first design such an algorithm and present it in Section 3. In Section 4, we describe the parallel algorithm in detail and prove the claimed bounds on expected time and communication rounds. In Section 5, we discuss the implementation details of the parallel algorithm. Finally in Section 6, we discuss the results of our implementation on genomic and gene expression datasets.