hmmer/hmmrank: multi-model domtblout multiplies rank rows and unions coordinates across profiles
Nobody has claimed this yet.
- Dominant language
- Nextflow
- Stars
- 429
- Forks
- 1.1k
- Avg merge
- 1d 6h
- Merged PRs (30d)
- 153
Description
Description of the bug
Two things in hmmer/hmmrank, both from the coordinate columns added in #12781. The first is a real bug, the second is a test that contradicts intended behaviour.
1. Rank rows are multiplied when a domtblout holds several query models
The module derives its profile column from the file name of each *.tbl.gz / *.domtbl.gz. A single HMM file may contain several models, in which case one hmmsearch run produces one table holding rows for several query models, all sharing that one filename-derived profile. (accno, profile) is then not a unique key, and both the aggregation and the join misbehave:
domain_coords <- domains %>%
distinct(accno, profile, tlen, qlen) %>%
left_join(islands(domains, 'hmm'), by = c('accno', 'profile')) %>%
...
- Row multiplication.
qlendiffers per model, sodistinct()returns one row per model for the same(accno, profile). Theleft_join(domain_coords, by = c('accno', 'profile'))against the tblout side — which likewise has one row per model under the sameprofile— is then a cartesian product: a target hitting N models yields N × N rows instead of N, each paired with an arbitrary model'sqlenand coordinates. - Coordinates unioned across models.
islands()groups by(accno, profile)only, so domains belonging to different models are merged into one range.hmm_from/hmm_toare positions in the profile, so a union across two different profiles is not meaningful.
Neither fires without --domtblout, so output from before #12781 is unaffected.
Suggested fix. The query name is already present on both sides and distinguishes the models: column 3 of the tblout (parsed as profile_desc) and column 4 of the domtblout (currently discarded as the dummy d1). Carrying it through closes both problems at once — name d1 as the query and keep it in the transmute(), group by (accno, profile, query) in islands(), distinct(accno, profile, query, tlen, qlen), and join with by = c('accno', 'profile', 'profile_desc' = 'query').
2. The hmm_len assertion rules out an intended NA
hmmsearch reports a sequence in --tblout on the per-sequence reporting threshold, while a --domtblout row additionally needs the per-domain one. A hit whose sequence score clears while no individual domain does therefore has no domain rows, and the left_join leaves all 14 coordinate columns NA.
That is worth keeping — the empty values are the clearest available signal that a hit has no reportable domains behind it — but tests/main.nf.test asserts the opposite:
{ assert rows[1..-1].every { it.split('\t')[header.indexOf('hmm_len')].isInteger() } }
"NA".isInteger() is false, so such a row fails the test rather than being reported. Either the assertion should allow NA, or the intended behaviour should be settled the other way. Saying which in meta.yml would help downstream consumers either way.
While there: lines 78 and 145 read path(process.out.hmmrank[0][1]).linesGzip in the then block, before assert process.success runs inside assertAll. On a process failure process.out.hmmrank is empty and that line throws IndexOutOfBoundsException, hiding the real failure and its log.
Command used and terminal output
Not reproduced from a run — the module's own tests use one model per file, so the multi-model path is untested. It is reachable in nf-core/phyloplace, whose samplesheet makes the "extract one model from a multi-model HMM" field optional; without it the whole HMM file goes to hmmsearch as-is. barrnap's arc.hmm, bac.hmm, euk.hmm and mito.hmm — used in that pipeline's test data — carry 2–4 models each.
System information
n/a — logic bug, not environment-dependent.
Contributor guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Research direction
Start in the hmmer/hmmrank module at the domain_coords aggregation, islands() grouping, and tblout/domtblout parsing described in the issue. Run the module tests and inspect tests/main.nf.test, especially the hmm_len assertion and lines 78 and 145; check meta.yml for the intended NA behavior. Done means multi-model rows retain separate queries and the tests expose real process failures without rejecting intended NA values.
Written by the indexing model from the issue text.
Assessment
- Domain
- testing, tooling
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 55/100