A practical fpt algorithm for Flow Decomposition and transcript assembly

TL;DR

This paper introduces a practical fixed-parameter tractable algorithm for Flow Decomposition in DAGs, crucial for transcript assembly, and demonstrates its efficiency and exactness on RNA-sequence data.

Contribution

The paper presents the first practical linear FPT algorithm for Flow Decomposition, along with an implementation that outperforms heuristics in transcript assembly tasks.

Findings

The solver finds exact solutions efficiently on RNA-sequence data.

The algorithm's runtime is competitive with state-of-the-art heuristics.

Hardness results show no polynomial kernels and NP-hardness for related problems.

Abstract

The Flow Decomposition problem, which asks for the smallest set of weighted paths that "covers" a flow on a DAG, has recently been used as an important computational step in transcript assembly. We prove the problem is in FPT when parameterized by the number of paths by giving a practical linear fpt algorithm. Further, we implement and engineer a Flow Decomposition solver based on this algorithm, and evaluate its performance on RNA-sequence data. Crucially, our solver finds exact solutions while achieving runtimes competitive with a state-of-the-art heuristic. Finally, we contextualize our design choices with two hardness results related to preprocessing and weight recovery. Specifically, $k$-Flow Decomposition does not admit polynomial kernels under standard complexity assumptions, and the related problem of assigning (known) weights to a given set of paths is NP-hard.

Figures

Tables

Table 1. Table 1 . Summary of the RNA sequencing dataset  [ 15 ] . Statistics are for non-trivial instances. Columns 4 through 6 give averages; column 7 (8) reports the percent of non-trivial ground truth decompositions that are optimal (non-optimal) size. Because Toboggan timed out on some instances, these percentages do not sum to 100. The high percentage of instances with ground truth of minimum size supports the use of Flow Decomposition as a model for transcript assembly.

dataset	instances	non-trivial	nodes	deg	$k$	optimal	non-optimal
zebrafish	1,549,373	445,880	18.49	2.27	2.32	99.91%	0.053%
mouse	1,316,058	473,185	18.67	2.37	2.75	99.40%	0.074%
human	1,169,083	529,523	18.79	2.41	2.83	99.49%	0.043%
human-salmon	13,300,893	7,153,472	20.54	2.55	3.74	94.39%	0.035%
all	17,335,407	8,602,060	20.22	2.52	3.55	95.27%	0.039%

Equations

$$f(a)=\sum _{i=1}^k{w}_i\cdot 1_{p_i}(a)$$

$$\sum _{i\in g^{\prime -1}(a)}{w}_i=f(a).$$

$$\sqrt{2\pi k}{\left(\frac{k}{e}\right)}^k{e}^{1/(12k+1)}\le k!\le \sqrt{2\pi k}{\left(\frac{k}{e}\right)}^k{e}^{1/12k}.$$

$$\left\{\genfrac{}{}{0ex}{}{k}{\ell }\right\}\ell !$$

$$\le \frac{\sqrt{k}}{2\sqrt{k-\ell }}{\left(\frac{k}{e}\right)}^k{\left(\frac{e}{k-\ell }\right)}^{k-\ell }{e}^{\frac{1}{12k}-\frac{1}{12(k-\ell )+1}}\frac{{\ell }^k}{{\ell }^{\ell }}$$

$$\le \sqrt{k}\frac{k^k{\ell }^{k-l}}{(k-\ell {)}^{k-\ell }{e}^{\ell }}.$$

$$\sqrt{k}\frac{(\alpha k{)}^{k-\alpha k}}{(k-\alpha k{)}^{k-\alpha k}{e}^{\alpha k}}{k}^k\le \sqrt{k}{\left({\left(\frac{\alpha }{1-\alpha }\right)}^{(1-\alpha )}{e}^{-\alpha }\right)}^k{k}^k\le \sqrt{k}(0.649k{)}^k,$$

$$\frac{4^{k^2}}{k!k^k}\cdot \sqrt{k}(0.649k{)}^k=\sqrt{k}\frac{4^{k^2}0.649^k}{k!}.$$

$$\frac{4^{k^2}{k}^{2.5k}0.649^k}{k!}{k}^{o(k)}\cdot (n+\lambda )\le 4^{k^2}{k}^{1.5k}1.765^k{k}^{o(k)}\cdot (n+\lambda )$$

$$k\ge |\mathcal{F}(C_1)\cap \mathcal{F}(C_2)|+\frac{2}{3}(|F_1|+|F_2|),$$

Full text

1. Introduction

We study the problem Flow Decomposition [6, 10, 16, 19] under the paradigm of parameterized complexity [4]. Motivated by the principles of algorithm engineering, we design and implement an fpt algorithm that not only solves the problem exactly but also has a run- time that is competitive with heuristics [16, 19]. Furthermore, we characterize several key aspects of the problem’s complexity. Decomposing flows is the central algorithmic problem in a recent method for analyzing high-throughput transcriptomic sequencing data, and we benchmark our implementation on data from this use case [16].

Flow Decomposition asks for the minimum number of weighted paths necessary to exactly cover a flow on a directed acyclic graph with a unique source $s$ and sink $t$ (an $s$-$t$–DAG). More precisely, we say a set of $s$-$t$–paths $P=p_{1},\dots,p_{k}$ and corresponding weights $\mathbf{w}=(w_{1},\dots,w_{k})$ are a flow decomposition of an $s$-$t$–DAG $G$ with flow $f$ if

[TABLE]

for every arc $a$ in $G$, where $\mathbf{1}_{p_{i}}(a)$ is the indicator function whose output is 1 if $a\in p_{i}$, and [math] otherwise. Specifically, in this paper, we are concerned with $k$-Flow Decomposition, which uses the natural parameter of the number of paths.

Problem 1.1**.**

$k$*-Flow Decomposition ($k$-FD) \Input& ($G,f,k$) with $G$ an $s$-$t$–DAG, $f$ a flow on $G$, and $k$ a positive integer.

\Prob Is there an integral flow decomposition of $(G,f)$ using at most $k$ paths?*

Shao and Kingsford recently used flow decompositions to assemble unknown transcripts from observed subsequences of RNA called exons [16]. In this context, the input graph is a splice graph, which is constructed in one of two ways: either by contracting induced paths in a de Bruijn graph (de novo assembly) or by connecting exons inferred from an alignment process that co-occur in reads (reference-based). In both constructions, every vertex in the graph corresponds to a (partial) exon being transcribed and every arc indicates that two exons appear in sequence in one of the reads; the arcs are weighted based on abundance (frequency) in the data. If we affix a dummy source $s$ and a dummy sink $t$ to the vertices of in- and out-degree zero, respectively, the resulting labeling forms a flow111This assumes an idealized scenario with constant read coverage, but methods exist to rectify noisy data [18]. from $s$ to $t$. As such, a transcript corresponds to a weighted $s$-$t$–path and we aim to recover the collection of transcripts by decomposing the flow into a minimal number of paths in the above sense.

A precursory investigation of the data used by Shao and Kingsford [16] to evaluate the performance of Flow Decomposition heuristics in RNA transcript assembly left us with three guiding observations:

1. (1)

97% of instances admit decompositions into fewer than 8 paths. Our algorithm should thus exploit the small natural parameter. 2. (2)

