Faster algorithms for 1-mappability of a sequence

TL;DR

This paper introduces faster algorithms for the 1-mappability problem in sequences, significantly improving efficiency over previous methods with worst-case and average-case solutions.

Contribution

The authors develop two new algorithms with improved worst-case time complexities and an average-case algorithm for integer alphabets, advancing the state of the art in sequence mappability.

Findings

Algorithms with worst-case time O(mn) and O(n log^2 n)

Average-case algorithm with linear space and time for certain m

Significant reduction in computational complexity for 1-mappability

Abstract

In the k-mappability problem, we are given a string x of length n and integers m and k, and we are asked to count, for each length-m factor y of x, the number of other factors of length m of x that are at Hamming distance at most k from y. We focus here on the version of the problem where k = 1. The fastest known algorithm for k = 1 requires time O(mn log n/ log log n) and space O(n). We present two algorithms that require worst-case time O(mn) and O(n log^2 n), respectively, and space O(n), thus greatly improving the state of the art. Moreover, we present an algorithm that requires average-case time and space O(n) for integer alphabets if m = {\Omega}(log n/ log {\sigma}), where {\sigma} is the alphabet size.

Equations

$$L=\lfloor \frac{m}{3}\rfloor ,B=\lfloor \frac{n}{L}\rfloor$$

$${\text{EXT}}_{i,j}\cdot \text{Occ}$$

$$\mathcal{O}(n+m\cdot \frac{B\cdot n}{{\sigma }^L}).$$

$$m\cdot \frac{B\cdot n}{{\sigma }^L}=\frac{m\cdot \lfloor \frac{n}{\lfloor m/3\rfloor }\rfloor \cdot n}{{\sigma }^{\lfloor \frac{m}{3}\rfloor }}\le \frac{m\cdot (\frac{n}{m/3-1})\cdot n}{{\sigma }^{\frac{m}{3}-1}}\le \frac{12n^2}{n^{\frac{\log \sigma }{\log n}(\frac{m}{3}-1)}}=12n^{2-\frac{(m-3)\log \sigma }{3\log n}}.$$

$$m\ge 3\cdot \frac{\log n}{\log \sigma }+3$$

Full text

1 Introduction

The focus of this work is directly motivated by the well-known and challenging application of genome re-sequencing—the assembly of a genome directed by a reference sequence. New developments in sequencing technologies [15] allow whole-genome sequencing to be turned into a routine procedure, creating sequencing data in massive amounts. Short sequences, known as reads, are produced in huge amounts (tens of gigabytes); and in order to determine the part of the genome from which a read was derived, it must be mapped (aligned) back to some reference sequence that consists of a few gigabases. A wide variety of short-read alignment techniques and tools have been published in the past years to address the challenge of efficiently mapping tens of millions of reads to a genome, focusing on different aspects of the procedure: speed, sensitivity, and accuracy [8]. These tools allow for a small number of errors in the alignment.

The $k$-mappability problem was first introduced in the context of genome analysis in [5] (and in some sense earlier in [3]), where a heuristic algorithm was proposed to approximate the solution. The aim from a biological perspective is to compute the mappability of each region of a genome sequence; i.e. for every factor of a given length of the sequence, we are asked to count how many other times it occurs in the genome with up to a given number of errors. This is particularly useful in the application of genome re-sequencing. By computing the mappability of the reference genome, we can then assemble the genome of an individual with greater confidence by first mapping the segments of the DNA that correspond to regions with low mappability. Interestingly, it has been shown that genome mappability varies greatly between species and gene classes [5].

Formally, we are given a string $x$ of length $n$ and integers $m<n$ and $k<m$, and we are asked to count, for each length-$m$ factor $y$ of $x$, the number of other length-$m$ factors of $x$ that are at Hamming distance at most $k$ from $y$.

Example 1**.**

Consider the string $x=\texttt{aabaaabbbb}$ and $m=3$. The following table shows the $k$-mappability counts for $k=0$ and $k=1$.

[TABLE]

For instance, consider the position 0. The 0-mappability is 1, as the factor aab occurs also at position 4. The 1-mappability at this position is 3 due to the occurrence of aab at position 4 and occurrences of two factors at Hamming distance $1$ from aab: aaa at position $3$ and abb at position 5.

The [math]-mappability problem can be solved in $\mathcal{O}(n)$ time with the well-known LCP data structure [7]. For $k=1$, to the best of our knowledge, the fastest known algorithm is by Manzini [14]. This solution runs in $\mathcal{O}(mn\log n/\log\log n)$ time and $\mathcal{O}(n)$ space and works only for strings over a fixed-sized alphabet. Since the problem for $k=0$ can be solved in $\mathcal{O}(n)$ time, one may focus on counting, for each length-$m$ factor $y$ of $x$, the number of other factors of $x$ that are at Hamming distance exactly $1$ — instead of at most $1$ — from $y$.

