Skip to content

hmmsearch -T prefilter is global (shared HMM module) — loosens VOG/DRAM-v and sulfur/metals; make it per-DB #13

Open

Description

@tpall

Summary

PR #12 switched the HMM prefilter from -E ${e_value} to a fixed bit-score floor -T 10 --domT 10. That change lives in the shared HMM_SEARCH process (modules/local/annotate/hmmsearch.nf), so it applies to every HMM search — kofam, sulfur, vog, canthyd, fegenie, metals, camper-hmm — not only the pooled kofam/sulfur searches the PR targeted. The per-DB e_value inputs (e.g. vog_e_value = 1e-10) are now silently ignored.

-T is required for the pooled path (bit scores are search-space/Z-independent, so pooled == per-genome), so a global revert to -E isn't the answer — the prefilter needs to be per-DB.

Why it matters (benign for some DBs, not others)

Whether the swap is harmless depends on the downstream filtering in bin/hmm_parser.py:

  • kofam (plus canthyd via A/B rank, fegenie via its cutoffs) carry per-profile bit-score thresholds, so hmm_parser keeps score >= threshold and the -T 10 prefilter is a near no-op — the "+0.08%" the PR reported. ✅
  • vog, sulfur, metals pass the dummy sheet (bin/assets/forms/distill_sheets/dummy.tsv, contents NULL), so hmm_parser finds no threshold and falls through to only perc_cov >= 0.35. Their real significance gate was the hmmsearch prefilter itself. For VOG that was -E 1e-10 (stringent); -T 10 (bit score ≥ 10) is far looser. ⚠️

Measured impact (DRAM-v mode)

Ran DRAM-v on a 161-contig viral voTU catalog (geNomad gene input, catalog mode, --max_auxiliary_score 3, --use_dramv --use_kofam --use_merops) on two refs differing only by this prefilter line:

metric -E 1e-10 -T 10 change
VOG-annotated genes 2871 5872 +105% (doubles)
V-flag genes 999 1737 +74%
auxiliary_score histogram {1: 6848, 4: 1416, 5: 1212} identical 0
AMG calls (M & aux ≤ 3) 294 294 0

Raw VOG annotations and the V flag roughly double, but the AMG endpoint is unchanged — auxiliary_score is geNomad/VirSorter-derived (independent of VOG), and the AMG gate keys off amg_flags (M, lacking A/P/T/N) + auxiliary_score, none of which depend on VOG.

Net: AMG-identification output is safe, but any downstream use of raw VOG annotations / V-flag counts / "% viral genes" is now ~2× inflated with low-confidence hits.

Proposed fix

Make the prefilter per-DB instead of hardcoded in the shared module. HMM_SEARCH already reads def args = task.ext.args ?: "" but never uses it in the command — wire the prefilter through that (or an explicit process input), so:

  • pooled kofam/sulfur searches use -T 10 --domT 10 (needed for Z-independent pooling; provably benign there),
  • per-genome vog/sulfur/metals keep -E ${e_value} (e.g. VOG -E 1e-10).

This preserves PR #12's pooling-correctness win while restoring the intended significance filter for the threshold-less DBs.

Refs: modules/local/annotate/hmmsearch.nf, bin/hmm_parser.py, docs/dev/search-bundling.md.


🤖 Filed with Claude Code after an -E-vs--T DRAM-v validation run.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions