Computing $k$-mers in Graphs

In recent years, a sizeable body of literature in the field of computational biology has emerged on storing and indexing sets of $k$-mers (strings of length $k$). Typically, these $k$-mers are all the overlapping windows of length $k$ in some genome sequence. In strings, the problem of counting the number of distinct $k$-mers is easily solvable in linear time using a combination of suffix and LCP arrays [18, 16], which are computable in $O(n)$ time for polynomial alphabets [15]. There exist highly optimized practical tools for counting short $k$-mers on strings and encoding them efficiently by exploiting this overlap structure [19].

However, a significant effort is underway to move from linear reference genomes into a more rich graph-based reference genome that encodes all variation in a species. Very little attention has been paid to the problem of counting $k$-mers in more general structures such as trees, DAGs and general graphs, some of which can encode an exponential number of distinct $k$-mers. There is some practical work toward this direction in the field of pangenomics [20], but the theoretical limits and possibilities remain largely unexplored. In this paper, we take the first steps towards mapping this uncharted territory.

First, we show that it is #P-hard to compute the number of distinct $k$-mers in directed labeled graphs, even if we restrict ourselves to connected deterministic graphs without directed cycles. However, we find that in the special case of deterministic Wheeler graphs [10], we can compute the number of $k$-mers in $O(|W|k)$ time, where $|W|$ denotes the total number of vertices and edges of $W$, with a simple dynamic programming approach. Pushing further, we study what it takes to improve the time complexity with respect to $k$. Here, the problem turns out to be significantly more complicated than on strings, but we manage an $O(n^4\log k)$ algorithm by generalizing the technique of prefix doubling [18] to Wheeler graphs. We believe there is room to improve the $O(n^4)$ factor. These results pave the way toward counting $k$-mers in arbitrary graphs without explicit enumeration, by reducing them to Wheeler graphs that preserve the $k$-mer set.

As an application of $k$-mer counting, we present an algorithm that builds the node-centric de Bruijn graph of the $k$-mers of a Wheeler graph $W$ in $O(n_k+|W|k\log ({\max }_{1\le \mathrm{\ell }\le k}{n}_{\mathrm{\ell }})+\sigma \log m)$ time, where $n_{\mathrm{\ell }}$ is the number of $\mathrm{\ell }$-mers in $W$. Crucially, the multiplicative factor of $k$ is on the size of the Wheeler graph, rather than the size of the potentially exponential de Bruijn graph (which contains one node per $k$-mer). Previous algorithms to build de Bruijn graphs from graphs or sets of $k$-mers spend $\mathrm{\Omega }(k)$ time per $k$-mer inserted in the graph, which inevitably incurs a cost of $O(k{n}_k)$. In the bioinformatics literature, we find the GCSA2 construction algorithm [20], which builds the de Bruijn graph of a graph by building all distinct $k$-mers of the graph by prefix doubling. Other de Bruijn graph construction algorithms, such as dynamic versions of the data structure by Bowe et al. [5, 4], do not take a graph as input, and only allow the insertion of an explicit $k$-mer at a time, at a worst-case cost of $\mathrm{\Omega }(k)$ per insertion. This application showcases how our theoretical results translate into practical gains for graph algorithms.