Our contributions. Here we make the following threefold contribution:

(a)

We present an algorithm that, given a string of length $n$ over a fixed-sized alphabet and a positive integer $m$, solves the $1$-mappability problem in $\mathcal{O}(\min\{mn,n\log^{2}n\})$ time and $\mathcal{O}(n)$ space, thus improving on the algorithm of [14] that requires $\mathcal{O}(mn\log n/\log\log n)$ time and $\mathcal{O}(n)$ space.

(b)

We present an algorithm to solve the $1$-mappability problem in $\mathcal{O}(mn)$ time and $\mathcal{O}(n)$ space that works for strings over an integer alphabet.

(c)

We present an algorithm that, given a string $x$ of length $n$ over an integer alphabet $\Sigma$ of size $\sigma>1$, with the letters of $x$ being independent and identically distributed random variables, uniformly distributed over $\Sigma$, and a positive integer m=\Omega\Big{(}\frac{\log n}{\log\sigma}\Big{)}, solves the $1$-mappability problem for $x$ in average-case time $\mathcal{O}(n)$ and space $\mathcal{O}(n)$.

The paper is organised as follows. In Section 2, we provide basic definitions and notation as well as a description of the algorithmic tools we use to design our algorithms. In Sections 3 and 4, we provide the average-case and the worst-case algorithms, respectively. We conclude with some final remarks in Section 5.

2 Preliminaries

We begin with some basic definitions and notation. Let $x=x[0]x[1]\ldots x[n-1]$ be a string of length $|x|=n$ over a finite ordered alphabet $\Sigma$ of size $|\Sigma|=\sigma=\mathcal{O}(1)$. We also consider the case of strings over an integer alphabet, where each letter is replaced by its rank in such a way that the resulting string consists of integers in the range $\{1,\ldots,n\}$.

For two positions $i$ and $j$ on $x$, we denote by $x[i\mathinner{.\,.}j]=x[i]\ldots x[j]$ the factor (sometimes called substring) of $x$ that starts at position $i$ and ends at position $j$ (it is of length [math] if $j<i$). By $\varepsilon$ we denote the empty string of length 0. We recall that a prefix of $x$ is a factor that starts at position 0 ($x[0\mathinner{.\,.}j]$) and a suffix of $x$ is a factor that ends at position $n-1$ ($x[i\mathinner{.\,.}n-1]$). We denote the reverse string of $x$ by $\textsf{rev}(x)$, i.e. $\textsf{rev}(x)=x[n-1]x[n-2]\ldots x[1]x[0]$.

Let $y$ be a string of length $m$ with $0<m\leq n$. We say that there exists an occurrence of $y$ in $x$, or, more simply, that $y$ occurs in $x$, when $y$ is a factor of $x$. Every occurrence of $y$ can be characterised by a starting position in $x$. Thus we say that $y$ occurs at the starting position $i$ in $x$ when $y=x[i\mathinner{.\,.}i+m-1]$.

The Hamming distance between two strings $x$ and $y$ of the same length is defined as $d_{H}(x,y)=|\{i:x[i]\neq y[i],\,i=0,1,\ldots,|x|-1\}|$. If $|x|\neq|y|$, we set $d_{H}(x,y)=\infty$. If two strings $x$ and $y$ are at Hamming distance $k$, we write $x\approx_{k}y$.

The computational problem in scope can be formally stated as follows.

1-mappability

Input: A string $x$ of length $n$ and an integer $m$, $1\leq m<n$

Output: An integer array $C$ of size $n-m+1$ such that $C[i]$ stores the number of factors of $x$ that are at Hamming distance $1$ from $x[i\mathinner{.\,.}i+m-1]$

2.1 Suffix array and suffix tree

