I can confirm that the results I get with

e565a323e7

match the data on the SBDI FigShare site. Specifically, the set of archaeal 16S sequences in

correct/arc-ssu-r232-n20.correct.fna.gz

are the same as those in

sbdi-gtdb-sativa.r11rs232-1.20genomes.fna

. That said, it is unclear why some sequences are filtered out by the

e565a323e7

workflow. For example,

GB_GCA_964441565.1~CAYTWW010000004.1

is removed during initial filtering (i.e. is not in the

filtergappy/arc-ssu-r232-fwd.gapfiltered.fna

file), but 1470 bp in length on a ~200 kb contig from a species that only has a single genome. It would be nice if the SBDI pipeline produced a file indicating why sequences were removed. As is stands, this seems like a bug though perhaps it had lots of gaps in the alignment or some other property that results in it being filtered?

This might be due to our perhaps different aims. My pipeline aims to find a maximum of 20 sequences (parameter-driven) and it therefore selects up to 50 (hardcoded) sequences in the

FILTERGAPPY

module using a prioritization scheme that includes type status, isolate etc.. I suppose you're aim is to find all good sequences which makes this different of course. I think if one wants to use Sativa or a similar approach, this makes it even more important to reduce the number of sequences by clustering since some species have an enormous number of sequences.

Our aims certainly do differ in this respect. However, I'm still confused about some sequences that don't make it to the FILTERGAPPY step. Many of these are the only 16S sequence for a species and they appear fine on the surface. Perhaps SATIVA will find them to actually be contaminated, but any criterion used to prefilter them isn't obvious to me.

You're absolutely right. I'll trace the code to see where I lose those sequences. It's sort of vacation season here now, so I'm not in front of the computer many hours every day. But, I'm curious so I'll get this done soon