The data set contains over 17 million mostly small instances. Our implementation should therefore be able to handle a large throughput. 3. (3)

The flow decompositions computed by our implementation should reliably recover the domain-specific solution i.e. the true transcripts.

The first point corroborates the parameterized approach. Using dynamic programming over a suitable graph decomposition, a common algorithm design technique from parameterized complexity, we show the problem can be solved in linear-fpt time (Section 3):

Theorem 1.2**.**

There is an $2^{O(k^{2})}(n+\lambda)$ algorithm for solving $k$-Flow Decomposition, where $\lambda$ is the logarithm of the largest flow value.

To address the second point, we implement a Flow Decomposition solver, Toboggan [7], based on this algorithm and compare it with the state-of-the-art heuristic Catfish (Section 5). Our results show that Toboggan’s running time is comparable to Catfish and thus suitable for high-throughput applications. With respect to the third point, using Flow Decomposition for assembly implicitly assumes that the true transcripts correspond to paths in a minimum-size flow decomposition. Prior work focuses on heuristics which cannot be used to evaluate the validity of this assumption. With Toboggan, we validate that minimum-size flow decompositions accurately recover transcripts in most instances (Section 5.3).

Toboggan incorporates several heuristic improvements (discussed in Section 4) to keep the running time and memory consumption of the core dynamic programming routine as small as possible. In particular, the simple preprocessing we employ is highly successful in solving many instances directly. A provable guarantee for this preprocessing in the form of a small kernel, however, is unlikely: we show that unless $\mathsf{NP}$ is contained in {\mathchoice{\hbox{\mathsf{coNP}}}{\hbox{\mathsf{coNP}}}{\mathsf{coNP}}{\mathsf{coNP}}}/\emph{poly}, $k$-FD does not admit a polynomial kernel (Section 6). Further, our experimental evaluation hints at the following conjecture: for a given decomposition of an $s$-$t$–DAG into a minimum number of paths, there is a unique assignment of weights to those paths to achieve a flow-decomposition. If the conjecture holds, it implies a tighter bound on the running time of the algorithm in Theorem 1.2. Despite this, we prove the more general ‘weight recovery’ problem is $\mathsf{NP}$-hard given a decomposition into an arbitrary number of paths (Section 7).

2. Preliminaries

Related work.

Flow Decomposition is known to be $\mathsf{NP}$-complete, even when all flow values are in $\{1,2,4\}$, and does not admit a PTAS [6, 19]. The best known approximation algorithm for the problem is based on parity-balancing path flows [10], and guarantees an approximation ratio of $\lambda L^{\lambda}$ with a running-time of $O(\lambda|V|\cdot|E|^{2})$, where $L$ is the length of the longest $s$-$t$–path, and $\lambda$ is the logarithm of the maximum flow value. The variant in which the decomposition approximates the original flow values to minimize a specified penalty function was studied in [17], and the authors give an fpt algorithm and FPTAS parameterized by $k$, the largest $s$-$t$-cut, and the maximum flow value.

Due to its practical use for sequencing data, much of the prior work on Flow Decomposition has focused on fast heuristics, including two greedy algorithms which iteratively add the remaining path which is shortest (greedy-length) or of maximum possible weight (greedy-width) [19]. Both heuristics produce solutions with at most $k=|E|-|V|+2$ paths and have variants [6] which decompose all but an $\varepsilon$-fraction of the flow using at most $O(1/\varepsilon^{2})$ times the minimum number of paths, for any $\varepsilon>0$. Historically, greedy-width has provided the best performance [6, 10, 19], but the recent heuristic Catfish [16] showed significant improvements over greedy-width in accuracy and runtime.

Notation.

Given a directed acylic graph (DAG), $G=(V,A)$, we say $G$ is an $s$-$t$–DAG if $G$ has a single source, $s$, and a single sink, $t$. We denote by $A^{+}(v)$ the set of out-arcs of a node $v$, and by $A^{-}(v)$ the in-arcs. For a set of nodes, $S$, we define $A^{+}(S)=\{vu\mid v\in S,u\notin S,vu\in A^{+}(v)\}$. If $f$ is a flow on $G$, we write $f(a)$ for the flow value on an arc $a$ and $F$ for the total flow (the sum of flow values on the arcs in $A^{+}(s)$).

Terminology.

Given a DAG, a topological ordering on the nodes is a labeling such that every arc is directed from a node with a smaller label to one with a larger label. We label the nodes of an $s$-$t$–DAG $G$ as $v_{1},\dots,v_{n}$ corresponding to an arbitrary, fixed topological ordering of $G$; this implies that $s=v_{1}$ and $t=v_{n}$. We further define the sets $S_{i}=\{v_{j}\mid j\leq i\}$ based on the ordering. We refer to the $s$-$t$–cuts $A^{+}(S_{i})$ as topological ordering cuts (top-cuts). Our algorithm for computing a flow decomposition will ultimately depend on tracking how the paths cross top-cuts, a notion we refer to as a routing.

Definition 2.1** (Routings and extensions).**

A surjective function $g\colon[k]\to A^{+}(S_{i})$ is a routing out of $S_{i}$. A routing $g^{\prime}\colon[k]\to A^{+}(S_{i})$ is an extension of $g\colon[k]\to A^{+}(S_{i-1})$ if for each $j\in[k]$, $g^{\prime}(j)=xy$ implies $g(j)=xy$ if $x\neq v_{i}$, and $g(j)=zx$ for some $z\in S_{i-1}$ if $x=v_{i}$.

In other words, an extension of a routing takes all the integers that map to in-arcs of $v_{i}$ and maps them instead to the out-arcs of $v_{i}$ (because it is surjective) while leaving the rest of the mapping unchanged. It will be important in our analysis to differentiate arcs that occur in multiple paths from arcs that appear on only a single path: we say paths $p$ and $p^{\prime}$ coincide on arc $a$ if $a\in p$ and $a\in p^{\prime}$. An arc $a$ is saturated by a path $p$ if $p$ is the only path for which $a\in p$.

Parameterized complexity.

A parameterized decision problem $\Pi\in\Sigma^{*}\times\mathbb{N}$ is fixed parameter tractable if there exists an algorithm that decides it in time $g(k)\cdot n^{c}$ for some computable function $g$ and a constant $c$. When $c=1$, we call such an algorithm linear fpt. A polynomial kernel is an algorithm that transforms, in polynomial time, an instance $(I,k)\in\Sigma^{*}\times\mathbb{N}$ of $\Pi$ into an equivalent instance $(I^{\prime},k^{\prime})$ with $|I^{\prime}|,k^{\prime}\leq k^{O(1)}$, meaning that $(I,k)\in\Pi$ if and only if $(I^{\prime},k^{\prime})\in\Pi$. Some more advanced machinery pertaining to kernelization lower bounds will be defined in Section 6.

3. A linear fpt algorithm for $k$-FD

We solve $k$-FD via dynamic programming over the topological ordering: we enumerate all ways of routing $k$ $s$-$t$–paths over each top-cut $A^{+}(S_{i})$. Each such routing determines a set of constraints for the path weights, which we encode in linear systems. For example, if paths $p_{1}$ and $p_{2}$ are routed over an arc $a$, we add the constraint $w_{1}+w_{2}=f(a)$ to our system.