Let $x$ be a string of length $n>0$. We denote by SA the suffix array of $x$. SA is an integer array of size $n$ storing the starting positions of all (lexicographically) sorted non-empty suffixes of $x$, i.e. for all $1\leq r<n$ we have $x[\textsf{SA}{}[r-1]\mathinner{.\,.}n-1]<x[\textsf{SA}{}[r]\mathinner{.\,.}n-1]$ [13]. Let lcp$(r,s)$ denote the length of the longest common prefix between $x[\textsf{SA}{}[r]\mathinner{.\,.}n-1]$ and $x[\textsf{SA}{}[s]\mathinner{.\,.}n-1]$ for positions $r$, $s$ on $x$. We denote by LCP the longest common prefix array of $x$ defined by LCP$[r]=\textsf{lcp}{}(r-1,r)$ for all $1\leq r<n$, and LCP$[0]=0$. The inverse iSA of the array SA is defined by $\textsf{iSA}{}[\textsf{SA}{}[r]]=r$, for all $0\leq r<n$. It is known that SA, iSA, and LCP of a string of length $n$, over an integer alphabet, can be computed in time and space $\mathcal{O}(n)$ [16, 7]. It is then known that a range minimum query (RMQ) data structure over the LCP array, that can be constructed in $\mathcal{O}(n)$ time and $\mathcal{O}(n)$ space [4], can answer lcp-queries in $\mathcal{O}(1)$ time per query [13]. A symmetric construction on $\textsf{rev}(x)$ can answer the so-called longest common suffix (lcs) queries in the same complexity. The lcp and lcs queries are also known as longest common extension (LCE) queries.

The suffix tree $\mathcal{T}(x)$ of string $x$ is a compact trie representing all suffixes of $x$. The nodes of the trie which become nodes of the suffix tree are called explicit nodes, while the other nodes are called implicit. Each edge of the suffix tree can be viewed as an upward maximal path of implicit nodes starting with an explicit node. Moreover, each node belongs to a unique path of that kind. Thus, each node of the trie can be represented in the suffix tree by the edge it belongs to and an index within the corresponding path. The label of an edge is its first letter. We let $\mathcal{L}(v)$ denote the path-label of a node $v$, i.e., the concatenation of the edge labels along the path from the root to $v$. We say that $v$ is path-labelled $\mathcal{L}(v)$. Additionally, $\mathcal{D}(v)=|\mathcal{L}(v)|$ is used to denote the string-depth of node $v$. Node $v$ is a terminal node if its path-label is a suffix of $x$, that is, $\mathcal{L}(v)=x[i\mathinner{.\,.}n-1]$ for some $0\leq i<n$; here $v$ is also labelled with index $i$. It should be clear that each factor of $x$ is uniquely represented by either an explicit or an implicit node of $\mathcal{T}(x)$. In standard suffix tree implementations, we assume that each node of the suffix tree is able to access its parent. Once $\mathcal{T}(x)$ is constructed, it can be traversed in a depth-first manner to compute $\mathcal{D}(v)$ for each node $v$.

It is known that the suffix tree of a string of length $n$, over an integer alphabet, can be computed in time and space $\mathcal{O}(n)$ [6]. For integer alphabets, in order to access the children of an explicit node by the first letter of their edge label, perfect hashing [9] can be used.

3 Efficient Average-Case Algorithm

In this section we assume that $x$ is a string over an integer alphabet $\Sigma$. Recall that if two strings $y$ and $z$ are at Hamming distance $1$, we write $y\approx_{1}z$.

Fact 1** (Folklore).**

Given two strings $y$ and $z$ of length $m$, we have that if $y\approx_{1}z$, then $y$ and $z$ share at least one factor of length $\lfloor m/2\rfloor$.

Fact 2**.**

Given a string $x$ and any two positions $i,j$ on $x$, we have that if $x[i\mathinner{.\,.}i+m-1]\approx_{1}x[j\mathinner{.\,.}j+m-1]$, then $x[i\mathinner{.\,.}i+m-1]$ and $x[j\mathinner{.\,.}j+m-1]$ have at least one common factor of length $L=\lfloor m/3\rfloor$ starting at positions $i^{\prime}\in\{i,\ldots,i+m-L\}$ and $j^{\prime}\in\{j,\ldots,j+m-L\}$ of $x$, such that $i^{\prime}-i=j^{\prime}-j$ and $i^{\prime}=0\pmod{L}$.

Proof.

It should be clear that every factor of $x$ of length $m$ fully contains at least two factors of length $L$ starting at positions equal to [math] mod $L$. Then, if $x[i\mathinner{.\,.}i+m-1]$ and $x[j\mathinner{.\,.}j+m-1]$ are at Hamming distance $1$, analogously to Fact 1, at least one of the two factors of length $L$ that are fully contained in $x[i\mathinner{.\,.}i+m-1]$ occurs at a corresponding position in $x[j\mathinner{.\,.}j+m-1]$; otherwise we would have a Hamming distance greater than 1. ∎

