Reconstruction of viral population structure from next-generation sequencing data using multicommodity flows

AmpMCF algorithm

We consider an amplicon A as a multiset of reads such that each read r∈A has the same predefined starting and ending position in the genome start(A) and end(A), respectively. Two amplicons A ₁ , A ₂ are considered overlapped if and only if start(A₁) ≤ start(A ₂ ) < end(A ₁ ) ≤ end(A₂). A set of amplicons A = {A₁, ..., A _m } is said to be overlapping if and only if A _i and A_{i+ 1}overlap for i = 1...m-1. Given an overlapping set A, we define a partial order < on the set of reads R = A₁⋃...⋃A _m as follows: r <r' if and only if r∈A _i , r'∈A_i+1and r and r' are consistent over their overlap of length l_i,i+1 = end(A_i)- start(A_i+1)+1, i.e., the suffix of length l_i,i+1 of r coincides with the prefix of length l_i,i+1 of r'.

Given an overlapping set A = {A₁, ..., A _m }, we construct an (m+2)-staged directed vertex-weighted read-graph as follows: G = (V(G) = V₁ ⋃ ... ⋃ V _m ⋃ {s, t}, E(G), c), where each v∈ V _i , 1 ≤ i ≤ m corresponds to a distinct read r _v ∈A _i . An edge (u, v) ∈ E(G) if and only if either r _u <r _v or u = s, v∈ A₁ or u∈A _m , v = t. We also define the function c: V₁ ⋃ ... ⋃ V _m → [0,1], where c(v) denotes the frequency of the read represented by v ∈ V _i in amplicon A _i . It is evident that every full-size quasispecies that has a sequenced read from each amplicon A _i corresponds to an (s, t)-path in the graph G.

A bipartite clique of G is defined as a set of vertices C⊆V(G) such that C⊆V _i ⋃V_{i+ 1}for some i and every vertex from the set C⋂V _i is adjacent to every vertex from the set C⋂V_i+1.

Lemma 1. Consistent overlaps in amplicons A _i , A _i+1 correspond to disjoint bipartite cliques in G.

Proof. Suppose the contrary; then there exist vertices v, v' ∈ A _i and u, u' ∈A_i+1, such that r _v <r _u , r _v <r_u', r _v' <r _u , but it is not true that r _v' <r_u'. Since r_v'and r _u are comparable but r_v'and r_u'are not, the prefixes of length l_{i,i+ 1}of r _u and r_u'must not be consistent. This implies a contradiction with r _v <r _u and r _v <r_u'.

Using this simple finding, we transform the read graph G into a new "forked" edge-weighted directed read-graph H = (V(H), E(H), d) as follows. Consider each p × q-bipartite clique C = K _p,q of G not containing vertices s,t. C ⊆ A _i ⋃A_i+1for some i∈{1, ..., m-1}. Add a new "fork" vertex v _fork , delete all edges of the bipartite clique C and add edges from the sets {(u,v _fork ): u∈C⋂A _i } and {(v _fork , v): v∈C⋂A_i+1}. Define a new edge weight function d : E(H) → N as follows: d(uv _fork ) = c(u), d(v _fork v) = c(v), d(su) = d(vt) = 0. Figure 1 illustrates this transformation.

As for G, every full-size quasispecies corresponds to (s,t)-path in the forked read graph H. However, H is (2m+1)-staged directed graph with much fewer edges than G: for every bipartite clique K _p,q pq edges in G are replaced by only p+q edges in H. Since in network flow problems variables usually are associated with edges, this reduction is highly useful for the construction of the fast network flows-based algorithm for the quasispecies spectrum reconstruction problem.

The quasispecies reconstruction problem may be restated as the following covering problem:

Problem 1. Given a forked read graph H, cover H with a set of unique (s,t)-paths P _i with frequencies g _i such that the total frequency of paths is minimal and for every directed edge (u,v)∈ E(H) the sum of frequencies of paths containing (u, v) is at least d(uv).

This condition is guaranteed by the constraints (5). Constraints (2) and (3) are covering and flow conservation constraints, respectively. Constraints (4) guarantee that flows g ⁱ are unsplittable for every i = 1, ..., k, i.e. the edges carrying each flow form a simple (s,t)- path P _i in the forked read-graph H. In particular, the constraints imply that for every i = 1, ..., k the values $g_{uv}^i$ are equal for all edges of P _i . Therefore $g_{uv}^i$ can be interpreted as values proportional to frequencies of quasispecies i.

ShotMCF algorithm

where ε is the sequencing error rate, i.e. the probability of error per nucleotide. Note that the analogous formula was used in the quasispecies theory for the calculation of the probability of mutation between two different quasispecies [15].

Here $g_{uv}^i$ are flow variables. $g_v^i={\sum }_{uv\in E}{g}_{uv}^i$ are auxiliary variables used for the simplicity of notations, which represent total flow through vertices v∈V. The resulted flow is fractional and can split, thus allowing for read errors and mutations. (11)-(14) is a variant of MCF, where vertex capacity constraints are used instead of edge capacity constraints. Once the problem is solved, the frequency of each candidate quasispecies could be estimated using (9).

ShotMCF algorithm

Figures 2, 3 illustrates the comparison of ShotMCF and EM algorithms. The difference in performance between two algorithms is statistically significant for all distributions and sizes of data. The p-values of a Kruskal-Wallis test are summarized in Table 1.