Concretely, we keep a table $T_{i}$ for $0\leq i\leq n$ whose entries are sets of pairs $(g,L)$, where $g\colon[k]\to A^{+}(S_{i})$ is a routing of the paths out of $S_{i}$ and $L$ is a system of linear equations that encodes the known path weight constraints. In particular, for every arc $a\in A^{+}(S_{i})$, we write $g^{-1}(a)$ for the set of $s$-$t$–paths that are routed over $a$. Each system of linear equations, $L$, is of the form $\mathbf{Aw}=\mathbf{b}$, where $\mathbf{A}$ is a binary matrix with $k$ columns, $\mathbf{w}$ is the solution weight vector, and $\mathbf{b}$ is a vector containing values of $f$. Each row $r$ of $\mathbf{A}$ corresponds to an arc $a_{r}$ and encodes the constraint that the weights of paths routed over $a_{r}$ sum up to $f(a_{r})$. The $j$th entry of row $r$ is 1 if and only if path $p_{j}$ was routed over $a_{r}$, and then $f(a_{r})$ equals the $r$th entry of $\mathbf{b}$. Figure 1 illustrates an example of an entry $(g,L)$.

We now describe how the dynamic programming tables are constructed. For ease of description, we augment $G$ to have two additional dummy arcs $a_{s}$ and $a_{t}$, where $a_{s}$ is an in-arc of $s$ and $a_{t}$ is an out-arc of $t$ and $f(a_{s})=f(a_{t})=F$. We begin with table $T_{0}$ that has a single entry. The routing of this entry routes all paths over $a_{s}$ and the linear system has a single row that constrains all path weights to sum to the total flow value.

For $i>0$, we construct the table $T_{i}$ from $T_{i-1}$. Conceptually, we “visit” node $v_{i}$ and “push” all paths routed over its in-arcs to be routed over its out-arcs. Formally, we require the routing out of $S_{i}$ to be an extension (Definition 2.1) of a routing out of $S_{i-1}$. For each entry $(g,L)\in T_{i-1}$ we compute an entry $(g^{\prime},L^{\prime})$ of $T_{i}$ for each extension $g^{\prime}$ of $g$. For a given $g^{\prime}$, we create $L^{\prime}$ by adding a row to $L$ for each arc $a\in A^{+}(v_{i})$, encoding the constraint

[TABLE]

At the conclusion of the dynamic programming, all entries of the final table $T_{n}$ will have the same routing, since $v_{n}=t$ has one (dummy) out-arc $a_{t}$. This table allows us to decide whether there is a solution; $(G,f)$ is a yes-instance of $k$-FD if and only if some $L\in T_{n}$ has a solution $\mathbf{w}$ whose entries are positive integers. Pseudocode for our algorithm can be found in Algorithm 1 in Appendix A.

Lemma 3.1**.**

An integer vector $\mathbf{w}$ is a solution to $L\in T_{n}$ if and only if there is a set of $s$-$t$–paths $P$ such that $(P,\mathbf{w})$ is a flow decomposition of $(G,f)$.

Proof 3.2**.**

We first prove the forward direction. If we keep backpointers in our DP tables i.e. pointers from entry $(g_{i},L_{i})$ to the entry in the previous table $(g_{i-1},L_{i-1})$ for which $g_{i}$ was the extension of $g_{i-1}$, we can obtain a sequence of routings $g_{n-1},\dots,g_{1}$ that correspond to backwards traversal of the backpointers. Let $h(j)=g_{1}(j),\dots,g_{n-1}(j)$. By Definition 2.1, removing the duplicate elements that appear in consecutive order from $h(j)$ yields a series of arcs that form an $s$-$t$–path $p_{j}$. Because each routing requires every arc to have a path routed over it, the system $L$ contains constraints corresponding to each arc. Thus an integer solution $\mathbf{w}$ to $L$ corresponds to a weighting of $P=p_{1},\dots,p_{k}$ such $(P,\mathbf{w})$ is a flow decomposition.

In the reverse direction, we first observe that the incidence of paths in $P$ on a top-cut $A^{+}(S_{i})$ corresponds to a routing out of $S_{i}$. Let $g_{i}$ be the corresponding routing out of $S_{i}$. For each node $v_{i}$, the paths routed over $A^{-}(v_{i})$ must immediately be routed over an arc in $A^{+}(v_{i})$, meaning $g_{i}$ is an extension of $g_{i-1}$. Because $(P,\mathbf{w})$ is a flow decomposition, the paths routed over each arc will have weights summing to the flow value on that arc. Thus, any constraint in a linear system associated with the routings $g_{i}$ will have $\mathbf{w}$ as a solution.

To analyze the running time we now derive bounds on the dynamic programming table sizes. First, we bound the number of possible routings in Lemma 3.3, then we give an upper bound on the number of linear systems in Lemma 3.5.

Lemma 3.3**.**

There are at most $\sqrt{k}\,\left(0.649k\right)^{k}$ routings of $k$ paths over any top-cut.

Proof 3.4**.**

A routing of $k$ paths over a top-cut $C$ can be thought of as a partition of the paths into $\ell=|C|$ many non-empty sets along with a mapping of these sets to the arcs of $C$. We can assume that $\ell\leq k$, since a routing of $k$ paths requires every top-cut to have at most $k$ arcs. The number of ways to partition $k$ objects into $\ell$ non-empty sets is $\genfrac{\{}{\}}{0.0pt}{}{k}{\ell}$, the Stirling number of the second kind, and there are $\ell!$ ways to assign each of these partitions to a specific arc. Thus, the number of routings is $\genfrac{\{}{\}}{0.0pt}{}{k}{\ell}\ell!$. To proceed, we use the upper bound $\genfrac{\{}{\}}{0.0pt}{}{k}{\ell}\leq\frac{1}{2}{k\choose\ell}\ell^{k-\ell}$ due to Rennie and Dobson [12]. We also make use of the tighter version of Stirling’s approximation due to Robbins [13], which states that

[TABLE]

Hence, we have the upper bound

[TABLE]

Letting $\ell=\alpha k$ for $\alpha\in(0,1)$, the above expression becomes

[TABLE]

where the constant $0.649$ can be derived numerically by maximizing $g(\alpha)=\left(\frac{\alpha}{1-\alpha}\right)^{(1-\alpha)}e^{-\alpha}$ on the interval $[0,1]$.

Lemma 3.5**.**

Each table has at most $\frac{4^{k^{2}}}{k!\,k^{k}}$ distinct linear systems.

Proof 3.6**.**

Without loss of generality, we can ensure each linear system $L$ has at most $k$ rows by removing linearly dependent rows. We note that because there are only $2^{k}$ subsets of weights, if $f$ maps the arcs to more than $2^{k}$ unique flow values, there cannot be a flow decomposition of size $k$. Since the elements of $\mathbf{b}$ can thus take on at most $2^{k}$ many values, and $\mathbf{A}$ contains binary rows of length $k$, it follows that there are at most $4^{k}$ rows for $L$. Accordingly, we can bound the number of possible linear systems by $\binom{4^{k}}{k}$. By imposing an order on the rows we can remove a factor of $1/k!$. Because $k\leq\sqrt{4^{k}}$, we can apply the bound $\binom{n}{k}\leq(n/k)^{k}$ which holds [8] for $k\leq\sqrt{n}$; thus, the number of linear systems is at most $4^{k^{2}}/(k!\,k^{k})$.

