peterjc
commented
1 year ago Partly answering my own question: Looking at the code, it seems for its UNOISE implementation VSEARCH uses the Needleman-Wunsch aligner in https://github.com/torognes/vsearch/blob/master/src/align.cc and https://github.com/torognes/vsearch/blob/master/src/searchcore.cc as follows:
if (opt_cluster_unoise)
{
int d = hit->mismatches;
double skew = 1.0 * si->qsize / db_getabundance(hit->target);
double beta = 1.0 / pow(2, 1.0 * opt_unoise_alpha * d + 1);
if (skew <= beta || d == 0)
{
/* accepted */
...Quoting the man-page https://github.com/torognes/vsearch/blob/master/man/vsearch.1 explains much of the scoring:
When comparing sequences during chimera detection, dereplication, searching and clustering, T and U are considered identical, regardless of their case. When aligning sequences, identical symbols will receive a positive match score (default +2). If two symbols are not identical, their alignment result in a negative mismatch score (default -4). Aligning a pair of symbols where at least one of them is an ambiguous symbol (BDHKMNRSVWY) will always result in a score of zero. Alignment of two identical ambiguous symbols (for example, R vs R) also receives a score of zero. When computing the amount of similarity by counting matches and mismatches after alignment, ambiguous nucleotide symbols will count as matching to other symbols if they have at least one of the nucleotides (ACGTU) they may represent in common. For example: W will match A and T, but also any of MRVHDN. When showing alignments (for example with the --alnout option) matches involving ambiguous symbols will be shown with a plus character (+) between them while exact matches between non-ambiguous symbols will be shown with a vertical bar character (|).
This is consistent with https://github.com/torognes/vsearch/blob/master/src/align.cc which also gives the gap penalties:
default (interior) scores:
match: +2
mismatch: -4
gap open: 20
gap extend: 2i.e. This looks like it does NW alignment, and then considers the mismatch count of the best alignment(s) for the UNOISE distance.
torognes
The Edgar 2016 paper says UNOISE2 used the Levenshtein distance: https://www.biorxiv.org/content/10.1101/081257v1.full
Comment https://github.com/torognes/vsearch/pull/283#issuecomment-360850868 on the original pull request adding UNOISE functionality to VSEARCH suggests from empirical testing that the USEARCH implementation of UNOISE2 and UNOISE3 was not using Levenshtein distance, but counting substitutions.
Do you recall how you came to that conclusion @davidealbanese?
I've been trying to validate my own UNOISE implementation by cross checking against VSEARCH (and USEARCH), from which I have this minimal test case.
Consider the following pair of test cases, where I am expecting in both cases denoising to merge the low abundance variant with the high abundance centroid.
In the above example with ratio at exactly 1:128, the variant gets merged into the more common sequence. However, if we increase the variant's abundance by one it does not:
This ratio of 1:128 is a threshold in the UNOISE beta function with dist=3, alpha=2.
This is behaving as though VSEARCH thinks the distance between these two sequences is 3, giving a beta function value 1 / 2^(alphadist+1) = 1 / 2^(23+1) = 1 / 2^7 = 1 / 128.
However while this could be viewed as a difference of three substitutions, I think the Levenshtein distance is only 2, two single base inserts - see for example this EMBOSS needle alignment:
This does match the UNOISE2 behaviour, ... <-- Edit, I need to double check this
Note that in general the two sequences need not be the same length, so I am unclear what exactly is being used as the distance rule.