Parameters for matching ITSx exactly?

[sujaikumar]

Dear Aaron

Congratulations on re-implementing ITSx - you're right - it IS a bottleneck and having a fast rust version that also works on long reads is terrific!

I was just trying this out as a replacement for ITSx in the lotus3 pipeline.

Here is the original command in ITSx (installed using conda install -c bioconda -c conda-forge lotus3)

ITSx -i otu_seeds.fna -o otu_seeds.fna.itsX --cpu 36 --multi_thread T \
  --heuristics T -t F --silent T --fasta T --nhmmer F --save_regions ITS2 --partial 0 \
  --hmmBin /home/ubuntu/miniforge3/envs/lotus3/bin/hmmscan --preserve T

And this is the corresponding itsxrust command I came up with (hmmscan was already in the path, and ITSx -t F means look up the F.hmm model in the ITSx_db/HMMs installation folder) - installed using conda install -c bioconda -c conda-forge itsxrust

itsxrust extract --input otu_seeds.fna --output otu_seeds.fna.itsxrust --hmmer-cpu 36 \
  --hmm /home/ubuntu/miniforge3/envs/lotus3/bin/ITSx_db/HMMs/F.hmm --region its2

However, I am getting subtly different results:

• ITSx: 2106 ITS2 regions
• itsxrust: 1406 ITS2 regions

I am assuming that the difference is because of default values of other parameters which I did not set in either ITSx or itsxrust.

Could you help identify the corresponding parameters in itsxrust?

If it helps, the ITSx help says:

ITSx -- Identifies ITS sequences and extracts the ITS region
by Johan Bengtsson-Palme et al., University of Gothenburg
Version: 1.1.3
-----------------------------------------------------------------
Usage: ITSx -i <input file> -o <output file>
Options:
-i {file} : DNA FASTA input file to investigate
-o {file} : Base for the names of output file(s)
-p {directory} : A path to a directory of HMM-profile collections representing ITS conserved regions, default is in the same directory as ITSx itself
--stdin {T or F} : Use input from standard input instead of an input file, off (F) by default
--date {T or F} : Adds a date and time stamp to the output directory, off (F) by default
--reset {T or F} : Re-creates the HMM-database before ITSx is run, off (F) by default

Sequence selection options:
-t {character code} : Profile set to use for the search, see the User's Guide (comma-separated), default is all
-E {value} : Domain E-value cutoff for a sequence to be included in the output, default = 1e-5
-S {value} : Domain score cutoff for a sequence to be included in the output, default = 0
-N {value} : The minimal number of domains that must match a sequence before it is included, default = 2
--selection_priority {score, sum, domains, eval} : Selects what will be of highest priority when determining the origin of the sequence, default is score
--search_eval {value} : The E-value cutoff used in the HMMER search, high numbers may slow down the process, cannot be used with the --search_score option, default is to use score cutoff, not E-value
--search_score {value} : The score cutoff used in the HMMER search, low numbers may slow down the process, cannot be used with the --search_eval option, default = 0
--allow_single_domain {e-value,score or F} : Allow inclusion of sequences that only find a single domain, given that they meet the given E-value and score thresholds, on with parameters 1e-9,0 by default
--allow_reorder {T or F} : Allows profiles to be in the wrong order on extracted sequences, off (F) by default
--complement {T or F} : Checks both DNA strands against the database, creating reverse complements, on (T) by default
--cpu {value} : the number of CPU threads to use, default is 1
--multi_thread {T or F} : Multi-thread the HMMER-search, on (T) if number of CPUs (--cpu option > 1), else off (F) by default
--heuristics {T or F} : Selects whether to use HMMER's heuristic filtering, off (F) by default
--nhmmer {T or F} : Selects whether to use nhmmer instead of hmmsearch for HMMER searches, off (F) by default

