From 03f0f5d8c6cfe1f04dc40cb52f23d99e7dea7803 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Taavi=20P=C3=A4ll?= Date: Sun, 28 Jun 2026 09:36:53 +0300 Subject: [PATCH 1/3] experiment(perf): HMM pooling via -T bit-score prefilter (Z-independent) On branch perf/hmm-pool-bitscore (experiment, NOT for PR #10). Re-enables kofam/sulfur pooling + chunking (un-reverts 3457a801) and switches hmmsearch from the e-value prefilter (-E, scales with search-space size Z) to a bit-score prefilter (-T 10 --domT 10). Bit scores are Z-independent, so pooled multi-genome hmmsearch reports the SAME hits as per-genome; hmm_parser still applies per-profile bit-score thresholds. Floor 10 < lowest kofam threshold (15.07), so no real hit is pre-cut. Goal: validate per-genome(-T) == pooled(-T) on extraves (proves -T removes the Z-dependence that made -E pooling churn ~5% of kofam calls). NB changes results vs -E, so a future pipeline-wide decision, not mid-campaign. Co-Authored-By: Claude Opus 4.8 (1M context) --- modules/local/annotate/hmmsearch.nf | 8 ++- nextflow_schema.json | 5 +- subworkflows/local/db_search.nf | 102 ++++++++++++++++++++-------- 3 files changed, 82 insertions(+), 33 deletions(-) diff --git a/modules/local/annotate/hmmsearch.nf b/modules/local/annotate/hmmsearch.nf index 4e20eb62..fb6c7d51 100644 --- a/modules/local/annotate/hmmsearch.nf +++ b/modules/local/annotate/hmmsearch.nf @@ -22,8 +22,14 @@ process HMM_SEARCH { def ec_flag = ec_from_info ? "--ec_from_info" : "" """ + # Bit-score prefilter (-T/--domT) instead of the e-value prefilter (-E): bit scores + # are independent of search-space size (Z), so a pooled multi-genome search reports + # the SAME hits as per-genome (hmm_parser then applies the per-profile bit-score + # thresholds). Floor 10 sits below the lowest kofam threshold (15.07) so nothing a + # real threshold would keep is pre-cut. NB: this is the perf/hmm-pool-bitscore + # experiment — it changes results vs the -E ${e_value} prefilter. hmmsearch \\ - -E ${e_value} \\ + -T 10 --domT 10 \\ --domtblout ${input_fasta}_hmmsearch.out \\ --cpu ${task.cpus} \\ ${database_loc}/*.hmm \\ diff --git a/nextflow_schema.json b/nextflow_schema.json index bbdd81b7..f3766b21 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -553,9 +553,10 @@ "hidden": true }, "search_chunk_size": { - "type": "integer", + "type": ["integer", "string"], + "pattern": "^[0-9]+$", "default": 0, - "description": "Genes per pooled search chunk when pool_searches is set (e.g. 100000). 0 = one chunk per genome (current behaviour). Scaffolding — not wired yet.", + "description": "Genes per pooled HMM search chunk when pool_searches is set (e.g. 20000). 0 = single pooled search task. Splits the pooled FASTA for parallel hmmsearch; see docs/dev/search-bundling.md.", "hidden": true }, "camper_e_value": { diff --git a/subworkflows/local/db_search.nf b/subworkflows/local/db_search.nf index 637dc409..376091fb 100644 --- a/subworkflows/local/db_search.nf +++ b/subworkflows/local/db_search.nf @@ -46,6 +46,10 @@ include { HMM_SEARCH as HMM_SEARCH_CANTHYD } from "../../modules/lo include { HMM_SEARCH as HMM_SEARCH_SULFUR } from "../../modules/local/annotate/hmmsearch.nf" include { HMM_SEARCH as HMM_SEARCH_FEGENIE } from "../../modules/local/annotate/hmmsearch.nf" include { HMM_SEARCH as HMM_SEARCH_METALS } from "../../modules/local/annotate/hmmsearch.nf" +include { HMM_SEARCH as HMM_SEARCH_KOFAM_POOLED } from "../../modules/local/annotate/hmmsearch.nf" +include { HMM_SEARCH as HMM_SEARCH_SULFUR_POOLED } from "../../modules/local/annotate/hmmsearch.nf" +include { SPLIT_POOLED_HITS as SPLIT_POOLED_KOFAM } from "../../modules/local/annotate/split_pooled_hits.nf" +include { SPLIT_POOLED_HITS as SPLIT_POOLED_SULFUR } from "../../modules/local/annotate/split_pooled_hits.nf" include { CONCAT_HMM_HITS as CONCAT_HMM_HITS_KOFAM } from "../../modules/local/annotate/concat_hmm_hits.nf" include { CONCAT_HMM_HITS as CONCAT_HMM_HITS_VOG } from "../../modules/local/annotate/concat_hmm_hits.nf" @@ -147,32 +151,33 @@ workflow DB_SEARCH { def formattedOutputchannels = channel.of() def dbcanOutputChannels = channel.of() - // Pooled search bundling currently covers the merops/viral/methyl mmseqs searches. - // The other mmseqs DBs still reference the per-genome index (skipped when pooling), - // so guard against enabling them together with --pool_searches. HMM searches are - // unaffected (they always run per-genome). Generalizing the rest is mechanical — - // see docs/dev/search-bundling.md. + // Phase-2 search bundling (docs/dev/search-bundling.md): build the pooled, + // genome-prefixed protein FASTA + gene-locs ONCE when --pool_searches — shared by + // every pooled mmseqs AND hmm search. Gene ids are made globally unique by + // PREFIX_GENES_FOR_POOL (___; megahit k141_* scaffold names collide + // across bins); SPLIT_POOLED_HITS strips the prefix per genome. Pooling currently + // covers the mmseqs merops/viral/methyl and hmm kofam/sulfur searches; the other + // mmseqs DBs still use the per-genome index (skipped under pooling), so guard them. if (params.pool_searches && (use_kegg || use_pfam || use_camper || use_canthyd || use_uniref)) { - error("--pool_searches currently supports the merops, viral and methyl mmseqs searches only. Disable kegg/pfam/camper/canthyd/uniref or drop --pool_searches.") + error("--pool_searches currently supports the merops/viral/methyl mmseqs searches and the kofam/sulfur hmm searches only. Disable kegg/pfam/camper/canthyd/uniref or drop --pool_searches.") + } + if (params.pool_searches) { + PREFIX_GENES_FOR_POOL( ch_called_proteins.join(ch_gene_locs) ) + ch_pooled_faa = PREFIX_GENES_FOR_POOL.out.faa + .collectFile(name: 'pooled_genes.faa') + .map { faa -> tuple('pooled', faa) } + ch_pooled_locs = PREFIX_GENES_FOR_POOL.out.locs + .collectFile(name: 'pooled_gene_locs.tsv', keepHeader: true) + .map { locs -> tuple('pooled', locs) } + ch_pooled_hmm_in = ch_pooled_faa.join(ch_pooled_locs) // [ 'pooled', faa, gene_locs ] } // Here we will create mmseqs2 index files for each of the inputs if we are going to do a mmseqs2 database if (DB_channel_SETUP.out.index_mmseqs) { if (params.pool_searches) { - // Phase-2 search bundling (docs/dev/search-bundling.md): index the POOLED - // proteins ONCE and reuse the pooled query DB for every mmseqs search. - // Gene ids are made globally unique by PREFIX_GENES_FOR_POOL - // (___; megahit k141_* scaffold names collide across bins); - // SPLIT_POOLED_HITS strips the prefix per genome. The per-genome index is - // skipped entirely. ch_pooled_search_in = [ 'pooled', index_db, gene_locs ]. - PREFIX_GENES_FOR_POOL( ch_called_proteins.join(ch_gene_locs) ) - ch_pooled_faa = PREFIX_GENES_FOR_POOL.out.faa - .collectFile(name: 'pooled_genes.faa') - .map { faa -> tuple('pooled', faa) } + // index the pooled proteins once; reuse for every pooled mmseqs search. + // ch_pooled_search_in = [ 'pooled', index_db, gene_locs ]. MMSEQS_INDEX_POOLED( ch_pooled_faa ) - ch_pooled_locs = PREFIX_GENES_FOR_POOL.out.locs - .collectFile(name: 'pooled_gene_locs.tsv', keepHeader: true) - .map { locs -> tuple('pooled', locs) } ch_pooled_search_in = MMSEQS_INDEX_POOLED.out.mmseqs_index_out.join(ch_pooled_locs) } else { @@ -195,7 +200,35 @@ workflow DB_SEARCH { } // KOFAM annotation if (use_kofam) { - if (params.kofam_chunk_size && params.kofam_chunk_size > 0) { + if (params.pool_searches) { + // Pooled HMM. One hmmsearch over the pooled prefixed proteins is a single + // un-parallelised task; with search_chunk_size > 0 the pooled FASTA is split + // into chunks searched in parallel (restoring fan-out while keeping the task + // count O(#chunks), not O(#bins)), concatenated, then split per genome. + // hmm_parser filters on per-profile bit-score thresholds (DB-size independent); + // hmmsearch's -E ${kofam_e_value} prefilter scales with the searched DB size — + // validated on extraves. + def search_chunk_n = (params.search_chunk_size as Integer) // CLI passes a string + if (search_chunk_n > 0) { + ch_kofam_pooled_in = ch_pooled_faa + .splitFasta(by: search_chunk_n, file: true, elem: 1) + .map { pool, chunk -> tuple("pooled___${chunk.baseName}", chunk) } + .combine( ch_pooled_locs.map { it[1] } ) + } + else { + ch_kofam_pooled_in = ch_pooled_hmm_in + } + HMM_SEARCH_KOFAM_POOLED( ch_kofam_pooled_in, params.kofam_e_value, DB_channel_SETUP.out.ch_kofam_db, ch_kofam_list, true, kofam_name ) + ch_kofam_pooled_csv = HMM_SEARCH_KOFAM_POOLED.out.formatted_hits + .map { id, csv -> csv } + .collectFile(name: 'pooled_kofam_hits.csv', keepHeader: true) + .map { csv -> tuple('pooled', csv) } + SPLIT_POOLED_KOFAM( ch_kofam_pooled_csv, kofam_name ) + ch_kofam_formatted = SPLIT_POOLED_KOFAM.out.per_genome_hits + .flatten() + .map { csv -> tuple(csv.name.toString().split('___')[0], csv) } + } + else if (params.kofam_chunk_size && params.kofam_chunk_size > 0) { ch_kofam_chunks = ch_called_proteins .splitFasta(by: params.kofam_chunk_size, file: true, elem: 1) .combine(ch_gene_locs, by: 0) @@ -337,16 +370,25 @@ workflow DB_SEARCH { } // Sulfur annotation if (use_sulfur) { - ch_combined_proteins_locs = ch_called_proteins.join(ch_gene_locs) - HMM_SEARCH_SULFUR ( - ch_combined_proteins_locs, - params.sulfur_e_value, - DB_channel_SETUP.out.ch_sulfur_db, - default_sheet, - false, - sulfur_name - ) - ch_sulfur_formatted = HMM_SEARCH_SULFUR.out.formatted_hits + if (params.pool_searches) { + HMM_SEARCH_SULFUR_POOLED( ch_pooled_hmm_in, params.sulfur_e_value, DB_channel_SETUP.out.ch_sulfur_db, default_sheet, false, sulfur_name ) + SPLIT_POOLED_SULFUR( HMM_SEARCH_SULFUR_POOLED.out.formatted_hits, sulfur_name ) + ch_sulfur_formatted = SPLIT_POOLED_SULFUR.out.per_genome_hits + .flatten() + .map { csv -> tuple(csv.name.toString().split('___')[0], csv) } + } + else { + ch_combined_proteins_locs = ch_called_proteins.join(ch_gene_locs) + HMM_SEARCH_SULFUR ( + ch_combined_proteins_locs, + params.sulfur_e_value, + DB_channel_SETUP.out.ch_sulfur_db, + default_sheet, + false, + sulfur_name + ) + ch_sulfur_formatted = HMM_SEARCH_SULFUR.out.formatted_hits + } formattedOutputchannels = formattedOutputchannels.mix(ch_sulfur_formatted) } // MEROPS annotation From c9665817f1e9ef47e2fd235ee997ccb5827ce4e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Taavi=20P=C3=A4ll?= Date: Sun, 28 Jun 2026 12:02:04 +0300 Subject: [PATCH 2/3] docs(perf): -T bit-score prefilter makes HMM pooling byte-identical (validated) Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/dev/search-bundling.md | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/dev/search-bundling.md b/docs/dev/search-bundling.md index 4fa265fb..7bd88f72 100644 --- a/docs/dev/search-bundling.md +++ b/docs/dev/search-bundling.md @@ -35,6 +35,17 @@ at 56 bins: per-genome = 336 tasks (56 index + 168 search + 112 SQL) → pooled (3823 dropped / 3799 gained). Byte-identical HMM pooling is impossible without pinning `hmmsearch -Z/--domZ`, which would also change per-genome output. DECISION: keep HMM per-genome; `--pool_searches` is mmseqs-only. Fixed-`-Z` pooling is a future science decision. + - **UPDATE (branch `perf/hmm-pool-bitscore`): the bit-score-prefilter fix WORKS.** Replacing + `hmmsearch -E 1e-5` with a Z-independent bit-score prefilter `-T 10 --domT 10` (floor below + kofam's lowest profile threshold, 15.07) makes **pooled HMM byte-identical to per-genome**: on + extraves, per-genome(`-T`) vs pooled-chunked(`-T`) = IDENTICAL for all 5 DBs (141270 rows; kofam + fanned to 8 chunks). Costs are small: drift vs the `-E` baseline is **+63 kofam calls (~0.08%)** — + the hits `-E 1e-5` was silently dropping (bit-score is the intended KOfam filter, so arguably a + latent bug) — and the final annotations are the **same size** (34M, no bloat; only the transient + domtblout grows). So full-annotation pooling (mmseqs + HMM) is achievable + correct. Not in + PR #10 because `-T` changes HMM output vs `-E` → a deliberate pipeline-wide decision (re-run + finished cohorts for consistency), not mid-campaign. `-Z`/`--domZ` does NOT work: one pooled + `-Z` can't match the variable per-genome sizes. - **Single-task COMBINE — assessed, DEPRIORITISED.** `combine_annotations.py` is per-gene keyed by `[query_id, input_fasta]` (no cross-genome compute beyond the already-non-deterministic column order), so it parallelises cleanly — but only if **batched** (K genomes/batch → B From 0034dc9b29a320827a8d49fb7667fa706cd42713 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Taavi=20P=C3=A4ll?= Date: Sun, 28 Jun 2026 12:29:56 +0300 Subject: [PATCH 3/3] docs: reword O() task-count notation (reads as zeroes in monospace) Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/dev/search-bundling.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/dev/search-bundling.md b/docs/dev/search-bundling.md index 7bd88f72..3c045a22 100644 --- a/docs/dev/search-bundling.md +++ b/docs/dev/search-bundling.md @@ -16,7 +16,8 @@ only WITHIN a genome (megahit `k141_*` scaffolds collide across bins). Fixed by `PREFIX_GENES_FOR_POOL` (`___` on faa headers + gene-locs), stripped in `SPLIT_POOLED_HITS`. -The entire mmseqs annotation for a cohort is now O(#DBs) tasks. For merops+viral+methyl +The entire mmseqs annotation for a cohort is now #DBs tasks (one per database, not per +genome). For merops+viral+methyl at 56 bins: per-genome = 336 tasks (56 index + 168 search + 112 SQL) → pooled = **9** (1 index + 3 search + 3 split + 2 SQL). At 4621 bins: ~27,700 → 9. @@ -49,7 +50,7 @@ at 56 bins: per-genome = 336 tasks (56 index + 168 search + 112 SQL) → pooled - **Single-task COMBINE — assessed, DEPRIORITISED.** `combine_annotations.py` is per-gene keyed by `[query_id, input_fasta]` (no cross-genome compute beyond the already-non-deterministic column order), so it parallelises cleanly — but only if **batched** (K genomes/batch → B - batch-combines → one concat; per-genome would re-create the O(#bins) task explosion). It runs + batch-combines → one concat; per-genome would re-create the per-bin task explosion). It runs once, downstream, and is far cheaper than the (now-pooled) search fan-out, so the value is **robustness/memory** (de-risking the single big-memory COMBINE on ~4000-bin cohorts), not speed. Deferred: not worth the channel-batching + multi-dir staging + concat complexity vs the