Now that we have bounded the size of the dynamic programming tables, we analyze the complexity of solving the linear systems in the final table $T_{n}$. It has been shown that treating linear systems as integer linear programs (ILPs) produces integer solutions in fpt-time.

Proposition 3.7** ([9]).**

Finding an integer solution to a given system $L$ takes time at most $O(k^{2.5k+o(k)}|L|)$, where $|L|$ is the encoding size of the linear system.

These results give the following upper bound on the runtime of our algorithm, proving Theorem 1.2.

Theorem 3.8**.**

Algorithm 1 solves $k$-FD in time $4^{k^{2}}k^{1.5k}k^{o(k)}1.765^{k}\cdot(n+\lambda)$ where $\lambda$ is the logarithm of the largest flow value of the input.

Proof 3.9**.**

The correctness of the algorithm was already proven in Lemma 3.1, so all that remains is to bound the running time. By Lemmas 3.3 and 3.5, the total number of elements in DP table $T_{i}$ is bounded by

[TABLE]

Reducing the linear systems and checking for consistency is polynomial in the size of the matrix (which is $k\times k$). Finally, we need to find a feasible solution among all the linear systems left after the last DP step. We apply Proposition 3.7 to these systems, whose encoding size is bounded by $k^{O(1)}\lambda$. We arrive at the desired running time of

[TABLE]

where we use the well-known bound $\tfrac{k^{k}}{k!}\leq e^{k}$.

4. Implementation

To establish that our exact algorithm for $k$-FD is a viable alternative to the heuristics currently in use by the computational biology community, we implemented Algorithm 1 as the core of a Flow Decomposition solver Toboggan [7]. The solver iterates over values of $k$ in increasing order until reaching a yes-instance of $k$-FD. Toboggan also implements backtracking to recover the $s$-$t$–paths. Making Toboggan’s runtime competitive with existing implementations of state-of-the-art heuristics required non-trivial algorithm engineering. Our improvements broadly fall into three categories: preprocessing, pruning, and low-memory strategies for exploring the search space. The remainder of this section describes the most noteworthy techniques implemented in Toboggan that are not captured by Algorithm 1.

4.1. Preprocessing

Toboggan implements two key preprocessing routines. The first generates an equivalent instance on a simplified graph using a series of arc contractions. The second calculates lower bounds on feasible values of $k$ to reduce the number of calls to the $k$-FD solver.

Graph reduction.

We first reduce the graph by contracting all arcs into/out of nodes of in-/out-degree 1. We prove in Lemma 4.1 that given a flow decomposition of the contracted graph, we can efficiently recover a decomposition of the same size for the original graph.

Lemma 4.1**.**

Let $uv$ be an arc for which $|A^{+}(u)|=1$ or $|A^{-}(v)|=1$ and $G^{\prime}$ the graph created by contracting $uv$. Then $(G^{\prime},f)$ has a flow decomposition $(P^{\prime},\mathbf{w})$ of size $k$ iff $(G,f)$ has a flow decomposition $(P,\mathbf{w})$ of size $k$. Moreover, we can construct $(P,\mathbf{w})$ from $(P^{\prime},\mathbf{w})$ in polynomial time.

Proof 4.2**.**

We first note that if $G^{R}$ is the graph formed by reversing the directions of the arcs in $G$, any solution to $G$ can be transformed into a solution to $G^{R}$ by reversing the order of the paths and maintaining the same weights. Since this reversal is involutive, the correspondence between solutions to $G$ and $G^{R}$ is bijective, so it suffices to consider the case $|A^{+}(u)|=1$.

Given $G$ with a node $u$ that has $|A^{+}(u)|=1$, let $G^{\prime}$ be the graph $G$ with arc $uv$ contracted. Given a flow decomposition $(P,\mathbf{w})$ of $(G,f)$, we construct a corresponding decomposition $(P^{\prime},\mathbf{w})$ for $(G^{\prime},f)$ as follows. Every path $p\in P$ containing $u$ must have $v$ as the vertex succeeding $u$. Removing $u$ from each such $p$ will create a valid path in $G^{\prime}$, since $A^{-}(u)$ becomes part of $A^{-}(v)$ in $G^{\prime}$. Moreover, the incidence of paths on each other arc remains unchanged, so the solution covers the flows in $G^{\prime}$.

In the reverse direction, consider a labeling of the arcs of $G^{\prime}$ such that we can distinguish among the in-arcs of $v$ those which exist in $G$ from those that result from contracting $uv$. Let $A_{\text{new}}$ be the set of latter arcs. We construct $P$ from $P^{\prime}$ as follows. For each path $p\in P^{\prime}$, if $p$ is routed over $xv\in A_{\text{new}}$ we modify $p$ to include the arcs $xu$ and $uv$, rather than the arc $xv$, and add the modified path to $P$. If no such arc lies in $p$, we simply add $p$ to $P$.

As in the forward direction, it is clear that any arcs in $G$ without $u$ as an endpoint have their flow values covered by $P$, and that every path in $P$ is an $s$-$t$–path in $G$. Since $xv\in A_{\text{new}}$ is derived from an arc $xu\in G$, if $Q^{\prime}\subset P^{\prime}$ are the paths routed over $xv$ and $Q$ is the set of paths corresponding to $Q^{\prime}$ in $P$, then $Q$ exactly covers $xu$. Since every path in $P$ routed through $A^{-}(u)$ subsequently is routed over $uv$ and the flow values on the in-arcs of $u$ sum to $f(uv)$, $uv$ is also covered by $P$. Thus $(P,\mathbf{w})$ is a valid solution to $(G,f)$.

Lower bounds on $k$.

To reduce the number of values of $k$ that Toboggan considers before reaching a yes-instance, we compute a lower bound on the optimal value of $k$. One immediate lower bound is the size of the largest top-cut; we implement this in conjunction with the additional bound established in Lemma 4.3. Intuitively, if two top-cuts share few common flow values, many of the arcs must have multiple paths routed over them.

Lemma 4.3**.**

Given a flow $(G,f)$, let $C_{1}$ and $C_{2}$ be any two top-cuts with $|C_{1}|\geq|C_{2}|$. Letting $\mathcal{F}(S)$ be the multiset of flow values occurring in $S$, set $F_{1}=\mathcal{F}(C_{1})\setminus\mathcal{F}(C_{2})$ and $F_{2}=\mathcal{F}(C_{2})\setminus\mathcal{F}(C_{1})$. If $(G,f)$ has a flow decomposition of size $k$, then

[TABLE]

and this lower bound is tighter than the cutset size $|C_{1}|$ iff $|F_{1}|<2|F_{2}|$.

Proof 4.4**.**