We first initialize an array $C$ of size $n-m+1$, with [math] in all positions; for all $i$, $C[i]$ will eventually store the number of factors of $x$ that are at Hamming distance $1$ from $x[i\mathinner{.\,.}i+m-1]$. We apply Fact 2 by implicitly splitting the string $x$ into $B={\lfloor\tfrac{n}{\lfloor m/3\rfloor}\rfloor}$ blocks of length $L=\lfloor m/3\rfloor$—the suffix of length $n\bmod\lfloor m/3\rfloor$ is not taken as a block—starting at the positions of $x$ that are equal to [math] mod $L$. In order to find all pairs of length-$m$ factors that are at Hamming distance $1$ from each other, we can find all the exact matches of every block and try to extend each of them to the left and to the right, allowing at most one mismatch. However, we need to tackle some technical details to correctly update our counters and avoid double counting.

We start by constructing the SA and LCP arrays for $x$ and $\textsf{rev}(x)$ in $\mathcal{O}(n)$ time. We also construct RMQ data structures over the LCP arrays for answering LCE queries in constant time per query. By exploiting the LCP array information, we can then find in $\mathcal{O}(n)$ time all maximal sets of indices such that the longest common prefix between any two of the suffixes starting at these indices is at least $L$ and at least one of them is the starting position of some block.

Then for each such set, denoted by $P$, we have to do the following procedure for each index $i\in P$ such that $i=0\pmod{L}$.

For every other $j\in P$, we try to extend the match by asking two LCE queries in each direction. I.e., we ask an $\textsf{lcs}(i-1,j-1)$ query to find the first mismatch positions $\ell_{1}$ and $\ell^{\prime}_{1}$, respectively, and then $\textsf{lcs}(\ell_{1}-1,\ell^{\prime}_{1}-1)$ to find the second mismatch ($\ell_{2}$ and $\ell^{\prime}_{2}$, respectively). A symmetric procedure computes the mismatches $r_{1},r^{\prime}_{1}$ and $r_{2},r^{\prime}_{2}$ to the right, as shown in Figure 1. We omit here some technical details with regards to reaching the start or end of $x$.

Now we are interested in positions $p$ such that $\ell_{2}<p\leq\ell_{1}$ and $i+L-1\leq p+m-1<r_{1}$ and positions $q$ such that $\ell_{1}<q\leq i$ and $r_{1}\leq q+m-1<r_{2}$. Each such position $p$ (resp. $q$) implies that $x[p\mathinner{.\,.}p+m-1]\approx_{1}x[p^{\prime}\mathinner{.\,.}p^{\prime}+m-1]$, where $p^{\prime}=j-(i-p)$. Henceforth, we only consider positions of the type $p,p^{\prime}$.

Note that if $x[p\mathinner{.\,.}p+m-1]\approx_{1}x[p^{\prime}\mathinner{.\,.}p^{\prime}+m-1]$, we will identify the unordered pair $\{p,p^{\prime}\}$ based on the described approach $t_{p,p^{\prime}}$ times, where $t_{p,p^{\prime}}$ is the total number of full blocks contained in $x[p\mathinner{.\,.}p+m-1]$ and in $x[p^{\prime}\mathinner{.\,.}p^{\prime}+m-1]$ after the mismatch position. It is not hard to compute the number $t_{p,p^{\prime}}$ in $\mathcal{O}(1)$ time based on the starting positions $p$ and $p^{\prime}$ as well as $\ell_{1}$ and $r_{1}$ each time we identify $x[p\mathinner{.\,.}p+m-1]\approx_{1}x[p^{\prime}\mathinner{.\,.}p^{\prime}+m-1]$. To avoid double counting, we then increment the $C[p]$ and $C[p^{\prime}]$ counters by $1/t_{p,p^{\prime}}$.

By $\textsf{EXT}_{i,j}$ we denote the time required to process a pair of elements $i,j$ of a set $P$ such that at least one of them, $i$ or $j$, equals [math] mod $L$.

Lemma 1**.**

The time $\textsf{EXT}_{i,j}$ is $\mathcal{O}(m)$.

Proof.

Given $i,j\in P$, with at least one of them equal to [math] mod $L$, we can find the pairs $(p,p^{\prime})$ of positions that satisfy the inequalities discussed above in $\mathcal{O}(m)$ time. They are a subset of $\{(i-m+L,j-m+L),\ldots,(i-1,j-1)\}$. For each such pair $(p,p^{\prime})$ we can compute $t_{p,p^{\prime}}$ and increment $C[p]$ and $C[p^{\prime}]$ accordingly in $\mathcal{O}(1)$ time. The total time to process all pairs $(p,p^{\prime})$ for given $i,j$ is thus $\mathcal{O}(m)$. ∎

Theorem 1**.**

Given a string $x$ of length $n$ over an integer alphabet $\Sigma$ of size $\sigma>1$ with the letters of $x$ being independent and identically distributed random variables, uniformly distributed over $\Sigma$, the 1-mappability problem can be solved in average-case time $\mathcal{O}(n)$ and space $\mathcal{O}(n)$ if $m\geq 3\cdot\frac{\log n}{\log\sigma}+3$.