Output options:
--summary {T or F} : Summary of results output, on (T) by default
--graphical {T or F} : 'Graphical' output, on (T) by default
--fasta {T or F} : FASTA-format output of extracted ITS sequences, on (T) by default
--preserve {T or F} : Preserve sequence headers in input file instead of printing out ITSx headers, off (F) by default
--save_regions {SSU,ITS1,5.8S,ITS2,LSU,all,none} : A comma separated list of regions to output separate FASTA files for, 'ITS1,ITS2' by default
--anchor {integer or HMM} : Saves an additional number of bases before and after each extracted region. If set to 'HMM' all bases matching the corresponding HMM will be output, default = 0
--require_anchor {integer or HMM} : Requires the complete anchor to found in order to be included in the output sequences (see --anchor above). Cannot be used together with the --anchor option, default = 0
--only_full {T or F} : If true, output is limited to full-length regions, off (F) by default
--partial {integer} : Saves additional FASTA-files for full and partial ITS sequences longer than the specified cutoff, default = 0 (off)
--concat {T or F} : Saves a FASTA-file with concatenated ITS sequences (with 5.8S removed), off (F) by default
--minlen {integer} : Minimum length the ITS regions must be to be outputted in the concatenated file (see above), default = 0
--positions {T or F} : Table format output containing the positions ITS sequences were found in, on (T) by default
--table {T or F} : Table format output of sequences containing probable ITS sequences, off (F) by default
--not_found {T or F} : Saves a list of non-found entries, on (T) by default
--detailed_results {T or F} : Saves a tab-separated list of all results, off (F) by default
--truncate {T or F} : Truncates the FASTA output to only contain the actual ITS sequences found, on (T) by default
--silent {T or F} : Supresses printing progress info to stderr, off (F) by default
--graph_scale {value} : Sets the scale of the graph output, if value is zero, a percentage view is shown, default = 0
--save_raw {T or F} : Saves all raw data for searches etc. instead of removing it on finish, off (F) by default
--temp {directory} : Custom directory to put the temporary files in

-h : displays this help message
--help : displays this help message
--bugs : displays the bug fixes and known bugs in this version of ITSx
--license : displays licensing information

[ayobi]

Hi Sujai,

Thanks for the kind words!! and thanks especially for the careful side-by-side commands..that makes diagnosis much easier.

A few things that are not the cause, just to clear the it up: ITSxRust's HMMER E-value default (--inc-e 1e-5) matches ITSx's -E 1e-5, and ITSxRust searches both strands by default with reverse-complementation when the best chain is on the minus strand. So those line up..

The most likely contributors might be these ideas below::

1. Two-anchor requirement vs read-edge fallback. ITSxRust delimits ITS2 using two HMM anchors (5.8S 3′end + LSU 5′start) and requires both to be found. ITSx in your configuration has --allow_single_domain 1e-9,0 ON by default and falls back to the read boundary as a delimiter when a flank isn't found..so a read missing the LSU flank can still produce an ITS2 extract from ITSx but not from ITSxRust. This is a documented design trade-off (paper §4.1) and it's likely the dominant contributor to the ~700-sequence gap. For ITS1 it accounted for an ~18% extraction gap in our ONT benchmark; the equivalent effect is plausibly responsible here.

2. hmmscan vs nhmmer. Your ITSx command uses --nhmmer F, so ITSx calls hmmscan on F.hmm. ITSxRust uses nhmmer. Same profiles, different search algorithms., small per-hit sensitivity differences are possible, particularly for marginal hits at OTU-seed edges. Likely a second-order effect, but worth flagging.

3. maybe try the ONT preset. Even though lotus3 OTU seeds are Illumina-derived, --preset ont relaxes thresholds in ways that may help on short representative sequences: --inc-e 1e-3, --max-anchor-evalue 1e-3, --min-anchor-score 15, --min-its2 30 --max-its2 2500. Cheap to run as a bracket. :)

The cleanest way to pinpoint which mechanism is dropping the reads is to re-run with the built-in diagnostics:

itsxrust extract --input otu_seeds.fna --output otu_seeds.fna.itsxrust \
  --hmmer-cpu 36 --hmm /home/ubuntu/miniforge3/envs/lotus3/bin/ITSx_db/HMMs/F.hmm \
  --region its2 \
  --qc-json itsxrust_qc.json \
  --anchors-tsv itsxrust_anchors.tsv \
  --explain 20