Suppose $k$ paths cross cutsets $C_{1}$ and $C_{2}$. For every flow value in $\mathcal{I}=\mathcal{F}(C_{1})\cap\mathcal{F}(C_{2})$, it is possible a single path saturates the arc with that flow value in both cuts. Consider the remaining $h=k-|\mathcal{I}|$ paths. To maximize the number of distinct flow values these $h$ paths produce, let each path saturate an arc in $C_{1}$, yielding $h$ values in $F_{1}$, and then route those $h$ paths in pairs across distinct arcs in $C_{2}$. This produces at most $h/2$ new flow values in $F_{2}$. This yields $(3/2)h\geq|F_{1}|+|F_{2}|$, and substituting for $h$ proves Inequality (1).

To prove the relationship between this lower bound and $|C_{1}|$ there are two cases. If $|F_{1}|\geq 2|F_{2}|$, then by substituting we can upper bound $(2/3)(|F_{1}|+|F_{2}|)\leq|F_{1}|=|C_{1}|-|\mathcal{I}|\leq|C_{1}|$. If instead $|F_{1}|<2|F_{2}|$, substituting yields $|\mathcal{I}|+(2/3)(|F_{1}|+|F_{2}|)>|\mathcal{I}|+|F_{1}|=|C_{1}|$.

4.2. Search Space Strategies

To reduce the memory required by dynamic programming, our implementation diverges from the pseudocode of Algorithm 1 in two ways. Specifically, we solve a restricted weight variant of $k$-FD and use a separate phase to recover the $s$-$t$–paths.

Weight restriction.

Rather than making one pass through the dynamic programming that infers the solution weights from the linear systems, we employ a multi-pass strategy. Each pass restricts the potential weight vector by fixing a subset of its entries to weights chosen from flow values in the input instance while leaving the remaining values as variables to be determined by the DP routine. Since we observe that most instances in the data set (see Section 5.1) admit decompositions whose paths saturate at least one arc, we begin by enumerating potential weight-vectors whose values are all taken from flow values. If the DP routine finds no solution for any of these vectors, we enumerate vectors in which exactly one value is undetermined and all other entries are again taken from flow values, \etcIf all these attempts should fail, the final pass in our strategy will leave all weights to be determined by the DP and hence is guaranteed to find a solution, should one exist. For most instances, however, a solution vector is quickly guessed, vastly reducing the number of candidate linear systems in the DP and thus saving both time and memory.

Path recovery.

Computing the path weights requires storing only the current and previous dynamic programming tables in memory. In contrast, recovering the paths requires storing all dynamic programming tables. For this reason, we first determine the weights $\mathbf{w}$ and then recover the paths by running the dynamic programming again with weights restricted to $\mathbf{w}$. This is equivalent to solving the $k$-Flow Routing problem in Section 6.

4.3. Pruning

Within the dynamic programming we employ a number of heuristics that help recognize algorithmic states that cannot lead to a solution.

Weight bounds.

We augment the weight constraints imposed by the linear systems with a set of routing-independent constraints. These take the form of upper and lower bounds ($B_{i}$ and $b_{i}$, respectively) on the $i$th smallest weight $w_{i}$.

First, we compute $b_{k}$ by noting that for any top-cut $C$ and any arc $a\in C$, only $k-|C|+1$ paths can be routed over $a$, i.e. $b_{k}\geq f(a)/(k-|C|+1)$. Then, we compute $B_{k}$ by finding the largest weight of any $s$-$t$–path. This can be done via a simple dynamic programming algorithm: for any node $v$, an $s$-$v$–path with weight $w$ requires a path of equal weight to an in-neighbor $u$ such that $w\leq f(uv)$.

To compute the bounds on the other weights, $i<k$, we let $B_{1}$ be the smallest flow value and $B_{i}$ be the smallest222If no such flow value exists, set $B_{i}=B_{k}$. flow value greater than $\sum_{j=1}^{i-1}B_{j}$. In other words, if $f(a)>\sum_{j=1}^{i-1}w_{j}$, then $a$ must be part of a path of weight at least $w_{i}$. Finally, we use these upper bounds to derive the $b_{i}$s. If the weights greater than $w_{i}$ sum to $W$, then by the pigeonhole principle $w_{i}\geq(F-W)/i$, where $F$ is the total flow. Thus, $w_{i}\geq(F-\sum_{j=i+1}^{k}B_{j})/(i+1)$.

For simplicity of implementation, we force the values of our weight vector to appear in non-decreasing order. In each dynamic programming table, for each linear system $L$ we run a linear program to see if there are any (rational) weights within the bounds that satisfy $L$. If not, we remove the entry of the dynamic programming table containing $L$. We also use these bounds before dynamic programming to catch partially determined weight vectors that will not be able to yield a solution. For example, the predetermined vector333A $*$ indicates an undetermined value. $\mathbf{w}=(1,3,*,*,5,11)$ is incompatible with a lower bound on the third entry greater than 5.

Storing linear systems.

Within the dynamic programming we store linear systems in row-reduced echelon form (RREF). When a new row $r$ is introduced to the system, we perform Gaussian elimination to convert the new row to RREF, checking for linear dependence and inconsistency in $O(k^{2})$ time. The iterative row reduction also has the advantage of revealing determined path weights even if the system is not fully determined, and thus we can check whether any such values are not positive integers. Furthermore, once the system is full rank, no additional computation needs to be done to recover the weights.

5. Experimental Results

In this section we empirically evaluate the efficiency and solution quality of Toboggan by comparing with the current state-of-the-art, Catfish, on a corpus of simulated RNA sequencing data. Additionally, we use Toboggan to theoretically validate the $k$-FD problem as a model for the transcript assembly task.

5.1. Experimental setup

Our experiments were run on a corpus of data used in previous experiments by Shao and Kingsford [15, 16]. The full data set contains over 17 million instances, each generated by simulating RNA-seq reads from a transcriptome and building the $s$-$t$–DAG and flow using the procedure described in Section 1. Two software packages were used to simulate reads: Flux Simulator [5], which used transcriptomes from humans, mice, and zebrafish, and Salmon [11], which only used human transcriptomes. We will refer to these four subsets as human, mouse, zebrafish, and human-salmon, respectively. The simulated reads allow us to know the true transcripts and thus the corresponding flow decomposition a priori; we will refer to this particular decomposition as the “ground truth”.

In our three experimental analyses (model validation, algorithm timing, and solution-quality) we report results on just the three smaller datasets (human, mouse, zebrafish). We terminated Toboggan on any run whose weight computation took longer than 800 seconds on zebrafish, and 50s on human and mouse 444Under these constraints, Toboggan timed out on $5\,136$ of the 4M instances in these smaller data sets.. After performing more exhaustive computations on these smaller data sets, we then used the larger human-salmon data set (13.3M) to add one more data point for the model validation experiment. Because of the size of human-salmon, we terminated Toboggan after only 1 second. This smaller time limit resulted in Toboggan timing out on a significantly larger portion of human-salmon than the other three ($3\%$ on human-salmon compared to $0.1\%$ on the other data sets), which would bias the results on human-salmon toward its simpler instances. Because of this, we omit the human-salmon results from the timing and solution-quality experiments.

