Estimating Sequence Similarity from Read Sets for Clustering Next-Generation Sequencing data

TL;DR

This paper introduces a novel method to estimate sequence similarity directly from read sets in next-generation sequencing data, bypassing the need for sequence assembly and improving clustering accuracy especially at low coverage.

Contribution

It adapts the Monge-Elkan similarity for sequence read sets, providing a more effective and computationally feasible similarity measure for clustering without assembly.

Findings

Better approximation of true sequence similarities at low coverage

Improved clustering results compared to assembly-based methods

Avoids NP-hard sequence assembly problem

Abstract

To cluster sequences given only their read-set representations, one may try to reconstruct each one from the corresponding read set, and then employ conventional (dis)similarity measures such as the edit distance on the assembled sequences. This approach is however problematic and we propose instead to estimate the similarities directly from the read sets. Our approach is based on an adaptation of the Monge-Elkan similarity known from the field of databases. It avoids the NP-hard problem of sequence assembly. For low coverage data it results in a better approximation of the true sequence similarities and consequently in better clustering, in comparison to the first-assemble-then-cluster approach.

Tables

Table 1. Table 1: Overview of the datasets used in the experiments.

Name	Source	Read generation	Strand known	5’ to 3’ known	$n$	$\alpha$	$l$	Time-limit
Influenza	ENA	i.i.d., uniform distr.	✗	✗	$13$	$0.1$ to $100$	$3$ to $500$	$2$ hours
Various	ENA	i.i.d., uniform distr.	✗	✗	$18$	$0.1$ to $100$	$3$ to $500$	$2$ hours
Hepatitis	ENA	(art, )	✗	✓	$81$	$\{10,30,50\}$	$\{30,70,100\}$	$1$ day
Chromosomes	(1kgenomes, )	real-world	✗	✓	$23$	$4.32$	$76$	$1$ day

Table 2. Table 2: Runtime, Pearson’s correlation coefficient between distance matrices and Fowlkes-Mallows index for k = 4 𝑘 4 k=4 and k = 8 𝑘 8 k=8 . The ‘reference’ method calculates distances from the original sequences. We show only two assembly algorithms that gave the highest correlation and the better algorithm of pairs MES / MESS and MESSG / MESSGM .

Dataset	method	finished	$\frac{\text{assem.}}{\mathrm{ms}}$	$\frac{\text{distances}}{\mathrm{ms}}$	$\frac{\text{UPGMA}}{\mathrm{ms}}$	$\frac{\text{NJ}}{\mathrm{ms}}$	corr.	$\frac{\text{UPGMA}}{B_4}$	$\frac{\text{UPGMA}}{B_8}$	$\frac{\text{NJ}}{B_4}$	$\frac{\text{NJ}}{B_8}$
Influenza	reference	112/112	$0$	$4,893.33$	$4.13$	$11.29$	$1$	$1$	$1$	$1$	$1$
	$\max (\left|R_A\right|,\left|R_B\right|)$	112/112	$0$	$186.87$	$0.58$	$0.96$	$.801$	$.67$	$.32$	$.66$	$.32$
	${Dist}_{\text{MESS}}$	106/112	$0$	$854,821.17$	$0.63$	$0.41$	$.942$	$1$	$.85$	$1$	$.83$
	${Dist}_{\text{MESSG}}$	96/112	$0$	$859,725.65$	$0.23$	$0.33$	$.98$	$.99$	$1$	$1$	$.99$
	${Dist}_{\text{MESSGq}}$	112/112	$0$	$87,278.04$	$0.22$	$0.42$	$.976$	$1$	$.99$	$1$	$.99$
	SSAKE	51/112	$2,983.75$	$12,352.33$	$0.25$	$0.25$	$.72$	$.46$	$.59$	$.55$	$.66$
	SPAdes	31/112	$12,438.81$	$13,684$	$0.16$	$0.23$	$.722$	$.56$	$.68$	$.66$	$.72$
Various	reference	112/112	$0$	$71,019.2$	$4.81$	$11.1$	$1$	$1$	$1$	$1$	$1$
	$\max (\left|R_A\right|,\left|R_B\right|)$	112/112	$0$	$441.49$	$0.62$	$1.16$	$.901$	$.67$	$.66$	$.85$	$.92$
	${Dist}_{\text{MESS}}$	72/112	$0$	$1,212,284.96$	$0.33$	$0.53$	$.928$	$.63$	$.79$	$.87$	$.98$
	${Dist}_{\text{MESSGM}}$	64/112	$0$	$1,482,902.92$	$0.3$	$0.56$	$.934$	$.65$	$.89$	$.87$	$.94$
	${Dist}_{\text{MESSGq}}$	110/112	$0$	$660,983.1$	$0.4$	$0.54$	$.929$	$.66$	$.76$	$.84$	$.98$
	Edena	57/112	$1,112.84$	$156,249.28$	$0.39$	$0.61$	$.447$	$.47$	$.21$	$.48$	$.27$
	Velvet	109/112	$1,726.28$	$165,184.43$	$0.36$	$0.61$	$.39$	$.46$	$.19$	$.48$	$.31$