Proof.

The time and space required for constructing SA and LCP tables for $x$ and $\textsf{rev}(x)$ and the RMQ data structures over the LCP tables is $\mathcal{O}(n)$.

Let $B$ denote the number of blocks over $x$, and let $L$ denote the block length. We set

[TABLE]

to apply Fact 2. Recall that by $P$ we denote a maximal set of indices of the LCP table such that the length of the longest common prefix between any two suffixes starting at these indices is at least $L$ and at least one of them is the starting position of some block. Processing all such sets $P$ requires time

[TABLE]

where $\textsf{EXT}_{i,j}$ is the time required to process a pair $i,j$ of elements of a set $P$; and Occ is the sum of the multiples of the cardinality of each set $P$ times the number of the elements of set $P$ that are equal to [math] mod $L$. By Lemma 1 we have that $\textsf{EXT}_{i,j}=\mathcal{O}(m)$. Additionally, by the stated assumption on the string $x$, the expected value for Occ is no more than $\frac{Bn}{\sigma^{L}}$. Hence, the algorithm on average requires time

[TABLE]

Assuming that $m>3$, we have the following:

[TABLE]

Consequently, in the case when

[TABLE]

the algorithm requires $\mathcal{O}(n)$ time on average. The extra space usage is $\mathcal{O}(n)$. ∎

4 Efficient Worst-Case Algorithms

4.1 $\mathcal{O}(mn)$-time and $\mathcal{O}(n)$-space algorithm

In this section we assume that $x$ is a string over an integer alphabet $\Sigma$. The main idea is that we want to first find all pairs $x[i_{1}\mathinner{.\,.}i_{1}+m-1]\approx_{1}x[i_{2}\mathinner{.\,.}i_{2}+m-1]$ that have a mismatch in the first position, then in the second, and so on.

Let us fix $0\leq j<m$. In order to identify the pairs $x[i_{1}\mathinner{.\,.}i_{1}+m-1]\approx_{1}x[i_{2}\mathinner{.\,.}i_{2}+m-1]$ with $x[i_{1}+j]\neq x[i_{2}+j]$ (i.e. with the mismatch in the $j^{th}$ position), we do the following. For every $i=0,1,\ldots,n-m$, we find the explicit or implicit node $u_{i,j}$ in $\mathcal{T}(x)$ that represents $x[i\mathinner{.\,.}i+j-1]$ and the node $v_{i,j}$ in $\mathcal{T}(\textsf{rev}(x))$ that represents $\textsf{rev}(x[i+j+1\mathinner{.\,.}i+m-1])=\textsf{rev}(x)[n-i-m\mathinner{.\,.}n-i-j-2]$. In each such node $v_{i,j}$, we create a set $V(v_{i,j})$—if it has not already been created—and insert the triple $(u_{i,j},x[i+j],i)$.

When we have done this for all possible starting positions of $x$, we group the triples in each set $V(v)$ by the node variable (i.e., the first component in the triples). For each such group in $V(v)$ we count the number of triples that have each letter of the alphabet and increment array $C$ accordingly. More precisely, if $V(v)$ contains $q$ triples that correspond to the same node $u$, among which $r$ correspond to the letter $c\in\Sigma$, then for each such triple $(u,c,i)\in V(v)$ we increment $C[i]$ by $q-r$; we subtract $r$ to avoid counting equal factors in $C$. Before we proceed with the computations for the next index $j$, we delete all the sets $V(v)$. We formalize this algorithm, denoted by $1$-Map, in the pseudocode presented below and provide an example.

[TABLE]

Example 2**.**

Suppose we have V(v)=\{(u,\texttt{A},i_{1}),(u,\texttt{A},i_{2}),(u,\texttt{A},i_{3}),(u,\texttt{C},i_{4}),(u,\texttt{C},i_{5}),(u,\texttt{C},i_{6}),(u,\texttt{G},i_{7}),$$(u,\texttt{G},i_{8}),(u,\texttt{T},i_{9})\}, for some distinct positions $i_{1},i_{2},$ $\ldots,i_{9}$. We then increment $C[i_{1}]$, $C[i_{2}]$, $C[i_{3}]$, $C[i_{4}]$, $C[i_{5}]$, and $C[i_{6}]$ by $6$; $C[i_{7}]$ and $C[i_{8}]$ by $7$; and $C[i_{9}]$ by $8$.