ShotMCF statistically significantly outperforms EM on large data sets with geometric and skewed distributions, while the quality of EM is higher on small data sets. The quality of quasispecies reconstruction by EM, as implemented in ViSpA [8], declined with the increase in the dataset size for large numbers of reads, and was not significantly affected for ShotMCF. EM produced more accurate results on data sets with up to 300K reads generated using the uniform distribution. However, the trend of decrease in quality of EM estimations suggests that ShotMCF is more accurate on the larger data sets generated using the uniform distribution.

The accuracy of frequency estimation for variants with different abundances was analysed (Figure 4). All analysed sequences were partitioned into 5 groups according to their frequencies f: f ≤ 0.025, 0.025 <f ≤ 0.05, 0.05 <f ≤ 0.1, 0.1 <f ≤ 0.2 and f > 0.2. x-axis represents the groups and y-axis represents the average relative error of ShotMCF for each group. Frequencies of high-abundance variants were estimated more accurately. The accuracy of frequencies estimation increases monotonically with the abundance and stabilizes approximately at the abundance 0.1. The quality of frequency estimation increases, in general, with the number of reads in data set for all groups.

AmpMCF algorithm

The reconstructions obtained using AmpMCF (k = 12) and the Maximum Bandwidth (MB) algorithm proposed in [9] were compared. Maximum bandwidth is based on the packing formulation of the quasispecies spectrum reconstruction problem, and was shown to outperform the algorithm for quasispecies spectrum reconstruction from amplicon reads proposed in [10]. The following measures of quality of a solution were used:

1) RMSE

where P and S are probability distributions and $M=\frac{1}{2}\left(P+W\right)$.

Here, if CandQ is the set of quasispecies found by the algorithm and SimQ is the set of simulated quasispecies, then TruePositives = CandQ⋂SimQ, FalseNegatives = SimQ\CandQ and FalsePositives = CandQ\SimQ.

RMSE and JSD measure the quality of quasispecies frequencies estimation, and Sensitivity and PPV measure the quality of assembled quasispecies. Sensitivity is a measure of the positive identifications, which is defined as the percentage of correctly assembled quasispecies out of the population. PPV is a measure of the negative identification, which is defined as the percentage of correctly identified quasispecies over all assembled quasispecies.

Figures 5, 6, 7, 8 illustrate the comparison of AmpMCF and Maximum Bandwidth algorithms, and Table 2 summarizes statistical significance of the comparison of the algorithms.

According to RMSE, AmpMCF statistically significantly outperforms Maximum Bandwidth for all sizes of data sets with the geometric distributions, and for large data sets with the uniform distribution. Although AmpMCF exceeded in accuracy Maximum Bandwidth on the 5K and 20K datasets with the uniform distribution, the difference in performance was statistically insignificant, with p-value being slightly greater than the statistical significance threshold of 5%. For the skewed distribution the results were comparable without statistically significant advantage of one algorithm over the other.

According to JSD and PPV, ShotMCF statistically significantly outperforms Maximum Bandwidth on the 100K data sets with the geometric distribution, while Maximum Bandwidth had the lower JSD values on the 20K and 100K data sets with the skewed distribution. For all other measures, sizes and distributions the results were comparable with no statistically significant advantage of one algorithm over the other. The p-value for S could not be calculated for the 5K data sets with the skewed distribution, since both algorithms were equally sensitive on all test examples.

So AmpMCF outperformed Maximum Bandwidth in quasispecies frequencies estimation for populations with geometric and uniform distributions, while both algorithms showed a similar performance in quasispecies sequence reconstruction.

The low sensitivity of AmpMCF and Maximum Bandwidth on the 5K data set with the skewed distribution is associated with the erroneous reconstruction of low-abundance variants by both algorithms, with only a dominant variant being correctly identified. For larger data sets, populations with the skewed distributions were reconstructed much more successfully and variants with frequencies as low as 0.8% were detected. It should be also noted that low-frequency variants were detected with higher probability in populations with the geometric distribution (Figure 9). It suggests that the recoverability of low-frequency variants depends on the structure of a population and that the coverage provided by data sets of 5K reads is insufficient for low-frequency variants detection, if the population contains a dominant high-frequency variant.

In general, abundances of variants greatly affect their recoverability, with high-frequency variants being easier to detect (Figure 10). As above, all analysed sequences in Figure 10 were partitioned into 5 groups according to their frequencies f: f ≤ 0.025, 0.025 < f ≤ 0.05, 0.05 < f ≤ 0.1, 0.1 < f ≤ 0.2 and f > 0.2. x-axis represents the groups and y-axis the probability of variant recovery in each group. The probabilities of detection of variants within each group increase with the number of reads in a data set. While the probability of reconstruction of a variant with frequency less than 2.5% from the 5K data set was only 0.0092, all variants with frequencies greater than 20% were reconstructed from 20K and 100K data sets.

The accuracy of frequency estimation for detected variants with different abundances is illustrated by Figure 11. As for ShotMCF, the accuracy of frequency estimation increases with the abundance and stabilizes approximately at the abundance 0.1. In general, the accuracy of frequency estimation increases with the number of reads in a data set for all groups except for the group of low-frequency variants. The small value of RE for low-frequency variants from the 5K data sets can be explained with a low detection rate of such variants, which renders their RE undefined.