Hepatitis	reference	9/9	$0$	$2,300,822.78$	$33.22$	$51.89$	$1$	$1$	$1$	$1$	$1$
	$\max (\left|R_A\right|,\left|R_B\right|)$	9/9	$0$	$8,301.22$	$7.78$	$18.44$	$.181$	$.55$	$.37$	$.72$	$.83$
	${Dist}_{\text{MES}}$	9/9	$0$	$12,699,975.67$	$37.89$	$82.44$	$.833$	$1$	$.95$	$1$	$.92$
	${Dist}_{\text{MESSG}}$	9/9	$0$	$26,958,877.56$	$376.22$	$606.67$	$.964$	$1$	$.95$	$1$	$.92$
	${Dist}_{\text{MESSGq}}$	9/9	$0$	$616,507.11$	$3.78$	$7.89$	$.901$	$.94$	$.96$	$1$	$.94$
	Velvet	9/9	$22,599.89$	$3,016,592$	$1.67$	$4.11$	$.783$	$.8$	$.84$	$.96$	$.86$
	SPAdes	9/9	$97,626.44$	$1,882,311.56$	$304.89$	$119.11$	$.578$	$.82$	$.83$	$.95$	$.83$
Chromosomes	reference	1/1	$0$	$886,125$	$7$	$20$	$1$	$1$	$1$	$1$	$1$
	$\max (\left|R_A\right|,\left|R_B\right|)$	1/1	$0$	$868$	$1$	$2$	$.331$	$.64$	$.4$	$.61$	$.3$
	${Dist}_{\text{MES}}$	1/1	$0$	$14,006,230$	$0$	$1$	$.886$	$.42$	$.26$	$.6$	$.28$
	${Dist}_{\text{MESSG}}$	1/1	$0$	$29,676,206$	$0$	$1$	$.364$	$.54$	$.47$	$.64$	$.35$
	${Dist}_{\text{MESSGq}}$	1/1	$0$	$219,485$	$0$	$2$	$.366$	$.53$	$.47$	$.64$	$.35$
	SSAKE	1/1	$57,584$	$70,383$	$1$	$1$	$.718$	$.47$	$.19$	$.81$	$.15$
	SPAdes	1/1	$27,009$	$125,688$	$1$	$2$	$.444$	$.41$	$.18$	$.58$	$.23$

Equations

$$\mathrm{dist}(A,B)\le \max \{|A|,|B|\}.$$

$$|A|\gg l$$

$$|R_A|=\frac{\alpha }{l}|A|$$

$$b_j={\mathrm{argmin}}_{b_k\in R_B}\mathrm{dist}(a_i,b_k).$$

$${\mathrm{Dist}}_{\text{ME}}(R_A,R_B)=\frac{1}{|R_A|}\sum _{a_i\in R_A}\underset{b_j\in R_B}{\min }\mathrm{dist}(a_i,b_j).$$

$${\mathrm{Dist}}_{\text{MES}}(R_A,R_B)=\frac{1}{2}\left({\mathrm{Dist}}_{\text{ME}}(R_A,R_B)+{\mathrm{Dist}}_{\text{ME}}(R_B,R_A)\right).$$

$${\mathrm{Dist}}_{\text{MESS}}(R_A,R_B)=\max \{|R_A|,|R_B|\}{\mathrm{Dist}}_{\text{MES}}(R_A,R_B).$$

$$t<\frac{l}{2}$$

$$p={\left(1-\frac{1+2t}{|A|}\right)}^{\frac{|A|\cdot \alpha }{l}}\le \epsilon .$$

$$p=1-\frac{2t+1}{|A|}\frac{|A|\cdot \alpha }{l}+{\epsilon }^{\prime }=1-(2t+1)\frac{\alpha }{l}+{\epsilon }^{\prime }\le \epsilon .$$

$$t=\frac{1}{2}\left(\frac{l}{\alpha }-1\right).$$