The QC JSON gives a structured skip-reason breakdown.. :) how many reads were missing the 5.8S anchor, how many missing LSU, how many had anchors but failed length or chain constraints which is nice imo.. the anchors TSV gives per-read coordinates and confidence labels for reads that did extract. If you can share the QC JSON (and ideally the same for the --preset ont run), I can tell you straight away which mechanism accounts for the gap.

Big picture though: for ITSxRust to be a friction-free drop-in for ITSx in pipelines like lotus3, supporting a read-edge fallback option for the two-anchor case is probably what's needed. It's on the roadmap as a user-configurable option; your data point would help prioritise it. If the QC JSON confirms that missing-anchor reads are the bulk of the gap, I'd happily push that up the queue :))

Thanks again for trying it!! neat to see ITSxRust evaluated for lotus3!

Saludos,
Aaron

[sujaikumar]

Hi Aaron. Thank you for the SUPER clear explanation of what might be going on!

Attached are the outputs of:

itsxrust extract --input otu_seeds.uniq.fna --output otu_seeds.uniq.fna.itsxrust \
  --hmmer-cpu 36 --hmm /home/ubuntu/miniforge3/envs/lotus3/bin/ITSx_db/HMMs/F.hmm \
  --region its2 \
  --qc-json itsxrust_qc.json \
  --anchors-tsv itsxrust_anchors.tsv \
  --explain 20

itsxrust extract --input otu_seeds.uniq.fna --output ont.otu_seeds.uniq.fna.itsxrust \
  --hmmer-cpu 36 --hmm /home/ubuntu/miniforge3/envs/lotus3/bin/ITSx_db/HMMs/F.hmm \
  --region its2 \
  --preset ont \
  --qc-json ont.itsxrust_qc.json \
  --anchors-tsv ont.itsxrust_anchors.tsv \
  --explain 20

(I cleaned the input files as they had some duplicate IDs and sequences from the lotus pipeline. The final output of the first command now has 1407 ITS2 regions, and the --preset ont one has 1625)

The QC JSONs do confirm that the majority are due to missing gaps - good to know.

itsxrust_qc.json
itsxrust_anchors.tsv
ont.itsxrust_qc.json
ont.itsxrust_anchors.tsv

However the anchors.tsv files are both empty. Not sure if it's because of my installation...

Hope this helps!

• Sujai

[ayobi]

Hi Sujai,

Really useful info, thanks for attaching these. I'm glad the QC helped because I saw two things come out of it.

The empty anchors.tsv is a bug on my end sadly.. looking at the QC JSON, every kept read in your data is partial (0 confident, 0 ambiguous, all 1407/1625 partial). That makes sense for ITS2 amplicons amplified with primers inside 5.8S: the SSU and 5.8S 5′start anchors are never in the read, so the full 4-anchor chain can't form and every extraction goes through the 2-anchor fallback path. My writer is hooked into the full-chain code path but not the partial one, so it produces just the header when 100% of reads are partial. Small fix, will be in the next patch release when I get the time :)

On the gap ...the ONT preset's gain came almost entirely from loosening --inc-e (1e-5 → 1e-3), which let nhmmer find hits in 217 more reads. The per-anchor thresholds barely moved the needle haha so missing_anchors only dropped from 1125 to 1124. those 1124 reads aren't being filtered for being below threshold; nhmmer just isn't finding one of the two anchors in the read at all. Threshold tuning won't recover them, but read-edge fallback def would.

With ~1124 reads recoverable that way, ITSxRust would land around 2750/2947, comfortably above ITSx's 2106. I'm going to build it as an opt-in flag (--read-edge-fallback), defaulting off, with edge-called boundaries tagged in the output. Hopefully I can get to it soon!

The 198 no_hmm_hits reads in the ONT run are a separate thing: nhmmer found nothing at all in those. Might be the hmmscan-vs-nhmmer sensitivity difference, but smaller effect than the structural gap.

Thanks again for running this and the feedback!