The remainder of this paper is organized as follows. In Section 2, we lay down some basic notation and concepts. In Section 3, we discuss the computational complexity of the problem, and prove the problem to be $\mathrm{\#}$P-hard. In Section 4, we solve the problem on Wheeler graphs in time $O(|W|k)$ or $O(n^4\log k)$. In Section 5, we discuss ways of turning an arbitrary graph into a Wheeler graph with the same $k$-mers in order to apply the algorithms in Section 4. Section 6 describes the de Bruijn graph construction algorithm. Section 7 concludes the paper with some closing remarks.

We consider directed graphs with single character labels on the edges. Formally, we consider graphs $G$ formed by a vertex set $V(G)$ and an edge set $E(G)\subseteq V\times V\times \mathrm{\Sigma }$, where $\mathrm{\Sigma }$ is an alphabet of symbols that is totally ordered by a relation $\prec$. The indegree and outdegree of a node $v$ are denoted with $\mathrm{𝗂𝗇𝖽𝖾𝗀}(v)$ and $\mathrm{𝗈𝗎𝗍𝖽𝖾𝗀}(v)$, respectively. A $k$-mer $\alpha$ is a string of length $k$ from ${\mathrm{\Sigma }}^k$, indexed from $1$ to $k$. We compare strings by extending $\prec$ to ${\mathrm{\Sigma }}^{\ast }$ co-lexicographically, which is the right-to-left lexicographical order (we abbreviate this by colex order). We say that a labeled graph $G$ contains a $k$-mer $\alpha$ if there exists a walk $v_0,v_1,\mathrm{\dots },v_k$ in $G$ such that $(v_{i-1},v_i,\alpha [i])\in E(G)$ for all $1\le i\le k$. In this case, we say that $\alpha$ reaches $v_k$, and write $\alpha \rightsquigarrow v_k$. The minimum and maximum $k$-mer that reach a vertex $v$ are denoted ${\min }_v^k$ and ${\max }_v^k$, respectively. A Wheeler graph is a labeled directed graph with the following properties [10]:

We say that a graph is deterministic if no two outgoing edges from the same vertex have the same label. In deterministic Wheeler graphs, we have the following stronger property:

There exists an algorithm to transform a non-deterministic Wheeler graph into a deterministic one with the same $k$-mers, that in the worst case is only a linear factor larger [2]. In contrast, in general graphs the deterministic version could be exponentially larger. However, the running time of the algorithm is $O(n^3)$ [2], and improving this upper bound is an open problem. Thus, we assume that our input are deterministic graphs throughout the paper.

In a labeled graph we can identify the incoming infimum and supremum strings to every vertex, denoted ${inf}_v$ and ${sup}_v$ for all $v\in V$. We use the definition of [3], though note that they read the strings backwards, and use lexicographical order, while we read the strings forward and use colex order instead. Let ${\gamma }_1,{\gamma }_2,\mathrm{\dots },{\gamma }_{2n}$ be the colexicographically sorted infimum and supremum strings of the $n$ vertices of a Wheeler graph. The longest common suffix array is defined as an array $\mathrm{𝖫𝖢𝖲}[\mathrm{2..2}n]$, where $\mathrm{𝖫𝖢𝖲}[i]$ is the longest common suffix of ${\gamma }_{i-1}$ and ${\gamma }_i$. It was first defined in [7], where it was also proven that in a deterministic Wheeler graph, ${inf}_v⪯{sup}_v$, and ${sup}_u⪯{inf}_v$ for all in Wheeler order. Observe that we can find the $\mathrm{𝖫𝖢𝖲}$ value between any infima and suprema in $O(1)$ time, by building a range minimum query data structure [9, 13] over the $\mathrm{𝖫𝖢𝖲}$ array, and performing the appropriate query.

We find it convenient to divide the $\mathrm{𝖫𝖢𝖲}$ array into two parts: one is the external $\mathrm{𝖫𝖢𝖲}$ array, written $\mathrm{𝖤𝖫𝖢𝖲}$, that holds the $\mathrm{𝖫𝖢𝖲}$ between the supremum and infimum string of $v_{i-1}$ and $v_i$. The other one is the internal $\mathrm{𝖫𝖢𝖲}$ array, written $\mathrm{𝖨𝖫𝖢𝖲}$, that holds the $\mathrm{𝖫𝖢𝖲}$ value of the infimum and supremum of the same vertex. We compute them as $\mathrm{𝖤𝖫𝖢𝖲}(v_{i-1},v_i)=\mathrm{𝖫𝖢𝖲}[2i-1]$, and $\mathrm{𝖨𝖫𝖢𝖲}(v_i)=\mathrm{𝖫𝖢𝖲}[2i]$. In Section 4 we show how to use these arrays in the context of $k$-mers.

In this section we show that the $k$-mer counting problem on graphs is $\mathrm{\#}$P-hard even when the graph is acyclic, connected (i.e., there exists a node from which every other node can be reached) and deterministic (i.e., no two edges with the same label have the same origin).

Let $\mathrm{\#}SAT(F)$ denote the number of satisfying assignments to a boolean formula $F$ over $n$ variables $x_1,\mathrm{\dots },x_n$ in conjunctive normal form. Computing $\mathrm{\#}SAT(F)$ is the archetypical $\mathrm{\#}P$-complete problem [21]. It is also well known that the problem is $\mathrm{\#}P$-complete in disjunctive normal form since we can write the number of solutions to a CNF formula $F$ as $2^n-\mathrm{\#}SAT(\neg F)$, where $\neg F$ can be written in disjunctive normal form of the same size as $F$ by de Morgan’s laws.

In this section, we develop two polynomial-time algorithms for counting $k$-mers on deterministic Wheeler graphs. We abbreviate those by DWGs. Let $W=(V,E)$ be a DWG, with $|V|=n$, $|E|=m$ and $|W|=n+m$. The first algorithm has time complexity $O(|W|k)$. For the second algorithm, we reduce the dependency on $k$, but we end up with a time complexity of $O(n^4\log k)$.

These algorithms require, as an internal step, to obtain the LCS array up to a length of $k$, i.e. $\min \{\mathrm{𝖫𝖢𝖲}({\gamma }_i,{\gamma }_{i+1}),k\}$ for all values in the array. The current known LCS construction algorithms take $O(n\log \sigma +\min \{m\log n,m+n^2\})$ time on general deterministic graphs [3], and $O(n\log \sigma +m)$ if the graph is a Wheeler semi-DFA, which requires to have exactly one source. We do not make the assumption of a single source, thus we can only use the former algorithm, which is within the time bound for our $O(n^4\log k)$ counting algorithm. For our $O(|W|k)$ counting algorithm, we present a tailored LCS construction that is constructed during the $k$-mer counting algorithm at no extra cost, but with values capped to $k$, which is enough for our purposes.

The algorithms in Section 4 are of limited usability since they apply only to deterministic Wheeler graphs. In this Section we discuss how to transform some more general graphs into deterministic Wheeler graphs.

An layered DAG with $k$ layers is a directed acyclic graph that can be partitioned into $k$ disjoint subsets of vertices $V_0,\mathrm{\dots }{V}_{k-1}$ such that edges exists only from vertices in $V_i$ into $V_{i+1}$ for $i=0,\mathrm{\dots },k-2$. For any graph, we can build the $k$ times unfolded DAG, which preserves the $k$-mers in the graph:

We can turn this into a DFA by introducing a new initial state with epsilon-transitions into the first layer, and determinizing with the classic subset construction. Now, since every $k$-mer corresponds to a distinct path from the initial state to the final layer, we can count $k$-mers by counting the distinct paths to each node on the final layer of the DAG, which is easily done with a dynamic programming algorithm that processes the nodes in topological order.

However, this does not exploit overlaps between $k$-mers. For deterministic acyclic graphs that have a source node that can reach all nodes, we may use the Wheelerization algorithm of Alanko et al. [2] that computes the smallest equivalent Wheeler DFA in $O(n\log n)$ time in the size of the output. More generally, if the input is the state graph $G$ of any DFA recognizing a Wheeler language, the minimum equivalent Wheeler DFA can be constructed in $O({|G|}^{3}n\log n)$ time, where $n$ is the size of the output [8].

Another option is to build the de Bruijn graph of the $k$-mers, which is always Wheeler [10]. Unfortunately, this does not help much with counting the $k$-mers because building the de Bruijn graph takes time proportional to the number of $k$-mers, since there is one node per $k$-mer. That is, building the de Bruijn graph amounts to a brute force enumeration of all $k$-mers. However, the de Bruijn graph is not necessarily the smallest graph for a given set of $k$-mers: For example, the Wheeler DFA minimization algorithm of Alanko et al. [1], run on the de Bruijn graph, can reduce the size of the graph while maintaining Wheelerness. Building this compressed graph directly without first building the de Bruijn graph would be useful, but seems to be nontrivial. We leave this open for future work.

In this section, we define a representation of the node-centric de Bruijn graph of the $k$-mers of a DWG $W$, building on the definitions of black and white $k$-mers and the counting algorithms introduced in the previous section. We denote the de Bruijn graph by $dBg$. First, we describe the data structure contents, then we describe the algorithms to navigate the de Bruijn graph given this data structure, and then we analyse their time and space in bits. Finally we show how to build it in time and space $O(n_k+|W|k\log ({\max }_{1\le \mathrm{\ell }\le k}{n}_{\mathrm{\ell }})+\sigma \log m)$. The main technical difficulty of this section is that each vertex $v\in V(W)$ may have exponentially many white $k$-mers that reach only $v$. Thus, to support de Bruijn graph navigation we need to keep track of the rank of the $k$-mer among all the $k$-mers reaching $v$, and update it to the rank in another vertex as we traverse an edge.