For a small number of the 17M instances, the ground truth decomposition contains at least one path that appeared multiple times, corresponding to a subsequence repeated in different locations of the transcriptome. The data provides no means for mathematically or biologically distinguishing these; thus, we aggregated the duplicated paths into a single path, summing their weights. Additionally, we removed “trivial” instances in which the graph consisted of a single $s$-$t$–path; on such graphs Toboggan terminates during preprocessing without executing the $k$-FD algorithm described in this paper. We remark that this is a departure from the experimental setup of Catfish, which included such graphs, explaining some slight differences in our statistics. About $50\%$ of the 17M graphs in the corpus are trivial in this manner. The number of non-trivial graphs for each species is summarized in Table 1. All experiments were run on a dedicated system with an Intel i7-3770 (3.40GHz, 8 cores), an 8192 KB cache, and 32 GB of RAM.

5.2. Benchmarking

We start by analyzing the efficiency of Toboggan and Catfish, noting that this compares a Python implementation [7] with C++ code [15]. Their runtimes on the 1.4M non-trivial instances in the smaller data sets are shown side-by-side in Figure 2. We observe that the two implementations are both quite fast on the vast majority of instances: their median runtimes are $1.24$ milliseconds (ms) (Toboggan) and $3.47$ms (Catfish). However, the implementations have different runtime distributions—whereas Catfish is consistent, terminating between 2.3–4.6ms on $90\%$ of instances and never running longer than 1.3 seconds, Toboggan trades off faster termination, e.g. less than $2$ms on $80\%$ of instances, with a higher variance and a small chance of a much longer runtime, e.g. over $50$ seconds on $0.48\%$ of instances.

5.3. Model Validation

Previous papers that use flow decompositions to recover RNA sequences [16, 18] implicitly assume that the true RNA sequences correspond to a minimum size flow decomposition, as opposed to one with an arbitrary number of paths. Because Toboggan provably finds the minimum size of a flow decomposition, our implementation enables us to investigate exactly how often this assumption holds in practice.

Table 1 gives the percentage of all non-trivial instances whose ground truth decompositions are in fact minimum decompositions, as well as the percentage of ground truth decompositions we proved have non-optimal size. We conclude from this table that the Flow Decomposition problem is in fact a useful model for transcript assembly, which underscores the need for efficient algorithms to compute minimum decompositions.

5.4. Ground Truth Recovery

Though the ground truth flow decompositions are almost always of minimum size, it is biologically desirable to find a particular minimum size decomposition rather than an arbitrary one. In this section we investigate how often the decompositions output by Toboggan and Catfish are identical to the ground truth decomposition, restricting our attention to non-trivial instances in which the ground truth decomposition is of minimum size. Additionally, for those instances where each algorithm does not exactly recover the ground truth, we analyze the similarity of the (imperfect) path set to the ground truth, using the Jaccard index.

The table in Figure 3 summarizes the performance of the Toboggan and Catfish implementations in exactly computing the ground truth decompositions. Their behavior was quite similar, with slight differences at each decomposition size, and a $0.2\%$ advantage for Toboggan overall.

For instances where an algorithm’s output is not identical to the ground truth, an output can still recover some useful information if it is highly similar to the ground truth decomposition. With this in mind, we evaluate how similar each algorithm’s output is to the ground truth when they do not exactly match555There are 43,817 such instances for Catfish and 41,783 for Toboggan.. The plot in Figure 3 shows the distribution of the Jaccard index of each algorithm’s output compared to the ground truth paths.

6. Kernelization lower bounds

As described in Section 4.1, our implementation employs a graph reduction algorithm that rapidly identifies any trivial graph and immediately solves Flow Decomposition. Out of the 17M total instances, this preprocessing solves all of the roughly 8.7M trivial instances. On non-trivial instances, Toboggan then tries to identify the correct number $k$ of paths in an optimal solution. The naïve lower bound from the largest edge cut is equal to the correct value of $k$ in $96.5\%$ of the 8.6M non-trivial instances; incorporating the bound from Lemma 4.3 brings this up to $98.4\%$.

In the framework of parameterized complexity, it is therefore natural to ask whether $k$-FD admits a polynomial kernel. Below we answer the question in the negative, despite the real-world success of our preprocessing. Our proof strategy requires hardness reductions involving the following restricted variants of Flow Decomposition.

Problem 6.1**.**

$\mathcal{U}$-$k$-Flow Decomposition* ($\mathcal{U}$-$k$FD) \Input& An $s$-$t$–DAG, $G$ with an integral flow $f$, an integer $k$, and a set $\mathcal{U}\subset\mathbb{Z}$.

\Prob Does $(G,f)$ have a flow decomposition into $k$ paths whose weights are all members of $\mathcal{U}$?*

Problem 6.2**.**

$\mathcal{U}$-$k$-Flow Routing* ($\mathcal{U}$-$k$FR) \Input& An $s$-$t$–DAG, $G$, integral flow $f$, and $k$ integers $\mathbf{w}=(w_{1},\ldots,w_{k})$ taken from $\mathcal{U}\subset\mathbb{Z}$.

\Prob Is there a flow decomposition of $(G,f)$ into $k$ paths with respective weights $\mathbf{w}$?*

Lemma 6.3**.**

$\{1,2,4\}$-$k$-Flow Routing* is $\mathsf{NP}$-complete.*

Proof 6.4**.**

Consider an instance $(G,f,k)$ of $\{1,2,4\}$-$k$FD, which is $\mathsf{NP}$- complete [6]. Note that every $s$-$t$–path in a potential solution to such an instance can only take weights in $\{1,2,4\}$. This enables us to Turing-reduce $(G,f,k)$ to at most $k^{3}$ instances of $\{1,2,4\}$-$k$FR: we simply guess how many of the $k$ $s$-$t$–paths use each of the three possible values. It follows that $\{1,2,4\}$-$k$FR is $\mathsf{NP}$-complete.

In order to show that $k$-FD does not admit polynomial kernels, we will provide a cross-composition from $\{1,2,4\}$-$k$FR. We first need the following technical definition to set up the necessary machinery:

Definition 6.5** (Polynomial equivalence relation [2]).**

An equivalence relation $\mathcal{R}$ over $\Sigma^{*}$ is called a polynomial equivalence relation if the following hold:

1. (1)

there exists an algorithm that decides for $x,y\in\Sigma^{*}$ whether $x$ and $y$ are equivalent under $\mathcal{R}$ in time polynomial in $|x|+|y|$, and 2. (2)

for any finite set $S\subseteq\Sigma^{*}$ the index $|S/\mathcal{R}|$ is bounded by a polynomial in $\max_{x\in S}|x|$.

The benefit of a polynomial equivalence relation is that we can focus on collection of instances which share certain characteristics, as long as these characteristics do not distinguish too many instances. A simple example is that we can ask for input instances $(G_{i},f_{i},k_{i})$ in which all graphs $G_{i}$ have the same number of vertices and the values $k_{i}$ are the same.

Definition 6.6** (AND-cross-composition [2]).**

Let $L\subset\Sigma^{*}$ be a language, $\mathcal{R}$ a polynomial equivalence relation over $\Sigma^{*}$, and let $\Pi\subseteq\Sigma^{*}\times\mathbb{N}$ be a parameterized problem. An AND-cross composition of $L$ into $\Pi$ (under $\mathcal{R}$) is an algorithm that, given $\ell$ instances $x_{1},\ldots x_{t}\in\Sigma^{*}$ of $L$ belonging to the same equivalence class of $\mathcal{R}$, takes time polynomial in $\sum_{i=1}^{\ell}|x_{i}|$ and outputs an instance $(y,k)\in\Sigma^{*}\times\mathbb{N}$ such that

