Hi everyone, I'm Auguste, PhD student at IRD Paris and Sorbonne University. Daniel and I work together on a European project and he suggested I join. I'm curating a large ITS reference database and i'm running SATIVA on it, which led me to swap its placement engine for EPA-ng few months ago, as Alexey suggested here in the discussion. The diff: github.com/Aaramis/sativa-epang/compare/master...epa-ng The benchmark against unmodified v0.9.3: github.com/Aaramis/sativa-epang-benchmark Two hunks in

sativa.py

, one in

epac/config.py

, one new module. The files holding the decision rule are byte identical to upstream. Three findings: • Daniel's S. aureus case is reproducible. Relabelling k sequences of one species with a distant genus name, unmodified SATIVA recovers 28 of 33 injected mislabels, this version 33 of 33. Every miss happens when at least two sequences share the wrong label: the strict leave-one-out removes one and the remaining copies confirm the error, a fold removes several at once. • Speed: 257 s to 15 s at 1600 sequences, 1536 s to 87 s at 5402, mostly from parallelism over queries. Agreement with 25 folds is 0.98 recall / 0.90 precision at 800 sequences, 0.85 / 0.78 at 5402. For scale, unmodified SATIVA against itself gives 0.92 / 0.90 when only its substitution model changes. • A trap for anyone doing the same: EPA-ng reads its model from the

RAxML_info

SATIVA writes, and above 500 taxa that file comes from a GTRCAT run where

alpha: 1.000000

is a placeholder. It was running GAMMA with an unfitted shape on every large alignment. Building the reference tree under GTRGAMMA fixes it for a few seconds of compute. One thing that is not in the benchmark but matters as much, and connects to the clustering discussed above: before alignment I cluster with MMseqs2 and merge clusters back together, so each alignment stays small and homogeneous enough for SATIVA. With ITS, aligning across distant taxa gives an alignment so sparse that the placements degrade and SATIVA stops proposing anything, whatever the engine. Speed is not the only reason to cluster first. My open question for the 25th: is a k-fold holdout an acceptable stand-in for the per-sequence leave-one-out, or is a batched exact one realistic at the libpll level? One sequence at a time through EPA-ng costs more than the RAxML it replaces, so batching is where the gain lives.

When you say "recovers 28 of 33" and "33 of 33" you mean your code identifies all 33 sequences as mislabelled? If so, great! 🙂

Hi Daniel, yes, all 33. They're sequences I relabelled on purpose with a distant genus. The k-fold version flags every one, unmodified SATIVA flags 28

Yes, please open issues both for actual bugs and any tips like what you say about cpu usage of steps. What Nextflow can bring is: 1. parallelization beyond what individual tools can do, 2. the capacity to run on different hardwares, including HPC (Slurm and other schedulers) and cloud, 3. software provision through container technologies and Conda (although this seems to have slipped for you), 4.

-resume

to allow restarts when something breaks or parameter choices change, 5. easier (IMO) to add steps to the pipeline. nf-core can, in addition: 1. provide standards for code and docs, 2. peer review (a double-edged sword since it can take time), 3. better chances for long-term maintenance, 4. a module repository. (I probably forgot things there.) The code you have would fit perfectly if modified into standalone programs that can be called by Nextflow. (Currently, my implementation has inline Python code. Would be much better with individual scripts in a repo under a release scheme.) And yes, if we decide to go ahead with Nextflow/nf-core it would be great to have large scale test datasets. Let's discuss further on Tuesday.

Thanks, I'd be up for that 😄 Where would you want it? My guess is a subworkflow selected by a --aligner param (mmseqs → mafft → trimal, reusing existing nf-core modules) rather than something inside the sativa module itself ? The honest blocker is that I still can't say where SATIVA stops working as a function of MSA quality. I've been using coverage, number of sequences, number of taxa, alignment width, but it's all pretty empirical and two of those already misled me. A 294-column alignment found almost no cross-kingdom contaminants while a "worse" 180-column one found about half of them, what mattered was whether all clades shared the same column frame, not how many columns survived. Occupancy didn't predict false positives either, SATIVA's mistakes on clean sequences happen where it's confident, not where the alignment is gappy. So if you know of a more principled way to score an MSA for this, I'd really like to hear it, that's the weakest part of my reasoning right now.

It could be packaged as a subworkflow on its own and added as an alternative to the current. Subworkflow and module calls basically work the same. Regarding how to automate this, I'm currently as unsure as you -- actually more -- since I haven't looked into it. However, if this requires running the Sativa subworkflow on more than one alignment to empirically score them, I suppose that could be wired in. If we decide to go ahead with this, I'd be happy to discuss.

Hi Daniel, I think I've landed on something that works, at least on my simulated benchmark. It's basically the UPP/WITCH idea from the Warnow lab, reimplemented to fit my pipeline. Instead of aligning each batch from scratch, I align one small stratified backbone once (~1200 seqs spanning every kingdom and phylum), keep the well-occupied columns as a fixed frame, then build one HMM per clade over that frame and project everything else in with nhmmer and hmmalign. The differences with UPP are that I split the backbone by taxonomy rather than by tree decomposition, and I use nhmmer to get the strand as well as the best profile. The nice side effect is that every sequence ends up in the same coordinate system, so the batches I feed to Sativa become a free slice of an already-aligned matrix instead of a separate alignment job. Re-cutting a batch costs seconds, which also makes the "run Sativa on more than one alignment and compare them" idea you mentioned a lot cheaper to do. On my simulated ITS set it detects more mislabels than the per-batch alignment, with fewer false positives on clean sequences, and the Sativa runs are faster too. Curious what you think, particularly whether this fits what you had in mind for the alternative alignment subworkflow ?

Well, the alternative alignment subworkflow was basically something I proposed when you mentioned this. The above sounds good. There's an assumption that taxonomic clades have more similar ITSs I suppose. You have thought this through much more than I and in my naivety I had been thinking of whether you had tried profile alignment. I see this is part of it, but clearly not as simple as just making a final HMM and going with that. My overall idea with the Nextflow implementation is that it should make it easier to get to a complete taxonomy-cleaning tool by including schemes like yours as well as preprocessing like raxtax, length-filtering etc. that is already included.

On the assumption you raised, it's narrower than it sounds. The taxonomy is only used to decide how to split the backbone into profiles. Each sequence is then assigned to a profile by nhmmer score, so by similarity, not by its declared label. That matters here, because a mislabelled sequence should score best against the profile of its true clade, which is exactly what we want. What I haven't tested is whether splitting the backbone by taxonomy beats splitting it by sequence clustering, the way UPP does. Worth checking. The complete tool idea makes a lot of sense to me. Having raxtax and length filtering upstream would fit well, since most of what fails at my alignment stage is short fragments that probably shouldn't reach it in the first place.