We now analyze the time complexity of this algorithm. The algorithm iterates $j$ from [math] to $m-1$. In the $j^{th}$ iteration, we need to compute $\{u_{i,j},v_{i,j}\mid i=0,\ldots,n-m\}$. When $j=0$, $u_{i,0}$ for every $i$ is the root of $\mathcal{T}(x)$ and we can find $v_{i,0}$ for all $i$ naïvely in $\mathcal{O}(mn)$ time. For $j>0$, $v_{i,j}$ can be found in $\mathcal{O}(1)$ time from $v_{i,j-1}$ by moving one letter up in ${\cal T}(\mathtt{rev}(x))$ for all $i$, while $u_{i,j}$ can be obtained from $u_{i,j-1}$ by going down in ${\cal T}(x)$ based on letter $x[i+j]$. We then include $(u_{i,j},x[i+j],i)$ in $V(v_{i,j})$.

This requires in total $\mathcal{O}(mn)$ randomized time due to perfect hashing [9] which allows to go down from a node in ${\cal T}(x)$ (or in ${\cal T}(\textsf{rev}(x))$) based on a letter in $\mathcal{O}(1)$ randomized time. We can actually avoid this randomization, as queries for a particular child of a node are asked in our solution in a somewhat off-line fashion: we use them only to compute $v_{i,0}$ ($m$ times) and $u_{i,j}$ (from $u_{i,j-1}$).

Observation 1**.**

For an integer alphabet $\Sigma=\{1,\ldots,n\}$, one can answer off-line $\mathcal{O}(n)$ queries in $\mathcal{T}(x)$ asking for a child of an explicit or implicit node $u$ labelled with the letter $c\in\Sigma$ in (deterministic) $\mathcal{O}(n)$ time.

Proof.

A query for an implicit node $u$ is answered in $\mathcal{O}(1)$ time, as there is only one outgoing edge to check. All the remaining queries can be sorted lexicographically as pairs $(u,c)$ using radix sort. We can also assume that the children of every explicit node of $\mathcal{T}(x)$ are ordered by the letter (otherwise we also radix sort them). Finally, all the queries related to a node $u$ can be answered in one go by iterating through the children list of $u$ once. ∎

Lastly, we use bucket sort to group the triples for each $V(v)$ according to the node variable (recall that the nodes are represented by the edge and the index within the edge) and update the counters in $\mathcal{O}(n)$ time in total (using a global array indexed by the letters from $\Sigma$, which is zeroed in $\mathcal{O}(|V(v)|)$ time after each $V(v)$ has been processed). Overall the algorithm requires $\mathcal{O}(mn)$ time.

The suffix trees require $\mathcal{O}(n)$ space and we delete the sets $V(v_{i,j})$ after the $j^{th}$ iteration; the space complexity of the algorithm is thus $\mathcal{O}(n)$. We obtain the following result.

Theorem 2**.**

Given a string $x$ of length $n$ over an integer alphabet and a positive integer $m$, we can solve the 1-mappability problem in $\mathcal{O}(mn)$ time and $\mathcal{O}(n)$ space.

4.2 $\mathcal{O}(n\log^{2}n)$-time and $\mathcal{O}(n)$-space algorithm

In this section we assume that $x$ is a length-$n$ string over an ordered alphabet $\Sigma$, where $|\Sigma|=\sigma=\mathcal{O}(1)$. Consider two factors of $x$ represented by nodes $u$ and $v$ in $\mathcal{T}(x)$; the first observation we make is that the first mismatch between the two factors is the first letter of the labels of the distinct outgoing edges from the lowest common ancestor of $u$ and $v$ that lie on the paths from the root to $u$ and $v$. For $1$-mappability we require that what follows this mismatch is an exact match.

Definition 1**.**

Let $T$ be a rooted tree. For each non-leaf node $u$ of $T$, the heavy edge $(u,v)$ is an edge for which the subtree rooted at $v$ has the maximal number of leaves (in case of several such subtrees, we fix one of them). The heavy path of a node $v$ is a maximal path of heavy edges that passes through $v$ (it may contain 0 edges). The heavy path of $T$ is the heavy path of the root of $T$.

Consider the suffix tree $\mathcal{T}(x)$ and its node $u$. We say that an (explicit or implicit) node $v$ is a level ancestor of $u$ at string-depth $\ell$ if $\mathcal{D}(v)=\ell$ and $\mathcal{L}(v)$ is a prefix of $\mathcal{L}(u)$. The heavy paths of $\mathcal{T}(x)$ can be used to compute level ancestors of nodes in $\mathcal{O}(\log n)$ time. However, a more efficient data structure is known.

Lemma 2** ([2]).**