1. a)

the parameter $k$ is polynomially bounded in $\max_{1\leq i\leq\ell}|x_{i}|+\log\ell$, and 2. b)

we have that $(y,k)\in\Pi$ if and only if all instances $x_{i}\in L$.

We will now use the following theorem (abridged to our needs here) and the subsequent AND-cross-composition to prove that $k$-FD does not admit small kernels.

Proposition 6.7** (Bodlaender, Jansen, Kratsch [2]).**

If an $\mathsf{NP}$-hard language $L$ AND-cross-composes into a parameterized problem $\Pi$, then $\Pi$ does not admit a polynomial kernelization unless {\mathchoice{\hbox{\mathsf{NP}}}{\hbox{\mathsf{NP}}}{\mathsf{NP}}{\mathsf{NP}}}\subseteq{\mathchoice{\hbox{\mathsf{coNP}}}{\hbox{\mathsf{coNP}}}{\mathsf{coNP}}{\mathsf{coNP}}}/\emph{poly} and the polynomial hierarchy collapses.

Theorem 6.8**.**

$k$-Flow Decomposition* does not admit a polynomial kernel unless {\mathchoice{\hbox{\mathsf{NP}}}{\hbox{\mathsf{NP}}}{\mathsf{NP}}{\mathsf{NP}}}\subseteq{\mathchoice{\hbox{\mathsf{coNP}}}{\hbox{\mathsf{coNP}}}{\mathsf{coNP}}{\mathsf{coNP}}}/\emph{poly} and the polynomial hierarchy collapses.*

Proof 6.9**.**

Let $\mathcal{R}_{w}$ be the equivalence relation on instances of $\{1,2,4\}$-$k$FR where $(G_{1},f_{1},\mathbf{w}_{1})\equiv(G_{2},f_{2},\mathbf{w}_{2})$ if and only if $\mathbf{w}_{1}=\mathbf{w}_{2}$. Since each entry of $\mathbf{w}_{i}$ is in $\{1,2,4\}$, $\mathcal{R}_{w}$ has at most $k^{3}$ equivalence classes, and is a polynomial equivalence relation.

Let $x_{1},\ldots,x_{\ell}$ be instances of $\{1,2,4\}$-$k$FR all contained in the same equivalence class of $\mathcal{R}_{w}$, with $x_{i}=(G_{i},f_{i},\mathbf{w})$ and $\mathbf{w}=(w_{1},\ldots,w_{k})$ the common prescribed path weights. We denote the source and sink of $G_{i}$ by $s_{i}$ and $t_{i}$, respectively.

We construct an additional instance $x_{\ell+1}$ as follows: $G_{\ell+1}$ consists of two vertices $s_{\ell+1},t_{\ell+1}$, and $k$ arcs $a_{1},\ldots,a_{k}$ from $s_{\ell+1}$ to $t_{\ell+1}$. The flow $f_{\ell+1}$ has value $w_{i}$ on arc $a_{i}$. By construction $x_{\ell+1}=(G_{\ell+1},f_{\ell+1},\mathbf{w})$ is a positive instance of $\{1,2,4\}$-$k$FR; moreover, it has a unique decomposition into $k$ $s$-$t$–paths (up to isomorphism).

Before we describe the composition, we treat a technicality. If the total flow $F_{i}$ for any $f_{i}$ is different from $\sum_{j=1}^{k}w_{j}$, then $x_{i}$ is a negative instance. In this case, instead of a composition we output a trivial negative instance $y$ for $k$-FD. Otherwise, we compose the instances $x_{1},\ldots,x_{\ell+1}$ into a single composite instance $y=(G,f,k)$ of $k$-FD. To form $G$, we chain the $G_{i}$s together by identifying the vertex $t_{i}$ with the vertex $s_{i+1}$ for $1\leq i\leq\ell$. The resulting $G$ is an $s$-$t$–DAG with source $s=s_{1}$ and sink $t=t_{\ell+1}$. We define $f$ to label each arc in $G$ with the flow value from its original instance. Since each $x_{i}$ has the same total flow, $f$ is a flow on $G$.

We point out that property a) of Definition 6.6 is trivially satisfied. To see that property b) holds, first assume that all instances $x_{i}$ are positive. Since all these solutions consist of $k$ $s$-$t$–paths with the same prescribed values $w_{1},\ldots,w_{k}$, we can chain the individual flow decompositions together into $k$ $s$-$t$–paths in $G$. Accordingly, $y$ is then a positive instance. In the other direction, assume that $y$ has a solution, \ie $f$ can be split into exactly $k$ $s$-$t$–paths. Due to our inclusion of the instance $x_{\ell+1}$ in the construction of $y$, the respective values of the $s$-$t$–paths must be exactly $w_{1},\ldots,w_{k}$. But then restricting this global solution to each individual instance $x_{i}$ (since all $s$-$t$–paths meet at the identified source/sink cut vertices) produces a solution. We conclude that therefore all $x_{i}$ must have been positive instances, as required

Finally, our construction clearly takes time polynomial in $\sum_{i=1}^{\ell}|x_{i}|$ making it an AND-cross-composition of $\{1,2,4\}$-$k$FR into $k$-FD. Invoking Theorem 6.7, this proves that $k$-FD does not admit a polynomial kernel unless {\mathchoice{\hbox{\mathsf{NP}}}{\hbox{\mathsf{NP}}}{\mathsf{NP}}{\mathsf{NP}}}\subseteq{\mathchoice{\hbox{\mathsf{coNP}}}{\hbox{\mathsf{coNP}}}{\mathsf{coNP}}{\mathsf{coNP}}}/\emph{poly}.

We note that our construction in the proof of Theorem 6.8 produces an instance of $\mathcal{U}$-$k$FD, which can easily be Turing reduced into an instance of $\mathcal{U}$-$k$FR.

Corollary 6.10**.**

Unless {\mathchoice{\hbox{\mathsf{NP}}}{\hbox{\mathsf{NP}}}{\mathsf{NP}}{\mathsf{NP}}}\subseteq{\mathchoice{\hbox{\mathsf{coNP}}}{\hbox{\mathsf{coNP}}}{\mathsf{coNP}}{\mathsf{coNP}}}/\emph{poly}, the problems $\mathcal{U}$-$k$-Flow Decomposition and $\mathcal{U}$-$k$-Flow Routing do not admit polynomial kernels.

7. Hardness of assigning weights

Every solution computed by Toboggan in our experiments corresponded to a fully determined linear system in the final dynamic programming table. This means that we never had to run the expensive ILP solver to determine the weights; instead, they were computed in polynomial time with respect to $k$ using row reduction. This raises the question: is the linear system obtained from a decomposition into paths always fully determined? In the following we show that the answer must be ‘no’ in case of non-optimal decompositions. Not only can there be multiple weight-assignments for the same set of paths, recovering these weights is actually $\mathsf{NP}$-hard. Formally, we consider the problem $k$-Flow Weight Assignment.

