hmmer/hmmrank: multi-model domtblout multiplies rank rows and unions coordinates across profiles

[erikrikarddaniel]

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. qlen differs per model, so distinct() returns one row per model for the same (accno, profile). The left_join(domain_coords, by = c('accno', 'profile')) against the tblout side — which likewise has one row per model under the same profile — is then a cartesian product: a target hitting N models yields N × N rows instead of N, each paired with an arbitrary model's qlen and 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_to are 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.