Pangenome report enrichment: InterPro links, Pfam2GO, island loci, outlier flag - #100
Merged
Merged
Conversation
…1e-5 The retried SLURM run completed cleanly and the regenerated presence_matrix.tsv confirms the fix through the real pipeline: ada-1/ham-5 flip to absent in Mcir, recall 0.8/fp_rate 0.0 matching the earlier standalone script-level test. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
DIAMOND_SEARCH previously ran one Nextflow task (one SLURM job) per (query, target) genome pair. Real timing (cluster-vs-pairwise-sensitivity investigation) showed each pair's actual diamond compute is tiny (~1-3s even for this study's largest proteomes) -- job COUNT, not diamond runtime, is the likely cost driver on a shared SLURM queue with a 50/min submission rate limit, the same pattern nextflow.config's hmm_search_chunk_size comment already documents for FAMILY_HMMSEARCH. pairs_ch is now grouped by query genome (groupTuple), so one job searches every target for that query in a loop -- for pezizo_set1 (5 ingroup x 10 others) this cuts ~50 jobs to 5 for the gain direction, similarly for losses. Output files keep their exact original names/format, so PARSE_HITS needs no changes; the [meta_pair, file] tuple shape it expects is reconstructed downstream from each output file's self-describing name. Switched DIAMOND_SEARCH from storeDir to normal resume-based caching + publishDir: storeDir's skip-if-exists check isn't well-defined against a variable-length glob output, so partial-cache correctness with the batched job shape wasn't guaranteed. This also fixes a second problem the same investigation found: storeDir made this step's real cost invisible to every nextflow trace file, since a cache hit skips the task (and its trace row) entirely. Validated: a standalone test confirmed the groupTuple/flatten/filename- reconstruction channel logic produces correct per-pair tuples (2 batched jobs instead of 4 individual ones, correctly paired); a real diamond run through the new process produced output byte-identical to the existing published Spom_vs_Nirr result. Only modules/diamond.nf's shared DIAMOND_SEARCH process and its two callers (workflows/search.nf, workflows/loss_search.nf) are touched -- CONTEXT_DIAMOND_SEARCH and SINGLETON_DIAMOND_SEARCH are separate, differently-named processes, unaffected. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
… on Agaricales controls Added 5 real negative (housekeeping) controls to Agaricales.controls.csv (previously positive-only) and swept min_domain_evalue against agaricomycetes_mmseqs's existing family_hmmsearch domtblout data. Recall 1.0/fp_rate 0.0 at every threshold tested, including disabled -- a no-regression result rather than a replication of pezizo_set1's breakpoint, since these controls don't include a promiscuous-domain case. Closes the second-clade item in todo/validate-hmm-presence-coverage-broader-sweep.md for this specific parameter. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CHANGES.md documents the job-count reduction (~110 -> ~11 for an 11-genome study), the storeDir -> resume+publishDir caching switch and why, and the real end-to-end validation results (byte-identical output on a fresh, cache-cleared, containerized SLURM rerun of pezizo_set1). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Points readers at the published audit comparing the two presence models, right where the comparison table already lives. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Reusable DSL2 pangenome-profiling module for NI (ported from the one-off Afumigatus_pangenome study script chain in NovInvenio_ Investigations), plus a full round of fixes from an independent Opus review of the initial port. Status as of this commit: DONE and verified: - Core chain implemented: PREFIX_PROTEOME/GENOME -> CONCAT_PROTEOMES -> CLUSTER_TIER1 (mmseqs) -> PRESENCE_MATRIX -> per-strain rescue (EXTRACT_ABSENT_QUERIES -> TBLASTN_PER_STRAIN -> RESCUE_PASS -> EXTRACT_RESCUE_POSITIONS) -> FREQUENCY_BINS -> COOCCURRENCE -> PAIR_CLASSIFICATION, plus MASH-based DEREPLICATE/ASSIGN_CLADES and an optional captain-gene (Starship) classification step. - mash/scipy/zstd now provisioned in pixi.toml/Dockerfile/novinvenio.def (were missing entirely -- a real crash-on-first-run blocker). - Zero-ABSENT-cells edge case no longer crashes (--tblastn_tsv made optional in rescue_pass.py / extract_rescue_positions.py). - params.pangenome_project derived from --pangenome_project or the samplesheet basename, so two species runs can no longer collide on outdir/storeDir. - --help works; diamond-backend clustering hard-errors (unvalidated ID- matching path, no equivalent to restore_mmseqs_cluster_ids.py yet) instead of silently running; COOCCURRENCE's memory floor raised to match the real 32GB working value documented for the source study; MASH_SKETCH_ALL gated behind --pangenome_dereplicate; duplicate DNA basenames now raise instead of silently dropping a strain; missing GFF3s fail fast per-row instead of surfacing deep inside GENE_POSITIONS. - 463 tests passing (39 new for this module), 3 skipped, one pre-existing unrelated failure deselected (mmseqs AVX2 SIGILL on a non-AVX2 node -- environment issue, not this branch). - Real 10-strain smoke test (8 Afumigatus ingroup + 2 outgroup, real data) ran the full chain through rescue-pass DB setup as SLURM job 28428780 -- confirms mash+scipy+mmseqs all work end-to-end on real data on an AVX2 node. Outcome of the remaining steps (frequency bins -> co-occurrence -> pair classification) for that job not yet confirmed as of this commit. KNOWN OPEN (not done): - Full end-to-end confirmation of FREQUENCY_BINS/COOCCURRENCE/ PAIR_CLASSIFICATION on real data -- most recently added, least tested code path. - Real GFF3-parsing gap found testing against a second species (Coccidioides, funannotate-style GFF3s use CDS Parent= not protein_id=) -- bin/pangenome_build_gene_positions.py's parse_gff3_protein_positions() needs a fallback + a protein-FASTA cross-check before this module works on non-NCBI-style GFF3 input. Fix proposed, not yet applied as of this commit. - Diamond-backend clustering path unvalidated (hard-disabled, not fixed). - Mash sketch-once-per-genome optimization not implemented (currently sketches twice: full set + ingroup-only). - Tier-2/superfamily clustering not wired into the subworkflow. Not pushed to origin -- local checkpoint only.
Job 28428780 confirmed the full core chain (clustering through co-occurrence and pair classification) runs successfully on real data -- 57/57 processes, 0 failures. This was the last part of the module an independent review had flagged as unconfirmed; it is now confirmed. See CHANGES.md for the real output stats and what's still open (the Coccidioides GFF3-parsing fix in progress, diamond backend still hard-disabled). Not pushed to origin -- local checkpoint only.
… absent
Coccidioides' funannotate GFF3s (NovInvenio_Investigations'
coccidioides_pangenome study) carry no protein_id= attribute at all --
0/535 -- every CDS row uses Parent= instead. Previously this produced empty
position data silently (no error), cascading to empty family-position data
and collapsing every co-gain/co-loss candidate to insufficient_data
downstream.
- Prefer protein_id= when present (unchanged NCBI-style behavior); fall
back to CDS Parent= otherwise, split on comma for the rare
multiple-transcripts-share-one-CDS case (GFF3 spec allows Parent=A,B).
- Cross-check every resolved ID against that strain's actual protein FASTA
headers (new required --protein_dir arg) -- an attribute matching
syntactically doesn't prove it's the ID the FASTA actually uses.
- Hard-error if a strain's resolved-and-FASTA-matching fraction falls below
50% (genuine dialect mismatch); warn above 2% unresolved (stray
annotation-tool quirks, not structural).
- Wires protein_dir_abs ("<data_dir_abs>/pep") through
modules/pangenome/positions.nf and workflows/pangenome_profile.nf.
7 new tests in tests/test_pangenome_build_gene_positions.py (funannotate
dialect incl. Parent=A,B, mismatched-dialect hard-error path, warn-but-continue
path). Full suite: 471 passed, 3 skipped, no regressions (pixi run pytest).
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Deferred by user during the 2026-09-16 island/Pfam pipeline design discussion -- not scoped, just captured for later. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Generalizes Afumigatus study's synteny_windows.py/find_accessory_islands.py into a permanent pipeline script. Named marker tblouts (--marker_tblout NAME=PATH, repeatable) replace the study's hardcoded captain/sm_backbone flags, producing one has_<name> column per marker instead of a single column that would collapse independent signals. Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Ported from studies/fungi/Afumigatus_pangenome/bin/ summarize_island_functions.py's domain_enrichment(), split out of that study's combined island-annotation-plus-enrichment script into a standalone step (island annotation moved to pangenome_report_tables.py). Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…FAMILY_PFAM_SCAN/DOMAIN_ENRICHMENT modules MARKER_HMMSEARCH is new, tuple-based code (not an alias of the existing captain.nf modules) -- a repeated alias invocation is illegal Nextflow DSL2, caught by Opus review before implementation started. Not yet wired into pangenome_profile.nf (Task 8). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
per_strain_summary added per Opus report-completeness review: automates the kind of by-hand outlier discovery the Afumigatus study's own QC section relied on. Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
An earlier draft of this task silently dropped 3 of the 5 spec-required figures (frequency histogram, presence/absence raster, accumulation curve) and never added the spec-required top-enriched-domains figure -- caught by Opus review before implementation started. This restores all five, adds a domain-enrichment figure Afumigatus itself only ever made by hand (no checked-in generator), and adds a Heaps' law + core-genome-decay fit as a genuine statistical upgrade over that study's own report, which only asserted openness from the curve's visual shape. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Wires presence_matrix/frequency_table through for the restored figures (Task 5/6 corrections) alongside the original island/enrichment inputs. Not yet wired into pangenome_profile.nf (Task 8). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
NovInvenio_Investigations has its own bin/ni (resolve/fetch/discover/run, 7 commits, its own design doc) for operating on studies inside an already-scaffolded NII repo. This repo's bin/ni only ever did `init` (scaffold a brand-new NII-style repo, issue #76) -- a different tool for an earlier lifecycle stage, not a fork or duplicate. Same short command name for two disjoint CLIs was a real naming collision risk (whichever bin/ ends up first on $PATH silently shadows the other) and a documentation gap (NII's README calls its own bin/ni "the one entrypoint" without ever mentioning this one exists). Renamed the less-frequently-used tool (run once per new deploy repo, vs. NII's used continuously per study) rather than NII's, and added a cross-reference in this script's own docstring pointing to NII's bin/ni for the next lifecycle stage. Investigation and recommendation: NovInvenio_Investigations' notes/superpowers/research/ 2026-09-16-duplicate-ni-command-name.md. Not pushed, per this branch's established local-commit convention while the coccidioides-testing session has parallel uncommitted work in progress on the same branch.
Flag-gated behind the new --pangenome_island_pfam_hmm (distinct from the
pre-existing --pangenome_pfam_hmm, which stays load-bearing for the
captain-by-name branch only). Named marker searches
(--pangenome_marker_names/--pangenome_marker_hmm_paths) run via the new
MARKER_HMMSEARCH process, invoked once over a Channel.fromList fan-out --
not a loop, which DSL2 forbids. Verified via a real 10-strain execution
(not just --help) in both the zero-marker and one-marker code paths.
Two fixes surfaced by that real-execution smoke test, folded in here since
they were required to make the mandated test pass:
- bin/pangenome_{build_islands,domain_enrichment,select_background_reps,
report_tables,report_render}.py were missing their executable bit from
Tasks 4/7 (chmod +x only, no content change).
- FAMILY_PFAM_SCAN's hmmscan needs its target HMM database hmmpress-indexed;
Nextflow's path() input only stages the named file, not sibling .h3*
index files, so the process now presses the DB in its own work dir on
demand if the indices aren't already staged alongside it.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…re link Review findings from the whole-branch review, fixed before first production run: - F1: FAMILY_PFAM_SCAN's hmmscan used -E (sequence-level) instead of --domE (domain-level) for --pangenome_pfam_domain_evalue, and pangenome_domain_enrichment.py's parse_domtblout() ignored the param entirely (hardcoded max_ievalue=1e-3). Added --domain_evalue to both pangenome_domain_enrichment.py and pangenome_report_tables.py (which also calls parse_domtblout), threaded from params.pangenome_pfam_domain_evalue in modules/pangenome/pfam_enrichment.nf and modules/pangenome/report.nf. - F2: pangenome_report_tables.py wrote a blank/headerless islands_with_domains.tsv when there were zero islands (empty fieldnames list); added a fixed fallback header. Also added lineterminator="\n" to both DictWriter calls (islands_with_domains.tsv, per_strain_summary.tsv), which were emitting CRLF unlike every other TSV this script writes. - F3: pangenome_report_render.py's report.md unconditionally embedded figures/pair_classification_summary.png even though that plot is only drawn `if classification_counts_dict:`, producing a broken image link at small-N/zero-classification scale. Wrapped the embed line in the same conditional, matching the pattern used for the report's other guarded figures. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
select_background_reps and domain_enrichment's parse_domtblout could silently produce 0 records/0 hits from a nonempty input (e.g. a rep-FASTA-header vs family-ID convention drift), letting the pipeline report "no enriched domains" with exit 0 and no diagnostic. Both now warn loudly (matching pangenome_build_islands.accessory_islands's existing warnings.warn convention) when the join/parse fails despite a nonempty input, while staying silent for the legitimate empty-input cases. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…gh (F6) pangenome_build_islands.py's load_hit_families and pangenome_pair_classification.py's load_captain_families were near-verbatim duplicate hmmsearch/hmmscan --tblout parsers, each with a hardcoded id_sep="|" default and no CLI passthrough -- a study configured with a different --pangenome_id_sep silently got zero marker hits in both places (a pre-existing bug in pair_classification that build_islands copy-pasted). Consolidates the shared parsing loop into lib/pangenome_matrix.iter_tblout_family_hits (yields (short, family) pairs; each caller builds its own return shape -- a flat set for load_hit_families, a per-strain dict for load_captain_families). Adds --id_sep to both scripts' CLIs (default "|", unchanged behavior unless set) and threads params.pangenome_id_sep (already defined in nextflow.config) through modules/pangenome/islands.nf's BUILD_ISLANDS and modules/pangenome/pair_classification.nf's PAIR_CLASSIFICATION. Adds tests/test_pangenome_pair_classification.py (did not exist before) and new id_sep-threading coverage in tests/test_pangenome_matrix.py and tests/test_pangenome_build_islands.py. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…(F5) significant_islands.tsv / islands_with_domains.tsv carry a has_<marker> column per named marker search (e.g. has_captain), but neither REPORT_TABLES nor REPORT_RENDER ever summarized or mentioned them -- the design spec explicitly calls out marker-vs-island co-occurrence as a headline result, but the reporting half was never tasked. Adds pangenome_report_tables.py's marker_summary() + a new marker_summary.tsv output (marker_name, n_islands_with_marker, n_islands_total, pct_islands_with_marker; header-only when no marker search was run), wires it into REPORT_TABLES's output/emit and into REPORT_RENDER's input/CLI in modules/pangenome/report.nf and the workflow wiring in workflows/pangenome_profile.nf, and adds a "## Marker co-occurrence" section to pangenome_report_render.py's render_report_markdown() (skipped entirely when marker_summary.tsv has no data rows). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
fit_core_decay() only caught RuntimeError around scipy.optimize.curve_fit, but curve_fit raises TypeError when there are fewer data points than free parameters (e.g. a 2-strain cohort with a 3-parameter model) -- previously uncaught, would crash REPORT_RENDER, the pipeline's last step. Widened to except Exception, keeping the same fallback (core_inf=core_mean[-1], tau=nan, fit_ok=False). fit_heaps_law() computed np.log(pan_mean) and np.polyfit on it, which degenerates below ~3 strains (2 points fit a line trivially, 1 point can't fit at all). Added an early short-circuit below 3 strains returning a not-available sentinel (matching fit_core_decay's fit_ok convention: nan kappa/gamma/r_squared, is_open=False, fit_ok=False) instead of attempting the fit. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…ng, Pfam accession (F9) Three small, same-shape fixes flagged in review: 1. per_strain_summary()'s n_genes field actually counts families present (paralogs collapse to one family, not raw gene counts). Renamed to n_families throughout: the TSV column header, internal variable names in pangenome_report_tables.py, and pangenome_report_render.py's read of the column + its "Genes per strain" label (now "Families per strain"). 2. plot_domain_enrichment()'s docstring said results are sorted by fdr_q; the actual code (main()) sorts top_domains by fisher_p. Fixed the docstring to match the actual (correct) behavior. 3. The domain-enrichment output table (island_pfam_enrichment.tsv) reported only the Pfam model name (e.g. SnoaL_2), dropping the more stable Pfam accession (e.g. PF13577.9) needed for any future cross-reference (e.g. Pfam2GO, as the sibling Afumigatus study already does by accession). Added parse_domtblout_accessions() (domtblout column index 1 -- verified against the real hmmscan --domtblout column layout, not the column index named in the original review note) and a new pfam_accession column alongside the existing domain column. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
- Add bare_pfam_accession() helper to strip version suffix from Pfam accessions (e.g. PF13577.9 -> PF13577) - Add pfam_url column to island_pfam_enrichment.tsv output, positioned after pfam_accession - Generate InterPro links using bare accession: https://www.ebi.ac.uk/interpro/entry/pfam/{ACCESSION}/ - Update report rendering to display domain names as markdown links when pfam_url is available - Add comprehensive tests for bare_pfam_accession() and link rendering with fallback behavior Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Add bin/pangenome_pfam2go.py, a standalone CLI that maps Pfam accessions already found by hmmscan (island_pfam_enrichment.tsv) to GO terms via a standard pfam2go mapping file, using Task 1's bare_pfam_accession() to strip version suffixes before lookup. Extend pangenome_report_render.py's Pfam domain enrichment table to show a GO terms column when present on any row, omitted otherwise -- no new render function parameters yet (Task 5 wires the Nextflow process). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Adds add_island_locus() to bin/pangenome_report_tables.py, resolving each accessory island's member_families (rep-protein IDs) down to its example_strain's own protein coordinates (via lib/pangenome_matrix's read_cluster_tsv inverted rep->member, plus a --gene_positions TSV) into a citable strain:contig:start-end locus_id, with n_members_with_coordinates and n_contigs_in_locus as resolution/paralog-pull-in diagnostics. Wires it into main() with new --cluster_tsv/--gene_positions/--id_sep args, ahead of the islands_with_domains.tsv write, and extends its zero-island fallback header with the 6 new columns. render_report_markdown() gains an islands_with_domains_rows kwarg and a "Top islands (by size)" table in the Accessory islands section, shown only when at least one island resolved a real locus. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Two blocking bugs from the original task brief, both traced to prior review: 1. rep_to_members values are always Short<id_sep>protein_id-prefixed in real pipeline data, but gene_positions is keyed on the bare protein_id -- add_island_locus compared the prefixed string directly and never matched, silently sentinel-ing every island's locus. Fixed by splitting each candidate's prefix via bin/pangenome_build_presence_matrix.py's split_member_id() and restricting candidates to members whose Short prefix matches the island's own example_strain. 2. member_families in significant_islands.tsv is always comma-joined (bin/pangenome_build_islands.py:243), never id_sep-joined -- add_island_locus split on id_sep instead of "," (id_sep is still correct for the per-member prefix split in fix 1, just wrong here). Added a realistic-fixture regression test (test_add_island_locus_restricts_to_example_strain) using Short<id_sep>protein_id-prefixed member_to_rep keys, comma-joined member_families, and bare-keyed gene_positions -- the shape that let both bugs through the original test suite. Updated the other add_island_locus tests to the same realistic shapes. pixi run pytest tests/test_pangenome_report_tables.py tests/test_pangenome_report_render.py -v: 31 passed. pixi run pytest tests/ -q: 529 passed, 3 skipped (clean, matches prior 528+1 new test). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
Adds add_outlier_flags() to bin/pangenome_report_tables.py, applying the Iglewicz-Hoaglin modified z-score to each strain's singleton-family count to flag assembly-fragmentation artifacts (anomalously high singleton counts). Wired into per_strain_summary() and the per_strain_summary.tsv writer (singleton_z, is_outlier columns). render_report_markdown() gains a per_strain_rows parameter and lists flagged outlier strains in the Per-strain summary section. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…ty net Adds bin/verify_diamond_cluster_ids.py, the diamond-side counterpart to restore_mmseqs_cluster_ids.py: fails loud if any id in diamond's raw tier1_cluster.tsv doesn't match a real FASTA header verbatim, and gives a distinct, actionable error if an id only matches after stripping an mmseqs-style sp|/tr|/... wrapper (i.e. diamond behaved like mmseqs when it wasn't expected to). Verify-only, never rewrites the file -- unlike mmseqs, diamond isn't expected to need restoration, only confirmation. Wired into CLUSTER_TIER1's diamond branch (modules/pangenome/prefix_and_cluster.nf) immediately after diamond cluster runs, before the rep-seq reconstruction or anything downstream reads the file. This closes the gap that had --pangenome_cluster_backend diamond hard-blocked in pangenome.nf, so that guard is relaxed to a doc comment. Test coverage specifically includes the adversarial case flagged as untested in docs/adr/0003: a Short-prefixed, multi-pipe UniProt-style id (Afum|sp|ACC|NAME), which is the pangenome subworkflow's own header convention colliding with a raw UniProt defline. Not done: a real diamond-binary concordance benchmark against mmseqs on real pangenome data (no diamond/mmseqs binaries were available in this environment) -- tracked in todo/diamond-tier1-cluster-backend.md. Treat --pangenome_cluster_backend diamond as opt-in/experimental until that lands. bin/pangenome_cooccurrence.py and the report layer needed no changes -- confirmed backend-agnostic by inspection. Closes #98. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Ran a real diamond cluster invocation (pixi env, diamond v2.2.0.180) against an adversarial FASTA including the multi-pipe, Short-prefixed UniProt-style header (Afum|sp|O74225|YCF1_SCHPO) flagged in ADR-0003 as never having been tested against a real binary. The header came through *_cluster.tsv unchanged, and verify_diamond_cluster_ids.py passed against that real output -- confirms the header-preservation assumption for real, not just via unit tests against synthetic data. Still open: a real multi-strain biological dataset concordance benchmark (this smoke test used 5 small toy sequences, so it validates ID fidelity only, not clustering quality at scale) -- todo/diamond-tier1-cluster-backend.md updated accordingly. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…nt into the workflow Adds a new PFAM2GO module that annotates island Pfam enrichment results with GO terms, wires it into pangenome_profile.nf behind --pangenome_pfam2go (falling back to DOMAIN_ENRICHMENT's raw output when unset), and threads the new cluster_tsv/gene_positions inputs into REPORT_TABLES for island-locus columns. Verified by two real end-to-end smoke-test runs (with and without --pangenome_pfam2go) plus work-dir inspection confirming REPORT_TABLES consumes PFAM2GO's annotated file (not DOMAIN_ENRICHMENT's raw output) when the flag is set, and correctly bypasses it when the flag is unset. Full test suite: 535 passed, 3 skipped. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
…x Top-islands defensive reads (F3) REPORT_RENDER's main() was parsing islands_with_domains.tsv only to count n_islands, discarding the row data, and never passed islands_with_domains_rows to render_report_markdown() -- so the Task-3 "Top islands" table (and locus_id) could never appear in a real report.md despite being correctly implemented and tested at the function level. main() now keeps the parsed rows and passes them through; n_islands is derived from the same list instead of a separate discarding count. Also switches the Top-islands render block's row reads from direct dict indexing to .get(...) with fallbacks, matching this file's established convention for the other report sections. Adds a main()-level integration test (argv + tmp_path fixtures, matching tests/test_pangenome_report_tables.py's existing pattern) asserting a real locus_id reaches the written report.md. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
… coordinates (F2) add_island_locus() computed locus_start=min(starts)/locus_end=max(ends) across ALL resolved members even when they spanned more than one contig, fabricating a coordinate range that mixed positions from different contigs under a single contig label (e.g. S1:contigA:100-999999 where 999999 is actually on contigB) -- a scientifically wrong locus string. When n_contigs_in_locus > 1, locus_id/locus_contig/locus_start/locus_end now all fall back to the existing "-" sentinel (matching the zero-resolved-members case), while n_members_with_coordinates and n_contigs_in_locus keep their real computed counts so the multi-contig condition stays visible rather than silently collapsed. Updates test_add_island_locus_flags_multi_contig to assert the sentinel values, not just the count. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
'${params.pangenome_id_sep}' inside REPORT_TABLES's script block breaks
out of its single-quoting if the param ever contained a literal single
quote. Low severity (operator-controlled pipeline param, not untrusted
input) but cheap to make robust: escape embedded single quotes with the
standard bash close-quote/escaped-quote/reopen-quote technique.
Scope note: modules/pangenome/islands.nf has the identical
--id_sep '${params.pangenome_id_sep}' pattern, pre-existing and already
shipped from an earlier, unrelated plan/review round. Left untouched here
per the fix-round's scope (fix only the occurrence this plan added); a
human can decide separately whether to sweep islands.nf's copy too.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
… bug Found by a fix-round implementer while writing an unrelated test in the report-enrichment plan's final-review fix round; predates this plan and is out of its scope, filed here for separate pickup instead. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87
pangenome: unblock diamond tier-1 clustering with an ID-fidelity safety net
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds 4 independently-gated report-enrichment features to the pangenome accessory-island + Pfam-enrichment pipeline step (ported from the Afumigatus_pangenome study's newest analysis work):
--pangenome_pfam2go(operator-supplied mapping file; skipped when unset)strain:contig:start-endcoordinates per accessory island, folded into the existingREPORT_TABLESprocessProcess
Built via
superpowers:brainstorming→ spec (2 review rounds, Fable + this repo's own design conventions) →superpowers:writing-plans→superpowers:subagent-driven-development(6 tasks, each independently implemented and reviewed; Task 3 needed one fix round for two real bugs found in code review — a member-ID prefix mismatch and a wrong list-separator assumption, both traced to mistakes in the plan/brief rather than implementer error).A final whole-branch review (Opus) found one blocking issue (the island-locus table's data was parsed but never actually passed to the render function — the same class of bug as an earlier per-task fix) plus a real scientific-correctness bug (multi-contig islands were getting a fabricated coordinate mixing two different contigs' positions). Both fixed and independently re-reviewed clean.
Verified via 2 real end-to-end Nextflow executions (10-strain Coccidioides subset, real Pfam-A.hmm database) — one with
--pangenome_pfam2go, one without — with work-directory-level inspection proving the channel actually routes through the newPFAM2GOprocess when the flag is set, and correctly bypasses it when unset (not just unit-test coverage).536 tests passing, 3 skipped, 0 failed.
Test plan
--pangenome_pfam2go) on a 10-strain Coccidioides subset, 0 failed tasks in either runNote
Local
mainin this repo is currently 24 commits ahead oforigin/main(unrelated prior work, not part of this PR) — this PR's diff againstorigin/mainwill include that prior history in addition to this plan's 10 commits, per explicit instruction to open the PR as-is rather than syncmainfirst.🤖 Generated with Claude Code
https://claude.ai/code/session_01JocL68YDR4Rpfb5bqnrJ87