Problem 7.1**.**

$k$-Flow Weight Assignment* ($k$FWA) \Input& An $s$-$t$–DAG $G$ with an integral flow $f$, and a prescribed set of $s$-$t$–paths, $P=p_{1},\cdots,p_{k}$.

\Prob Identify integral weights $\mathbf{w}=(w_{1},\dots w_{k})$ such that $(P,\mathbf{w})$ is a flow decomposition of $G$.*

Our proof that $k$FWA is $\mathsf{NP}$-complete relies on a reduction from Exact 3-Hitting Set, which is equivalent to monotone 1-in-3-SAT and known to be $\mathsf{NP}$-hard [14].

Problem 7.2**.**

Exact 3-Hitting Set* (X3HS) \Input& A finite universe $U=\{u_{1},\cdots,u_{n}\}$, and a collection $\mathcal{S}\subseteq\binom{U}{3}$ of subsets of $U$ of size 3.

\Prob Find a subset $X\subseteq U$ such that each element of $\mathcal{S}$ intersects (“hits”) $X$ exactly once, i.e. $\forall S\in\mathcal{S},|S\cap X|=1$.*

Theorem 7.3**.**

$k$-Flow Weight Assignment* is $\mathsf{NP}$-complete.*

Proof 7.4**.**

Given an instance $(U,\mathcal{S})$ of X3HS, we will construct an equivalent instance of $k$FWA. Our approach is to encode each element of $U$ using two “partner” $s$-$t$–paths in an $s$-$t$–DAG, $G$, whose weights sum to 3, and each triad in $\mathcal{S}$ as a different top-cut in $G$. We will route the partner paths and assign flow values so that the set of $s$-$t$–paths with weight 2 exactly correspond to elements in a solution to the hitting set problem.

We first construct $G$. Let $V(G)=\{s,v_{1},\ldots,v_{|\mathcal{S}|},t=v_{|\mathcal{S}|+1}\}$, where $v_{i}$ is associated with the $i$th triad in $\mathcal{S}$ for $1\leq i\leq|\mathcal{S}|$, and $s$/$t$ are source/sink vertices. Our construction will be such that each vertex other than $t$ has exactly one out-neighbor: $s$ has out-neighbor $v_{1}$ and $v_{i}$ has out-neighbor $v_{i+1}$. Create one $sv_{1}$ arc with flow value 3 for each element of $U$, and $(|U|-1)$ $v_{i}v_{i+1}$ arcs—one each of flow values 4 and 5, and $|U|-3$ with flow value 3. We now define a set of prescribed paths $P=\{p_{Y},\bar{p}_{Y}\}_{Y\in U}$. For each element $Y\in U$, the corresponding partner $s$-$t$–paths $p_{Y}$ and $\bar{p}_{Y}$ are routed over the same arc out of $s$. This guarantees that each pair of partner paths must receive weights summing to 3. At $v_{i}$, we route these paths to encode the corresponding triad $S_{i}=\{u_{1},u_{2},u_{3}\}$ as follows. The paths $p_{u_{1}}$, $p_{u_{2}}$, and $p_{u_{3}}$ are routed over the arc with flow 4 and $\bar{p}_{u_{1}}$, $\bar{p}_{u_{2}}$, and $\bar{p}_{u_{3}}$ go over the arc of flow 5. Now assign each element $u^{\prime}\in U\setminus S_{i}$ to a unique $v_{i}v_{i+1}$ arc $a_{u^{\prime}}$ of flow value 3 and route $p_{u^{\prime}}$ and $\bar{p}_{u^{\prime}}$ together over $a_{u^{\prime}}$. This construction is illustrated in Figure 4.

To complete the proof, we describe how a solution to this instance of $k$FWA yields a solution to the original instance of X3HS. Consider the possible values of the weights of the paths in $P$. By design, the out-arcs of $s$ force each pair of partner paths to have one of weight 1 and one of weight 2. Each triad $S_{i}$ is represented by two arcs out of $v_{i}$: one with flow value 4 ($a_{i}$), and the other with 5 ($\bar{a_{i}}$). Because exactly three prescribed $s$-$t$–paths are routed over $a$, and all our paths must have weight 1 or 2, this guarantees that exactly one $s$-$t$–path routed over $a$ has weight 2 (and the other two must have weight 1). Thus, finding a set of weights that solve this instance of $k$FWA is equivalent to choosing a set of paths to have weight 2 such that exactly one selected path is routed over each arc with flow value 4. But this is equivalent to choosing a set of elements (paths) such that each triad (arc) is hit by (incident to) exactly one of the chosen elements ($s$-$t$–paths of weight 2).

In the above reduction, we looked for a decomposition with $2\cdot|U|$ prescribed paths. However, if we do not prescribe paths, there is a flow decomposition with $|U|+1$ paths ($|U|-1$ of weight three, one of weight two, and one of weight one) and a corresponding linear system of full rank666This is true regardless of whether $(U,\mathcal{S})$ is a a yes-instance.. As such, it is possible that $k$FWA is not difficult when the prescribed paths on a yes-instance are part of an optimal decomposition.

Conjecture 7.5**.**

If $k$ is the minimum value for which $(G,f)$ has a flow decomposition of size $k$, then every integer-weighted solution has a corresponding linear system $L$ of rank $k$.

A direct consequence of this conjecture would be that the running time in Theorem 3.8 immediately improves to $\frac{1}{k!}\cdot 4^{k^{2}}0.649^{k}k^{O(1)}\cdot n$ since the ILP-solving step would never occur. In this context, we note that Vatinlen et al. [19] proved that when $k$ is minimum, every solution with real-valued weights has a corresponding linear system of full rank. However, their proof does not hold when the path weights are restricted to the integers.

8. Conclusion

We presented a holistic treatment of Flow Decomposition from the perspectives of parameterized complexity and algorithm engineering, resulting in a competitive solver, Toboggan. Our approach verifies that parameterized algorithms can (with sufficient engineering) be applied to real-world problems even in high-throughput situations. Our work also naturally leads to several practical and theoretical questions for further investigation.

On the practical side, we would like to understand the cases in which a minimal flow decomposition does not match the assembly problem’s ground truth, and how we might improve the recovery rate. The similarity in performance of Toboggan and Catfish in our experiments suggests that we need to refine either the problem formulation or our notion of minimality; in either case, more domain-specific knowledge is needed.

On the theory side, we ask whether there exists an fpt algorithm for Flow Decomposition with running time $k^{O(k)}n$ or better. In particular, it will be interesting to see whether the established techniques used to improve dynamic programming algorithms [1, 3] are applicable to our algorithm. Furthermore, if Conjecture 7.5 holds, it immediately implies a tighter upper bound on the running time of our algorithm, and might lead to further optimizations.

Acknowledgements. This work supported in part by the Gordon & Betty Moore Foundation’s Data-Driven Discovery Initiative through Grant GBMF4560 to Blair D. Sullivan.

Appendix A Pseudocode for $k$-Flow Decomposition

Appendix B Experimental results organized by species

Following the experimental setup of [16], in this section we report our results on each of the species data sets individually. Figures 5 and 6 give this breakdown for the aggregated results shown in Figure 3.