Sintax confidence

[kabilov]

Hi Robert!

My question concerns the Sintax in Usearch11 (usearch11.0.667_i86linux64).
The database I'm using (16S V3-V4) contains six identical sequences, but five of them are Archaea, and one is a Bacteria. The latter is probably an artifact, but that's not the point.

d:Archaea,p:Methanobacteriota,c:Methanobacteria,o:Methanobacteriales,f:Methanothermobacteraceae,g:Methanothermobacter,s:Methanothermobacter_sp022014235;
d:Archaea,p:Methanobacteriota,c:Methanobacteria,o:Methanobacteriales,f:Methanothermobacteraceae,g:Methanothermobacter,s:Methanothermobacter_sp000828575;
d:Archaea,p:Methanobacteriota,c:Methanobacteria,o:Methanobacteriales,f:Methanothermobacteraceae,g:Methanothermobacter,s:Methanothermobacter_thermautotrophicus;
d:Archaea,p:Methanobacteriota,c:Methanobacteria,o:Methanobacteriales,f:Methanothermobacteraceae,g:Methanothermobacter,s:Methanothermobacter_wolfei;
d:Archaea,p:Methanobacteriota,c:Methanobacteria,o:Methanobacteriales,f:Methanothermobacteraceae,g:Methanothermobacter,s:Methanothermobacter_thermautotrophicus_A;
d:Bacteria,p:Pseudomonadota,c:Gammaproteobacteria,o:Pseudomonadales,f:Porticoccaceae,g:CAYXXZ01,s:CAYXXZ01_sp050175615;
GCGCGAAACCTCCGCAATGCACGCAAGTGCGACGGGGGAACCCCAAGTGCCACTCTTAACGGGGTGGCTTTTCAGAAGTGTAAAAAGCTTCTGGAATAAGGGCTGGGCAAGACCGGTGCCAGCCGCCGCGGTAACACCGGCAGCTCAAGTGGTAGCCGCTTTTATTGGGCCTAAAGCGTCCGTAGCCGGTCTGATAAGTCTCTGGTGAAATCCCACAGCTTAACTGTGGGAATTGCTGGAGATACTATCATGACTCGAGGTCGGGAGAGGCTGGAGGTACTCCCAGGGTAGGGGTGAAATCCTGTAATCCTGGGAGGACCACCTGTGGCGAAGGCGTCCAGCTGGAACGAACCTGACGGTGAGGGACGAAAGCCAGGGGCGCGAACCGG

I expected the sintax analysis to yield a confidence level of 83.3% (5/6) across all levels (except species). But the actual result was different:
d:Archaea(0.8700) p:Methanobacteriota(0.7569) c:Methanobacteria(0.6585) o:Methanobacteriales(0.5729) f:Methanothermobacteraceae(0.4984) g:Methanothermobacter(0.4336) s:Methanothermobacter_thermautotrophicus_A(0.0781)
0.87 != 5/6 = 0.83

d = (0.87)^1 = 0.87
p = (0.87)^2 = 0.76
...
g = (0.87)^6 = 0.43

The result was roughly the same when the database consisted only of these six sequences.
d:Archaea(0.8800) p:Methanobacteriota(0.7744) c:Methanobacteria(0.6815) o:Methanobacteriales(0.5997) f:Methanothermobacteraceae(0.5277) g:Methanothermobacter(0.4644) s:Methanothermobacter_sp022014235(0.1022)
d = (0.88)^1 = 0.88
p = (0.88)^2 = 0.77
...
g = (0.88)^6 = 0.46

If I removed the bacterial sequence, the result became almost as expected (it's unclear why s = 1/4 and not 1/5).
d:Archaea(1.0000) p:Methanobacteriota(1.0000) c:Methanobacteria(1.0000) o:Methanobacteriales(1.0000) f:Methanothermobacteraceae(1.0000) g:Methanothermobacter(1.0000) s:Methanothermobacter_sp022014235(0.2500)

Could you comment on these inconsistencies?

Best wishes,
Marsel

[rcedgar]

Confidence is calculated by bootstrapping, as follows. Selects a subset (call it S) of k-mers from the query. The database is sorted in decreasing order of U = (number of unique k-mers in common between S and the database sequence). Take the top hit. Repeat 100 times, giving 100 top hits. For each rank, find the most abundant name among these 100 hits, its confidence is the fraction of times it occurs. Often, there are ties for the top hit, if so the tie is broken arbitrarily. Ties can happen without identical database sequences; they are quite common.

In your case, you should not expect a predictable or consistent result because of the randomness associated with selecting the subset and in tie-breaking.

You would get somewhat more consistent results by ensuring there are no identical sequences in the database; in principle they should be combined into one sequence which is annotated with the common taxonomy, deleting lower ranks where the taxon names are different. But this is probably not worth the trouble unless duplicates are common.

[kabilov]

Robert,

Thank you for the quick reply. Perhaps I did not explain my question clearly enough.

I am interested in the following scenario: the SINTAX database contains six completely identical sequences that differ only in their taxonomic annotations. As I understand it, in such a case the confidence score should no longer depend on k-mer frequencies, but rather on the probability of assigning the query to a particular taxon.

In my example, five of the six sequences share the same taxonomy down to the genus level, whereas the sixth differs at every taxonomic rank because it belongs to Archaea rather than Bacteria. Therefore, I would expect the confidence value at each taxonomic rank to be approximately 5/6 (≈0.83). However, the reported confidence values differ substantially from this expectation.

Could you please explain how SINTAX calculates confidence scores in this situation?

[rcedgar]

It's been a long time since I worked on that code, I did my best to explain what happens as best I remember it. To me, the result looks reasonable even if I'm not sure exactly where the numbers come from. Note that the confidence estimate is a heuristic which is at best a rough guideline. If you want to dig into this further, I suggest running with -threads 1 and adding some printf statements to the code or step through with a debugger. The relevant function is void SintaxSearcher::Classify() in sintaxsearcher.cpp.