After $\mathcal{O}(n)$-time preprocessing on $\mathcal{T}(x)$, level ancestor queries of nodes of $\mathcal{T}(x)$ can be answered in $\mathcal{O}(\log\log n)$ time per query.

Definition 2**.**

Given a string $x$ and a factor $y$ of $x$, we denote by $\textit{range}(x,y)$ the range in the SA of $x$ that represents the suffixes of $x$ that have $y$ as a prefix.

Every node $u$ in $\mathcal{T}(x)$ corresponds to an SA range $I_{u}=\textit{range}(x,\mathcal{L}(u))=(u_{\min},u_{\max})$. We can precompute $I_{u}$ for all explicit nodes $u$ in $\mathcal{T}(x)$ in $\mathcal{O}(n)$ time while performing a depth-first traversal of the tree as follows. For a non-terminal node $v$ with children $u^{1},\ldots,u^{q}$, we set $v_{\min}=\min_{i}\{u^{i}_{\min}\}$ and $v_{\max}=\max_{i}\{u^{i}_{\max}\}$. If $v$ is a terminal node (with children $u^{1},\ldots,u^{q}$), representing the suffix $x[j\mathinner{.\,.}n-1]$, we set $v_{\min}=\textsf{iSA}[j]$ and $v_{\max}=\max\{\textsf{iSA}[j],\max_{i}\{u^{i}_{\max}\}\}$. When a considered node $v$ is implicit, say along an edge $(p,q)$, then $I_{v}=I_{q}$.

Our algorithm relies heavily on the following auxiliary lemmas.

Lemma 3**.**

Consider a node $u$ in $\mathcal{T}(x)$ with $p=\mathcal{L}(u)$. Let $\textit{suf}(u,\ell)$ be the node $v$ such that $\mathcal{L}(v)=p[\ell\mathinner{.\,.}|p|-1]$. Given the SA and the iSA of $x$, $v$ can be computed in $\mathcal{O}(\log\log n)$ time after $\mathcal{O}(n)$-time preprocessing.

Proof.

The SA range of the node $u$ is $I_{u}=(u_{\min},u_{\max})$; $u_{\min}$ corresponds to the suffix $x[\textsf{SA}[u_{\min}]\mathinner{.\,.}n-1]$. By removing the first $\ell$ letters, the suffix becomes $x[\textsf{SA}[u_{\min}]+\ell\mathinner{.\,.}n-1]$. The corresponding SA value is $v_{\min}=\textsf{iSA}[\textsf{SA}[u_{\min}]+\ell]$.

Let $v_{1}$ be the node of $\mathcal{T}(x)$ such that ${\cal L}(v_{1})=x[\textsf{SA}[v_{\min}]\mathinner{.\,.}n-1]$. The sought node $v$ is the ancestor of $v_{1}$ located at string-depth $|p|-\ell$. It can be computed in $\mathcal{O}(\log\log n)$ time using the level ancestor data structure of Lemma 2. ∎

Lemma 4**.**

Let $u$ and $v$ be two nodes in $\mathcal{T}(x)$. We denote ${\cal L}(u)$ by $p_{1}$ and ${\cal L}(v)$ by $p_{2}$. We further denote by $\text{concat}(u,v)$ the node $w$ such that ${\cal L}(w)=p_{1}p_{2}$. Given the SA and the iSA of $x$, as well as $\textit{range}(x,p_{1})$ and $\textit{range}(x,p_{2})$, $w$ can be located in $\mathcal{O}(\log n)$ time after $\mathcal{O}(n)$-time preprocessing.

Proof.

We can compute $\textit{range}(x,p_{1}p_{2})=(w_{\min},w_{\max})$ in $\mathcal{O}(\log n)$ time using the SA and the iSA of $x$ [10, 11]; we can then locate $w$ in $\mathcal{O}(\log\log n)$ time using the level ancestor data structure of Lemma 2. ∎

We are now ready to present an algorithm for $1$-mappability that requires $\mathcal{O}(n\log^{2}n)$ time and $\mathcal{O}(n)$ space. The first step is to build $\mathcal{T}(x)$. We then make every node $u$ of string-depth $m$ explicit in $\mathcal{T}(x)$ and initialize a counter $\textit{Count}(u)$ for it. For each explicit node $u$ in $\mathcal{T}(x)$, the SA range $I_{u}=\textit{range}(x,{\cal L}(u))$ is also stored. We also identify the node $v_{c}$ with path-label $c$ for each $c\in\Sigma$ in $\mathcal{O}(\sigma)=\mathcal{O}(1)$ time.

[TABLE]