$$g(x)=\{\begin{array}{ll}0,& \text{mbox}if0\le x\le t-1,\\ 2\frac{x-t+1}{l+1-2t},& \text{mbox}ift-1<x\le l-t,\\ 2,& \text{mbox}ifl-t<x<l.\end{array}$$

$$\mathrm{dist}(a_i,b_j)\ge \theta$$

$$\mathrm{dist}(a,b)\ge \frac{1}{6}{\mathrm{dist}}_q(a,b).$$

$$\mathrm{dist}(a,b)$$

$$\ge 1+\frac{1}{6}\left({\mathrm{dist}}_q(a,b)-6\right)=\frac{1}{6}{\mathrm{dist}}_q(a,b).$$

Full text

Our approach is based on an adaptation of the Monge-Elkan similarity known from the field of databases. It avoids the NP-hard problem of sequence assembly. For low coverage data it results in a better approximation of the true sequence similarities and consequently in better clustering, in comparison to the first-assemble-then-cluster approach.

1 Introduction

Sequencing means reading the sequence of elements that constitute a polymer, such as the DNA. The human genome project (lander.nature.11, ) completed in 2003 was a prime example of sequencing, resulting in the identification of almost the entire genomic sequence (over 3 billion symbols) of a single human. Sequencing becomes technologically difficult as the length of the read sequence grows. The common principle of new-generation sequencing (NGS) is that only very short substrings (10’s to 100’s of symbols) are read at random positions of the sequence of interest. It is usually required that the number of such read substrings, called reads, is such that with high probability each position in the sequence is contained in multiple reads; the number of such reads is termed coverage. The complete sequence is determined by combinatorial assembly of the substrings guided by their suffix-prefix overlaps. For example, one possible assembly of reads $\{\textsf{AGGC},\textsf{TGGA},\textsf{GCT}\}$ is AGGCTGGA. Short reads imply low cost of wet-lab sequencing traded off with high computational cost of assembly. Indeed the assembly task can be posed as searching the Hamiltonian path in a graph of mutual overlaps.

One of the central tasks in computational biology is to infer phylogenetic trees, which typically amounts to hierarchical clustering of genomes. When they are represented only through sequences read-sets, the bioinformatician is forced to reconstruct the sequences from the read-sets prior to clustering. This of course entails the solution of the NP-hard assembly problem for each data instance with little guarantees regarding the quality of the resulting putative sequence. This motivates the question whether the assembly step could be entirely avoided. We address this question here by proposing a similarity function computable directly on the read sets, that should approximate the true similarities on the original sequences.

In the next section we provide a brief overview of the related work and relevant methods. In Sect. 3, we design the similarity (or, reversely distance) function. Then we provide a brief theoretical analysis of it. In Sect. 5 we compare it to the conventional approach on genomic data and then we conclude the paper.

2 Related Work

Related work includes studies on clustering NGS data (e.g. seed ; pipe ; kchouk2016clustering ). They however deal with clustering reads and we are not aware of a previous attempt to cluster read-sets. The paper (distance, ) proposes a similarity measure for NGS data, but again it operates on the level of reads. The previous work (zelezny.14, ) and (jalovec.14, ) also aims at avoiding the assembly step in learning from NGS data but these studies concern supervised classification learning and they do not elaborate on read-set similarity.

This work extends (ida2016, ), which corresponds to Sects. 3.1, 3.3, 3.4 and 3.5. However this work fixes several issues that made the original approach hard to apply on real data. In this paper we also provide a more detailed experimental evaluation.

Approximate string matching area is related to the method we provide. Here the goal is to find strings similar to a pattern in a dictionary. Namely we employ the $q$-gram distance (ukkonenQgram, ). Our method also includes the Monge-Elkan distance (mongeelkan, ) known from the field of databases.

3 Proposed method

The functor $|.|$ will denote the absolute value, cardinality and length (respectively) for a number, set, and string argument. Let $\mathop{\mathrm{dist}}(A,B)$ denote the Levenshtein distance (levenstein, ) between strings $A$ and $B$. The function measures the minimum number of edits (insertions, deletions, and substitutions) needed to make the strings identical, and is a typical example of a sequence dissimilarity measure used in bioinformatics. It is a property of the distance that

[TABLE]

We will work with constants $l\in\mathbb{N},\alpha\in\mathbb{R}^{+}$ called the read length and coverage, respectively, which are specific to a particular sequencing experiment. A read-set $R_{A}$ of string $A$ such that

[TABLE]

is a multiset of111Should the right hand side be non-integer, we neglect its fractional part.

[TABLE]

substrings sampled i.i.d. with replacement from the uniform distribution on all the $|A|-l+1$ substrings of length $l$ of $A$. Informally, the coverage $\alpha$ indicates the average number of reads covering a given place in $A$.

Throughout the paper we hold assumption that the coverage $\alpha$ and read length $l$ are constant. The constant read length assumption can be justified by the way how some of the sequencing machines work. They read one nucleotide of each read every iteration.

Our goal is to propose a distance function $\mathop{\mathrm{Dist}}(R_{A},R_{B})$ that approximates $\mathop{\mathrm{dist}}(A,B)$ for read-sets $R_{A}$ and $R_{B}$ of arbitrary strings $A$ and $B$. We also want $\mathop{\mathrm{Dist}}(R_{A},R_{B})$ to be more accurate and less complex to calculate than a natural estimate $\mathop{\mathrm{dist}}(\hat{A},\hat{B})$ in which the arguments represent putative sequences reconstructed from $R_{A}$ and $R_{B}$ using assembly algorithms such as (abyss, ; edena, ; ssake, ; spades, ; velvet, ).

Alphabet that we use in our model is $\Sigma=\{\textsf{A},\textsf{T},\textsf{C},\textsf{G}\}$. If $a_{i}$ is a read, then $\overline{a}_{i}$ denotes its complement, where A (resp. C) and T (resp. G) are interchanged. For example complement of read ATTCG is read TAAGC. Similarly we define reversed read, denoted

$a_{i}$

. For example reverse of ATTCG is GCTTA.

3.1 Base Case: Which Reads Belong Together

A natural approach to instantiate $\mathop{\mathrm{Dist}}(R_{A},R_{B})$ is to exploit the $|R_{A}||R_{B}|$ pairwise Levenshtein distances between the reads in $R_{A}$ and $R_{B}$. Most of those values are useless because they match reads from completely different parts of sequences $A$ and $B$. Therefore we want to account only for those pairs which likely belong together.

If we seek a read from $R_{B}$ that matches a read $a_{i}\in R_{A}$, we make the assumption that the most similar read $b_{j}\in R_{B}$ is the one that we look for (see Fig. 1), i.e.,

[TABLE]

To calculate the distance from $R_{A}$ to $R_{B}$, we average over all reads from $R_{A}$:

[TABLE]

This idea was presented in (mongeelkan, ) for searching duplicates in database systems. The method is known as the Monge-Elkan similarity222Here we alter the Monge-Elkan similarity into a distance measure. The standard way of using Monge-Elkan is as a similarity measure with $\min$ replaced by $\max$ and distance calculation by similarity calculation. (hence the ME label) and entails a simple but effective approximation algorithm.

$\mathop{\mathrm{Dist}}_{\textsf{ME}}(R_{A},R_{B})$ is non-symmetric in general, which is undesirable given that the approximated distance $\mathop{\mathrm{dist}}(A,B)$ is known to be symmetric. Therefore we define a symmetric version by averaging both directions

[TABLE]

3.2 Strand and orientation

In practical setting we do not know which DNA strand do the reads come from. If we match read $a_{i}$ with read $b_{j}$, there are two possible matchings. If reads come from the complementary strands, we need to calculate \mathop{\mathrm{dist}}(a_{i},\reflectbox{\overline{b}_{j}}) besides $\mathop{\mathrm{dist}}(a_{i},b_{j})$. We consider only the option that leads to lower distance.

Some sequencing techniques lead to loss of information which end of read belongs to 3’-end and which end belongs to 5’-end. In that case we have two other independent options. Reads $a_{i}$ and $b_{j}$ either match or their reverses match. We calculate both distances and the orientation that leads to lower distance is used as the most likely one.

Based on the sequencing setting we have up to four options how to match reads $a_{i}$ and $b_{j}$. The options that need to be considered are subset of $\mathop{\mathrm{dist}}(a_{i},b_{j})$, $\mathop{\mathrm{dist}}(a_{i},\overline{b}_{j})$, \mathop{\mathrm{dist}}(a_{i},\reflectbox{b_{j}}) and \mathop{\mathrm{dist}}(a_{i},\reflectbox{\overline{b}_{j}}).

3.3 Distance Scale

Consider duplicating a non-empty string $A$ into $AA$ and assume $R_{AA}=R_{A}\cup R_{A}$. Typically for a $B$ similar to $A$ we expect that $\mathop{\mathrm{dist}}(AA,B)>\mathop{\mathrm{dist}}(A,B)$ but the (symmetric) Monge-Elkan distance will not change, i.e. $\mathop{\mathrm{Dist}}_{\textsf{MES}}(R_{AA},R_{B})=\mathop{\mathrm{Dist}}_{\textsf{MES}}(R_{A},R_{B})$, indicating a discrepancy that should be rectified.

In fact, $\mathop{\mathrm{Dist}}_{\textsf{MES}}$ has the constant upper bound $l$, which is because it is the average (c.f. (4) and (5)) of numbers no greater than $l$ (see (1)). On the other hand, $\mathop{\mathrm{dist}}(A,B)$ has a non-constant upper bound $\max\{|A|,|B|\}$ as by (1).

To bring $\mathop{\mathrm{Dist}}_{\textsf{MES}}(A,B)$ on the same scale as $\mathop{\mathrm{dist}}(A,B)$, we should therefore multiply it by the factor $\max\{|A|,|B|\}/l=\max\{|A|/l,|B|/l\}$. By (3) we have $|A|=\frac{l}{\alpha}|R_{A}|$, yielding the factor $\max\{|R_{A}|/\alpha,|R_{B}|/\alpha\}$, in which $\alpha$ is a constant divisor which can be neglected in a distance function. Therefore, we modify the read distance into

[TABLE]

3.4 Margin Gaps

Consider the situation in Fig. 2 showing two identical sequences each with one shown read. The Levenshtein distance between the two reads is non-zero due to the one-symbol trailing (leading, respectively) gap of the top (bottom) read caused only by the different random positions of the reads rather than due to a mismatch between the sequences. Thus there is an intuitive reason to pardon margin gaps up to certain size $t$

[TABLE]

when matching reads. Here, $t$ should not be too large as otherwise the distance could be nullified for pairs of long reads with small prefix-suffix overlaps, which would not make sense.

To estimate a good value for $t$, consider sequence $A$ and its sampled read-set $R_{A}$. We now sample an additional read $a$ of length $l$ from $A$. Ideally, there should be a zero-penalty match for $a$ in $R_{A}$ as $a$ was sampled from the same sequence as $R_{A}$ was. This happens iff there is a read in $R_{A}$ sampled from the same position in $A$ as $a$, or from a position shifted by up to $t$ symbols to the left or right as then the induced gaps are penalty-free. Since $R_{A}$ is an uniform-probability i.i.d. sample from $A$, the probability that a particular read from $R_{A}$ starts at one of these $1+2t$ positions is333Strictly speaking, this reasoning is incorrect if read $a$ is drawn from a place close to $A$’s margins, more precisely, if it starts in fewer than $t$ ($t+l$, respectively) symbols from $A$’s left (right) margin, as then not all of the $2t$ shifts are possible. This is however negligible due to Ineq. (2). $\frac{1+2t}{|A|}$. We want to put an upper bound $\varepsilon>0$ on the probability that this happens for none of the $|R_{A}|=\frac{\alpha}{l}|A|$ reads in $R_{A}$:

[TABLE]

Consider the first-order Taylor approximation $(1+x)^{n}=1+nx+\varepsilon^{\prime}$ where the difference term $\varepsilon^{\prime}>0$ decreases with decreasing $|x|$. Due to Ineq. (6) and (2), $\frac{1+2t}{|A|}$ is small and we can apply the approximation on the above formula for $p$, yielding

[TABLE]

For simplicity, we choose $\varepsilon=\varepsilon^{\prime}$. The smallest gap size $t$ for which the inequality is satisfied is obtained by solving $1-(2t+1)\frac{\alpha}{l}=0$, yielding

[TABLE]

This choice of $t$ matches intuition in that with larger read-length $l$ we can allow a larger grace gap $t$ but with larger coverage $\alpha$, $t$ needs not be so large as there is a higher chance of having a suitably positioned read in the read-set. Another way to look at it is to realize that reads in a read-set are approximately $\frac{l}{\alpha}$ positions from each other. Consider matching read $a$ to reads from $R_{A}$. If there is a read $a_{1}\in R_{A}$ requiring gap larger than $\frac{l}{2\alpha}$ to match $a$, then there will typically be another read $a_{2}\in R_{A}$ requiring gap at most $\frac{l}{2\alpha}$ (see Fig. 3). From (6) and (7) we see that this method is applicable only when $\alpha>(\frac{1}{l}+1)^{-1}$, which is bit less than $1$. However the results start to be nonzero for $\alpha>(\frac{1}{l}+2)^{-1}$, which is bit less than $0.5$.

We implemented the grace margin gaps into a further version $\mathop{\mathrm{Dist}}_{\textsf{MESSG}}$ of the constructed distance function, which required only a small change to the standard Wagner-Fischer algorithm (wagnerfischer, ).444The dynamic programming algorithm for calculating the Levenshtein distance (levenstein, ) is commonly called Wagner-Fischer algorithm (wagnerfischer, ). When we refer to sequence alignment problem in bioinformatics, this algorithm is often called Needleman-Wunsch algorithm (needlemanwunsch, ). When the algorithm is filling the first or the last row and column of the table, margin gaps up to $t$ symbols are not penalized. Larger margin gaps are penalized in a way that satisfies the constraint that the distance between a word $a$ and an empty word is $|a|$. In particular, the standard linear gap penalty is replaced with a piecewise linear function that gives cost of margin gap at $x$-th position

[TABLE]

3.5 Missing Read

Sometimes there is no good match for read $a_{i}$ in $R_{B}$. During evolution the substring that contained $a_{i}$ may have been inserted into $A$ or may have vanished from $B$. Therefore if

[TABLE]

for some reads $a_{i}$ and $b_{j}$ and threshold $\theta$, we consider $a_{i}$ and $b_{j}$ to be dissimilar and we force their distance to be $l$. (See Fig. 4.)

Threshold $\theta$ should be a linear factor of the maximal distance between two sequences of length $l$, i.e. $\theta=\theta^{\prime}\cdot l$. Value of $\theta^{\prime}$ should reflect the probability that the read is in one sequence and not in the other. Because the true probability is hidden, it needs to be determined empirically.

The distance function equipped with the missing read detection as described gives rise to the version denoted as $\mathop{\mathrm{Dist}}_{\textsf{MESSGM}}$.

3.6 Sampling

Unlike the assembly, the approach in this paper does not require high coverage to produce good quality data. Therefore for data sequenced with high coverage we can randomly sample only a small amount of reads to improve runtime. From the experiments we will see that coverage between $2$ and $3$ provides a good balance between runtime and quality of results.

3.7 Embedding

Another approach that can further reduce the running time is based on using embedding for finding a good candidate for the best match. The minimum in (4) does not need to be exact, but may be approximate. We therefore embedd reads to a space of lower dimension, where we can calculate the distance faster. Then we find reads minimizing the embedded distance. For this small set of reads we evaluate the true Levenshtein distance with grace margin gap penalty $\mathop{\mathrm{Dist}}_{\textsf{MESSG}}$ to find better estimate for use in (4).

The embedding that we use is based on $q$-gram profile and $q$-gram distance as formalized in (ukkonenQgram, ). $q$-gram is any string from $\Sigma^{q}$. Then $q$-gram profile of string $a$ is a vector $\vec{Q}(a)\in\mathbb{N}_{0}^{|\Sigma|^{q}}$ that contains number or occurences for all possible $q$-grams in $a$. Finally $q$-gram distance $\mathop{\mathrm{dist}}_{q}(a,b)$ of strings $a$ and $b$ is Manhattan distance of their $q$-gram profiles, i.e. $\mathop{\mathrm{dist}}_{q}(a,b)=\|\vec{Q}(a)-\vec{Q}(b)\|_{1}$.

We use $q$-gram distance to approximate the Levenshtein distance to find candidates that minimize the right hand side of (4). The $q$-gram distance forms a lower bound on edit distance (see Theorem 3.1). Similarly to BLAST (blast, ) algorithm we use $q=3$, which provides a good balance between the runtime and quality of estimates. The dimension of embedded space of reads is equal to $|\Sigma|^{q}$, which is $64$ in our case. In the next sections we will refer to this modification of distance as $\mathop{\mathrm{Dist}}_{\textsf{MESSGq}}$

Theorem 3.1

For any reads $a$ and $b$ and $q=3$ holds

[TABLE]

Proof

We will prove the theorem by mathematical induction on Levenshtein distance between $a$ and $b$. Suppose that $\mathop{\mathrm{dist}}(a,b)=k$.

• •

If $k=0$, then both $\mathop{\mathrm{dist}}(a,b)=\mathop{\mathrm{dist}}_{q}(a,b)=0$ and theorem holds.

• •

Suppose that theorem holds for any reads with distance at most $k-1$. There exists a sequence of $k$ operations insert, delete and replace that transforms $a$ into $b$. Denote $b^{\prime}$ the result that we obtain after applying the first $k-1$ operations. Then $\mathop{\mathrm{dist}}(a,b^{\prime})=\mathop{\mathrm{dist}}(a,b)-1$. Consider the last operation and calculate the maximal value of $\mathop{\mathrm{dist}}_{q}(a,b)-\mathop{\mathrm{dist}}_{q}(a,b^{\prime})$.

In case of insertion (or vice versa deletion) we obtain in worst case instead of $q$-grams $abc$,$bcd$ $q$-grams $abX$,$bXc$,$Xcd$ (or vice versa). The $q$-gram distance grows by at most $5$. In case of mismatch $q$-grams $abX$,$bXc$,$Xcd$ are replaced by $q$-grams $abY$,$bYc$,$Ycd$. The $q$-gram distance grows at most by $6$. Therefore $\mathop{\mathrm{dist}}_{q}(a,b)-\mathop{\mathrm{dist}}_{q}(a,b^{\prime})\leq 6$. Putting this together with the induction assumption leads to

[TABLE]

By mathematical induction the theorem holds. ∎

4 Theoretical Analysis

4.1 Asymptotic Complexity

Calculating $\mathop{\mathrm{dist}}(A,B)$ for sequences $A$ and $B$ requires $\Theta(|A||B|)$ operations if we use the standard Wagner-Fischer dynamic programming algorithm (wagnerfischer, ). This algorithm also requires $\Theta(\min(|A|,|B|))$ memory as we are interested only in distance and not in the alignment. To calculate $\mathop{\mathrm{Dist}}_{\textsf{ME}}$ we need to know the distances between all pairs of reads, so we have to evaluate (see (3)) $\frac{\alpha}{l}|A|\frac{\alpha}{l}|B|$ distances where each one requires $l^{2}$ operations. Therefore $\alpha^{2}|A||B|$ operations are required. For the symmetric version $\mathop{\mathrm{Dist}}_{\textsf{MES}}$ we make $2\alpha^{2}|A||B|$ operations, which can be reduced to $\alpha^{2}|A||B|$ operations and $\Theta(l+\frac{\alpha}{l}(|A|+|B|))$ memory. Further modifications (MESS, MESSG, MESSGM) do not change the asymptotic complexity.

Sampling, as described in Sect. 3.6 reduces runtime by $\alpha^{2}$ factor to $\Theta(|A||B|)$, assuming that the final coverage is a constant. Method MESSGq does not give any theoretical guarantee on how many pairwise edit distances are explored. However if we assume that this number is a small constant, we get runtime of $\Theta(|R_{A}||R_{B}|+l^{2}(|R_{A}|+|R_{B}))$.

The constants $\alpha$ and $l$ are determined by the sequencing technology and the independent complexity factors are $|A|$ and $|B|$. To calculate the distance in the conventional way as $\mathop{\mathrm{dist}}(\hat{A},\hat{B})$ requires to reconstruct $\hat{A}$ and $\hat{B}$ from the respective read-sets through an assembly algorithm. This is an NP-hard problem which becomes non-tractable for large $|A|$ and $|B|$, and which is avoided by our approach.

4.2 Metric Properties

$\mathop{\mathrm{Dist}}_{\textsf{MES}}$ as well as the later versions are all symmetric and non-negative but none of the proposed versions satisfies the identity condition ($\mathop{\mathrm{dist}}(a,b)=0$ iff $a=b$) or the triangle inequality, despite being based on the Levenshtein distance $\mathop{\mathrm{dist}}$, which is a metric. For example, let $R_{A}=\{\textsf{ATC},\textsf{ATC},\textsf{GGG}\}$, let $R_{B}=\{\textsf{ATA},\textsf{GGG}\}$, and let $R_{C}=\{\textsf{CTA},\textsf{GGG}\}$. Then $\mathop{\mathrm{Dist}}_{\textsf{MES}}(R_{A},R_{B})=\frac{7}{12}$, and $\mathop{\mathrm{Dist}}_{\textsf{MES}}(R_{B},R_{C})=\frac{1}{2}$ but $\mathop{\mathrm{Dist}}_{\textsf{MES}}(R_{A},R_{C})=\frac{14}{12}>\frac{7}{12}+\frac{1}{2}$. While this might lead to counter-intuitive behavior of the proposed distances in certain applications, the violated conditions are not requirements assumed by clustering algorithms.

5 Experimental Evaluation

The purpose of the experiments is to compare different methods for estimating the Levenhstein distance $\mathop{\mathrm{dist}}(A,B)$ for various real DNA sequences $A,B$ from their read sets $R_{A},R_{B}$. The methods include (i) our newly proposed distances (MES, MESS, MESSG, MESSGM, MESSGq) applicable directly on $R_{A},R_{B}$ and implemented in Java with maximum of shared code, (ii) the conventional method based on assembling estimates $\hat{A},\hat{B}$ of the original sequences $A,B$ using 5 common de-novo gene assemblers (ABySS (abyss, ), Edena (edena, ), SSAKE (ssake, ), SPAdes (spades, ) and Velvet (velvet, )) and then estimating $\mathop{\mathrm{dist}}(A,B)$ as $\mathop{\mathrm{dist}}(\hat{A},\hat{B})$, (iii) a trivial baseline method estimating $\mathop{\mathrm{dist}}(A,B)$ as $\max\{|R_{A}|,|R_{B}|\}$. Our code was single-threaded. All the 5 assembly algorithms were configured with the default parameters and the current official C++ version was used. When a result of an assembly procedure consisted of multiple contigs, we selected the two best matching contigs to represent the distance.

The evaluation criteria consist of (i) the Pearson’s correlation coefficient measuring the similarity of the distance matrix produced by the respective methods to the true distance matrix, (ii) The Fowlkes-Mallows index (fmindex, ) measuring the similarity between the tree produced by a hierarchical clustering algorithm using the true distance matrix, and the tree induced from a distance matrix estimated by the respective method, (iii) runtime needed for assembly (if applicable), distance matrix calculation, and clustering time. For hierarchical clustering, we used the UPGMA algorithm (upgma, ) and the neighbor-joining algorithm (neighjoin, ). The Fowlkes-Mallows index shows how much the resulting trees differ in structure. Both trees are first cut into $k$ clusters for $k=2,3,\ldots,n-1$. Then the clusterings are compared based on the number of common objects among each pair of clusters. By this procedure we obtain a set of values $B_{k}$ that shows how much the trees differ at various levels.

The testing data contains four datasets. They are summarized in Table 1. The “influenza” dataset555AF389115, AF389119, AY260942, AY260945, AY260949, AY260955, CY011131, CY011135, CY011143, HE584750, J02147, K00423 and outgroup AM050555. The genomes are available at http://www.ebi.ac.uk/ena/data/view/<accession>. contains 12 influenza virus genome sequences plus an outgroup sequence. The “various” dataset666AB073912, AB236320, AM050555, D13784, EU376394, FJ560719, GU076451, JN680353, JN998607, M14707, U06714, U46935, U66304, U81989, X05817, Y13051 and outgroup AY884005 contains 17 genomes of different viruses. Furthermore, we used an independent third training dataset777CY011119, CY011127, CY011140, FJ966081, AF144300, AF144300, J02057, AJ437618, FR717138, FJ869909, L00163, KJ938716, KP202150, D00664, HM590588, KM874295, $\alpha=4$, $l=40$ to tune the value of $\theta^{\prime}$ (see Sect. 3.5). All the sequences were downloaded from the ENA repository http://www.ebi.ac.uk/ena.

Those two datasets contained artificial reads that were sampled under assumptions that reads are i.i.d. and that $\alpha$ and $l$ are constant. We sampled those two datasets several times with a high range of coverage and read length888$(\alpha,l)\in\{0.1,0.3,0.5,0.7,1,1.5,2,2.5,3,4,5,7,10,15,20,30,\allowbreak 40,\allowbreak 50,\allowbreak 70,\allowbreak 100\}\times\{3,5,10,15,20,25,30,40,50,70,100,150,\allowbreak 200,\allowbreak 500\}$ to obtain a representative information about algorithms behavior based on $\alpha$ and $l$. For averaging purposes we left out the three most outlying values of $\alpha$ and $l$.

The “hepatitis” dataset contains 81 Hepatitis A segments. This time we used ART (art, ) program to similate sequencing to obtain read data for $(\alpha,l)\in\{10,30,50\}\times\{30,70,100\}$. We sampled with higher coverage and in case of our methods we down-sampled to coverage of $2$.

The “chromosomes” dataset contains $23$ regions of human genome. The data were obtained from 1000 Genomes Project (1kgenomes, ). For each chromosome we selected one $20\,\mathrm{kbp}$ long region and obtained reads that were sequenced from this region. We used the ensembl (ensembl, ) reference human genome to calculate the reference distance matrix. In this dataset the average coverage per chromosome was $4.32$. However for particular chromosomes this value varied a lot.

In the preliminary tuning experiment, the value $\theta^{\prime}=0.35$ achieved the best Pearson’s correlation coefficient on the training dataset and we carried this value over to the testing experiments.

A time limit was applied on assembly step and distance matrix calculation. Whenever an algorithm did not finish in time or assembly step failed to produce any contig, we marked the attempt as unsuccessful and did not count it towards the average.

The main experimental results are shown in Table 2. The four partitions of the table correspond to the average results on the four datasets. The Pearson’s correlation coefficient (column ‘corr.’) demonstrates dominance of the MESSG and MESSGM methods (Sects. 3.4 and 3.5), which are the most developed versions of our approach. The MESSG differs from MESSGM only by not discarding poorly matching reads. This finding is generally supported also by the Fowlkes-Mallows index (last four columns) shown for two levels of trees learned by two methods. Fig. 5 provides a more detailed insight into the Fowlkes-Mallows values graphically for all the tree levels.

Fig. 6 shows how accurate are different algorithms based on coverage. We see that our approaches produce good results for coverage around $2$. Therefore on “hepatitis” dataset we used down-sampling to $\alpha=2$. The assembly algorithms however require several times higher coverage to produce results of the same accuracy. Fig. 7 shows that our method produces good estimates for shorter reads than the assembly methods.

Columns 4-5 of Table 2 indicate that the exact variants (MES, MESS, MESSG, MESSGM) of our approach were systematically slower in terms of absolute runtime than the approaches based on sequence assembly, despite the NP-hard complexity of the latter task. However the approximated version (MESSGq, including sampling) of our approach did produce results in time, which is comparable to the assembly time. The cost for this runtime improvement was only a small decrease of accuracy. Moreover the MESSGq approach was faster than calculating the reference distance matrix on larger datasets. The numbers also show that our asymptotic complexity estimate in Sect. 4.1 is generally correct: the ratio between the time spent on calculating the distances on one hand, and the runtime of the reference method on the other hand, is approximately $\alpha^{2}$.

The “chromosomes” dataset shows how our method is dependent on the assumption that estimates of the original sequences lengths are good. While method MES provides good estimates, the accuracy drops after applying the scaling from Sect. 3.3. The coverage is not constant on all read bags and therefore the results are poor knowing that all the original sequences had the same length.

6 Conclusion

We have proposed and evaluated several variants of a method for estimating edit distances between sequences only from read-sets sampled from them. In experiments, our approach produced better estimates than a conventional approach based on first estimating the sequences themselves by applying assembly algorithms on the read-sets.

From the experiments we see that the assembly algorithms require a higher coverage $\alpha$ and read length $l$ to produce estimates comparable to our approach. This fact provides a possibility to reduce the amount of wet-lab work, which is needed to apply machine learning methods based on similarity on read sets.

Furthermore, our approach offers many directions for technical improvements. For example, one may consider a partial assembly approach, in which sets of a few (up to a constant) reads would be pre-assembled and the Monge-Elkan distance would be applied on such partial assemblies. This view would open a ‘continuous’ spectrum between our approach on one hand, and the conventional assembly-based approach, on which the optimal trade-off could be identified.

Acknowledgements.

Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme ”Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042), is greatly appreciated.