We then call PerformCount$(\mathcal{T}(x),m)$, which does the following (inspect also the pseudocode above and Figure 2). At first, a heavy path HP of $\mathcal{T}(x)$ is computed. Initially, we want to identify the pairs of factors of $x$ of length $m$ at Hamming distance $1$ that have a mismatch in the labels of the edges outgoing from a node in HP. Given a node $u$ in HP, with $\mathcal{L}(u)=p_{1}$, for every side tree $S_{i}$ attached to it (say by an edge with label $c\in\Sigma$), we find all nodes of $S_{i}$ with string-depth $m$. For every such node $z$, with path-label $p_{1}cp_{2}$, we use Lemma 3 to obtain the node $w=\text{suf}(z,|p_{1}|+1)$; that is, ${\cal L}(w)=p_{2}$. We then use Lemma 4 to compute $\textit{range}(x,p_{1}c^{\prime}p_{2})$ for all $c^{\prime}\neq c$ such that there is an outgoing edge from $u$ with label $c^{\prime}$ and increment $\textit{Count}(z)$ by $|\textit{range}(p_{1}c^{\prime}p_{2})|$. Let the heavy edge from $u$ have label $d$; we also increment $\textit{Count}(z^{\prime})$, where $z^{\prime}=\text{concat}(u,\text{concat}(v_{d},w))$ is the node with path-label $p_{1}dp_{2}$, by $|I_{z}|$ while processing node $z$.

This procedure then recurs on each of the side trees; i.e. for side tree $S_{i}$, attached to node $u$, it calls PerformCount$(S_{i},m-\mathcal{D}(u))$. Finally, we construct array $C$ from array Count while performing one more depth-first traversal.

On the recursive calls of PerformCount in each of the side trees (e.g. $S_{i}$) attached to HP, we first compute the heavy paths (in $\mathcal{O}(|S_{i}|)$ time for $S_{i}$) and then consider each node of string-depth $m$ of $\mathcal{T}(x)$ at most once; as above, we process each node in $\mathcal{O}(\log n)$ time due to Lemmas 3 and 4. As there are at most $n$ nodes of string-depth $m$, we do $\mathcal{O}(n\log n)$ work in total. This is also the case as we go deeper in the tree. Since the number of leaves of the trees we are dealing with at least halves in each iteration, there at most $\mathcal{O}(\log n)$ steps. Hence, each node of string-depth $m$ will be considered $\mathcal{O}(\log n)$ times and every time we will do $\mathcal{O}(\log n)$ work for it. The overall time complexity of the algorithm is thus $\mathcal{O}(n\log^{2}n)$. The space complexity is clearly $\mathcal{O}(n)$. By applying Theorem 2 we obtain the following result.

Theorem 3**.**

Given a string $x$ of length $n$ over a fixed-sized alphabet and a positive integer $m$, we can solve the 1-mappability problem in $\mathcal{O}(\min\{mn,n\log^{2}n\})$ time and $\mathcal{O}(n)$ space.

5 Final Remarks

The natural next aim is to either extend the presented solutions to work for arbitrary $k$ without increasing the time and space complexities dramatically or develop fundamentally new algorithms if this is not possible. In fact, we already know that the fast average-case algorithm presented in Section 3 can be generalized to work for arbitrary $k$ in linear time. This adds, however, a multiplicative factor of $k$ on the condition for the value of $m$. An interesting generalization of this problem would be to consider the edit distance model instead of the Hamming distance model; i.e. apart from mismatches to also allow for letter insertions and deletions.

Furthermore, a practical extension of the aforementioned problem is the following. Given reads from a particular sequencing machine, the basic strategy for genome re-sequencing is to map a seed of each read in the genome and then try and extend this match [1]. In practice, a seed could be for example the first $32$ letters of the read—the accuracy is higher in the prefix of the read. It is reasonable to allow for a few (e.g. $k=2$) errors when matching the seed to the reference genome to account for sequencing errors and genetic variation. Hence a closely-related problem to genome mappability that arises naturally from this application is the following: What is the minimal value of $m$ that forces $\alpha\%$ of starting positions in the reference genome to have $k$-mappability equal to [math]?

A standard implementation of the algorithm presented in Section 3, when applied to a large sequence of length $n$, requires more than $20n$ bytes of internal memory. Such memory requirements are a significant hurdle to the mappability computation for large datasets on standard workstations. Another direction of practical interest is thus to devise efficient algorithms for the problems of $1$-mappability and $k$-mappability for the External Memory model of computation. Efficient algorithms for computing the suffix array and the longest common prefix array in this model are already known and shown to perform well in practical terms (see [12], for example). Since the average-case algorithm in Section 3 scans the longest common prefix array from left to right sequentially, it would be interesting to see whether it can be implemented efficiently in external memory.