diff --git a/CHANGELOG.md b/CHANGELOG.md index 7645791..33b20be 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,8 +6,20 @@ All notable changes to Tessera are recorded here. The format follows ## [Unreleased] +Fixes from the post-1.2.0 audit +(`docs/superpowers/specs/2026-10-01-post-1.2.0-audit-design.md`, items A1-A10 and B1-B11). +For a given alignment the callers find the same regions, with the same coordinates and +p-values, as before. Three things around them do change: the `donor_undercovered` flag and +donor-absent rows where a coverage gap is now recognised as a breakpoint artefact; runs +that select the barcode caller on a panel where it cannot run (it no longer counts toward +`--min-methods`, and a barcode-only run is refused); and what reaches the scan -- an +alignment built with the MAFFT or MAF-based backends, or a panel built with `--curate`, +may now be a different (corrected) one. + ### Fixed +Data safety and alignment: + - **The MAFFT backend ignored strand.** A genome, or one contig of a draft assembly, on the opposite strand to the backbone was aligned as given and came out at chance-level identity (about 0.40 against 0.97 for the same genome in forward orientation), which the scan then @@ -63,6 +75,77 @@ All notable changes to Tessera are recorded here. The format follows reference and an empty MAFFT result are reported as input/output errors that name the file, instead of "Unexpected error". +Reporting: + +- **The run provenance named the wrong callers.** MaxChi, Bootscan, GENECONV and the barcode + caller were each described as `heuristic (min ... / margin ... / merge ...)` in + `run_provenance.json` and the report, so a default run recorded + `hmm + 3seq + heuristic + heuristic`. Each caller is now described under its own name, and + the record gains the settings that change what is reported: the agreement gate + (`--min-methods`), sibling exclusion, lineage clustering and donor re-attribution. +- **A clean recombinant between divergent parents was reported as a possible missing + reference.** A window straddling a breakpoint matches neither parent well on its own, so + its best similarity fell below the coverage threshold, the stretch was called a `divergent` + coverage gap, and the region was marked `donor_undercovered` -- on the shipped + `divergent` example the headline read "low confidence" for a donor identical to the query. + A gap within one window of a called region boundary is now labelled `breakpoint` in + `coverage_gaps.tsv` and the report when every under-threshold window in it is matched, + to the coverage threshold, by a single switch between the region's two parents. Such a + gap does not caveat the region, is not turned into a donor-absent region, and is left out + of the headline. A stretch from a source outside the panel stays `divergent` even when it + is short and sits beside a breakpoint. **`donor_undercovered` and the confidence wording + change for regions whose only gaps were breakpoint artefacts.** One consequence is not yet + measured on the hybrid harness: a region whose donor is a close stand-in for an absent + lineage (above the coverage threshold) was previously caveated only through its + breakpoint gaps, and is now reported without a caveat. + Reference recruitment (`fill-references`, `find-references`) is unaffected. +- The report's "covering N kb (P %) of the query" added the lengths of overlapping regions + that name different donors; it now reports their union. +- After `--reattribute-donors`, `recombination_methods.tsv` and the report's method table + kept the donor from before re-attribution. +- `--lineage-map` pointing at a file that does not exist was ignored (exit 0, untyped + report) by `recomb`, `type-lineages`, `detect`, `fill-references` and `build-panel`. It + is now an error. +- **The barcode caller on an untyped panel read as a negative.** It named the first record + of the alignment as major parent and its column in the method table said `no`. It now + reports no major parent; in an ensemble it is logged and shown as `not run`, the + agreement gate counts only callers that ran (a gate lowered for that reason is logged and + recorded as, for example, `1 (requested 2; 1 caller ran)`), and a run that selected only + `barcode` is refused. The report and provenance give the actual reason, which on a typed + panel is that fewer than two clades carry enough markers. +- The report judged the PHI p-value at alpha 0.05 whatever `--alpha` was. +- **A PHI test that could not reject was reported as "no signal".** With no more + informative sites than the window can hold (`--phi-window`, default 100) the permutation + p-value is 1 for any data. This is now reported as `not testable` (`NA` in the + `recombination_profile.tsv` header, `phi_p = None` in the API). +- Plots labelled a donor-absent region "recombinant: "; `similarity_pair` showed + the two leading window winners rather than the major parent and the leading donor; with + `--top-n 1` the donor was drawn grey. +- The report footer listed `.pdf` plots under `--plot-format png` and files that were not + written; the methods text and references described a two-caller ensemble; `--method` + help omitted `geneconv`. +- `sibeliaz` was probed with `-v`, which it rejects, and the error line was recorded as the + aligner version. A failed probe is no longer recorded as a version. +- Input checks: `find-references --msa ` and `recomb -o ` failed + with "Unexpected error"; `--max-rounds 0` exited 0 having built nothing; `reassort` + accepted `--ani-floor 500`, `--margin -3` and a `--dataset` key matching no segment; + `type-lineages` rejected `.fasta.gz` / `.fas` collections (it now reads what the other + commands read, and like them rejects a file in the collection that is not FASTA, where it + used to skip it); a multi-line aligner error broke `segment_scan.tsv`. + +### Changed + +- `RecombinationSignal.phi_p` is `float | None`. +- `seaborn` is no longer a dependency; nothing imported it. + +### Documentation + +- `docs/detection-methods.md` describes lineage clustering (including its limit on panels + below about 1.5 % divergence), the `breakpoint` coverage kind and the untestable-PHI case. +- `validation/README.md` states that the specificity harness defaults to `--min-methods 2` + while the CLI defaults to 1, gives the measured rate at each, and records that + `hcv_clonal_1b` currently fails. + ## [1.2.0] - 2026-10-01 ### Fixed diff --git a/docs/detection-methods.md b/docs/detection-methods.md index 27089c4..b0b9d15 100644 --- a/docs/detection-methods.md +++ b/docs/detection-methods.md @@ -78,6 +78,23 @@ Recombination regions (major parent: cowpox_KC813504): variola cowpox_KC813504 66268 147150 80882 0.999 0.977 0.97 2e-300 66768 ``` +### Near-duplicate references (lineage clustering) + +A recruited panel often holds several near-identical genomes of one lineage. Competed +individually they tie in every window and fragment the call, so the HMM caller first +pools them (`--cluster-lineages`, on by default; hmm only). Two references are pooled +when their identity stays at or above 98.5 % in every window, with no region-sized run +below it. The pooled lineage competes as one state under the label of its best-covering +member, and a region names that member. Clustering is skipped for panels of fewer than 4 +or more than 200 references. + +One limitation follows from the absolute threshold. On a panel in which *every* pair of +references is at least 98.5 % identical in every window -- mpox, VZV and other sets below +roughly 1.5 % divergence -- all references pool into a single lineage and the HMM has +nothing left to compete, so it calls no region whatever the data. The site-based callers +(3SEQ, MaxChi) are unaffected. On such a panel, run with `--no-cluster-lineages` to keep +the HMM's vote. `run_provenance.json` records whether clustering was on. + ### Low-divergence panels (intra-species sets, DNA viruses) When the references are nearly identical -- e.g. mpox clades (~0.5 %), VZV (~0.2 %), @@ -181,8 +198,10 @@ and scans the query window by window for the clade whose local markers it carrie where a non-major clade dominates is a region attributed to that lineage. A marker is denoised across the clade's members and the comparison is multi-way, so the call is robust to a single near-identical adjacent-clade genome winning by chance -- the failure mode of -the genome-level callers at low divergence. It is silent on untyped panels, and opt-in -(needs typed references). It composes with `--pool-consensus` but needs only a typed panel. +the genome-level callers at low divergence. It is opt-in and needs typed references. On an +untyped panel it cannot run, which is not the same as finding nothing: in an ensemble it is +logged and shown as `not run` in the method comparison, and a run that selected only +`barcode` is refused. It composes with `--pool-consensus` but needs only a typed panel. ## Parent-free recombination signal (PHI + Rmin) @@ -232,6 +251,17 @@ localizes the signal rather than averaging it away -- are the more informative parent-free outputs. The diagnostic runs for every `--method`; disable with `--no-phi`, or widen its window with `--phi-window`. +The PHI p-value is judged at the run's `--alpha`, in the report and for the per-region +flag alike. The test is **not testable** when the alignment has too few informative sites +for the window (`--phi-window`, default 100 site ranks; at least window + 2 sites are +needed): every pair of sites then falls inside one window, so reordering the sites cannot +change the statistic and the permutation p-value would be 1 whatever the data. Tessera +reports this as `not testable` (`NA` in the `recombination_profile.tsv` header) rather +than as a non-significant result, and no region is flagged `parent_free_support`; lower +`--phi-window` to test such an alignment. Just above that limit the test is defined but +has little power: the shipped `divergent` example (159 informative sites) gives p = 1 at +the default window and p = 0.001 at `--phi-window 20`. + This all remains an **indicative screen**: the built-in HMM and 3SEQ tests are fast triplet/segmentation screens, not a full tree-based analysis (such as GARD). Treat regions as candidates to confirm. @@ -276,16 +306,16 @@ The practical readings: |---|---| | `report.html` | Self-contained report: run provenance, the region table, the per-dataset stats, and an embedded interactive plot | | `recombination_regions.tsv` | Called regions: minor/major parent, start/end in **both MSA columns and query bases**, `length_bp` / `length_msa`, `support`, `pvalue` / `qvalue` with the `test` and `statistic` that produced them, mean similarities, the calling `methods`, and `parent_free_support` | -| `recombination_methods.tsv` | Ensemble breakdown (only when several methods run): one row per region with a Y/n per method and the parent-free flag | -| `recombination_profile.tsv` | Parent-free signal: header with the PHI p-value and Rmin, then per-informative-site local incompatibility (the PHI profile) | +| `recombination_methods.tsv` | Ensemble breakdown (only when several methods run): one row per region with `yes` / `no` per method (`not run` for a selected caller that could not run) and the parent-free flag | +| `recombination_profile.tsv` | Parent-free signal: header with the PHI p-value (`NA` when not testable) and Rmin, then per-informative-site local incompatibility (the PHI profile) | | `similarity_windows.tsv` | Full per-window matrix: `msa_position`, `query_position`, `winner`, and one similarity column per dataset. Always base-pair windows (identity over all comparable columns) | | `informative_site_windows.tsv` | Only under informative-site windowing: the windows the HMM segmented -- `msa_position`, `query_position`, the window's `msa_start` / `msa_end` (end-exclusive), `winner`, and per dataset the query's identity **at polymorphic columns only**. Not comparable with similarity or ANI | | `similarity_stats.tsv` | Per-dataset similarity statistics (median, windows above identity thresholds) | | `window_winners.tsv` | Per-dataset count of windows won (ties included) | -| `coverage_gaps.tsv` | Stretches where even the closest reference is a poor match -- possible missing references | +| `coverage_gaps.tsv` | Stretches where even the closest reference is below the best-similarity threshold, with a `kind`: `divergent` (the query is far from every reference -- a possible missing reference), `low_information` (too few comparable bases to judge), or `breakpoint` (windows straddling a called breakpoint: each under-threshold window is matched, to the threshold, by a single switch between the region's two parents; not a missing reference, and it does not caveat the region) | | `similarity_top{N}.{fmt}` | Static plot of the nearest `--top-n` datasets, called regions shaded | | `informative_sites_top{N}.{fmt}` | Only under informative-site windowing: the same plot for identity at informative sites | -| `similarity_pair.{fmt}` | Static plot of the major vs leading minor parent, region shaded | +| `similarity_pair.{fmt}` | Static plot of the major parent against the donor of the longest called region, regions shaded (the two leading window winners when no donor was called) | | `run_provenance.json` | Machine-readable record of the run: Tessera version, parameters, caller description, and -- when the alignment came from `tessera msa` -- the aligner, its version and its arguments | ### Reading `pvalue`, `support` and the lengths diff --git a/docs/reference-panels.md b/docs/reference-panels.md index 791cd26..03c3e5c 100644 --- a/docs/reference-panels.md +++ b/docs/reference-panels.md @@ -244,9 +244,10 @@ masking sibling (e.g. SARS-CoV-2 sublineages), add `--seed-keep-siblings`. NCBI By default (`--auto-diversify`), BLAST seeding will **switch to the `ncbi-virus` diversity path automatically** when it finds only siblings -- i.e. when the query's lineage saturates `nt` and no parental lineage can be recruited by similarity. A broad -fetch is capped (`--fetch-limit`, default 2000) and dereplicated; for a heavily -sequenced taxon the capped sample may miss lineages, so a curated `--candidate-pool` -is recommended (and the run says so). Whether the diversity panel actually contains the +fetch is not truncated: the whole set is downloaded and dereplicated locally, and the run +logs a notice when it exceeds `--fetch-limit` (default 2000) because that step can take +a few minutes. For a heavily sequenced taxon a curated `--candidate-pool` is the faster +route. Whether the diversity panel actually contains the parents depends on the taxon: it works when they are genotype/lineage representatives, less so for fine genotype-specific recombinants. Disable with `--no-auto-diversify`. This complements the caller-side defence: `tessera recomb` excludes whole-genome diff --git a/docs/superpowers/specs/2026-10-01-post-1.2.0-audit-design.md b/docs/superpowers/specs/2026-10-01-post-1.2.0-audit-design.md index 9195746..0f123fd 100644 --- a/docs/superpowers/specs/2026-10-01-post-1.2.0-audit-design.md +++ b/docs/superpowers/specs/2026-10-01-post-1.2.0-audit-design.md @@ -327,6 +327,22 @@ this section disagree with the text above, the plan and this section hold. 7/12, positive control 3/3, 55 bp); `run_validation.py` 6 PASS, 1 FAIL (`hcv_clonal_1b`), 1 SKIP; the divergent example reports high confidence. +**Plan B, after its whole-branch review** (fix commit `0c9daf2`) + +- *B2 tightened again:* "the two parents together explain the gap" was measured over the + whole gap, where the well-matching flanks dilute a short stretch from a source outside + the panel (150 bp at 15 % divergence beside a breakpoint was relabelled). Each + under-threshold window in the gap must now be matched, to the threshold, by the best + single-switch mosaic of the two parents. +- *B2, unmeasured:* a region whose donor is a close stand-in for an absent lineage (above + the coverage threshold) used to be caveated only through its breakpoint gaps and is now + uncaveated. `run_hybrids.py` scores `panel_donor_absent` cases on that caveat and was + not run; it is a gate before merge. +- *B6:* the reason a caller did not run is carried into the provenance and report; a + `--min-methods` gate lowered because a selected caller could not run is logged and + recorded with the requested value. +- The changelog's "no region call changes" now lists its exceptions. + **Plan C** - *C2:* the sign test stays on the HMM segment and only the reported span is trimmed. diff --git a/example_data/README.md b/example_data/README.md index 56c4326..b1cd5d0 100644 --- a/example_data/README.md +++ b/example_data/README.md @@ -16,9 +16,10 @@ tessera recomb --msa divergent.msa.fasta --query query --output out_divergent \ --window-size 300 --window-step 30 ``` -Both callers call `parent_B` over the insert (q-value ~1e-29) with a sharp breakpoint, +The four default callers all call `parent_B` over the insert with a sharp breakpoint, so the region is flagged as agreeing (high confidence); the similarity plot shows an -obvious crossover. +obvious crossover. The two short stretches listed under reference coverage are the +windows straddling the breakpoints (kind `breakpoint`), not missing references. ## `cryptic_insert.msa.fasta` -- why the ensemble exists @@ -31,12 +32,14 @@ tessera recomb --msa cryptic_insert.msa.fasta --query query --output out_hmm \ --window-size 1000 --window-step 100 --method hmm # 0 regions ``` -The default ensemble also runs 3SEQ, which pools the discriminating sites into an exact -triplet test and recovers the event (q-value ~1e-12, `methods` = 3seq): +The default ensemble also runs the site-based callers, which pool the discriminating +sites into triplet tests and recover the event (q-value ~1e-12, `methods` = 3seq,maxchi): ``` tessera recomb --msa cryptic_insert.msa.fasta --query query --output out_cryptic \ --window-size 1000 --window-step 100 # finds parent_B insert ``` -Both runs also report the parent-free PHI / Rmin signal in `recombination_profile.tsv`. +Both runs also write the parent-free PHI / Rmin signal to `recombination_profile.tsv`. +The cryptic example has 54 informative sites, fewer than the default PHI window of 100, so +its PHI test is reported as `not testable`; add `--phi-window 5` to test it. diff --git a/pyproject.toml b/pyproject.toml index 4765364..3725599 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,6 @@ dependencies = [ "numpy>=1.26,<3.0", "pandas>=2.0,<4.0", "matplotlib>=3.8,<4.0", - "seaborn>=0.13,<1.0", "plotly>=5.18,<7.0", "biopython>=1.83,<2.0", ] diff --git a/src/tessera/aligners/sibeliaz.py b/src/tessera/aligners/sibeliaz.py index 7884c4c..614e49d 100644 --- a/src/tessera/aligners/sibeliaz.py +++ b/src/tessera/aligners/sibeliaz.py @@ -47,7 +47,9 @@ class SibeliazAligner(Aligner): capabilities = ToolCapabilities( name="sibeliaz", conda=("bioconda::sibeliaz",), - required_binaries=(BinarySpec("sibeliaz", version_args=("-v",)),), + # The sibeliaz wrapper has no version option ("-v" is rejected as an illegal + # option), so there is nothing to probe; its version is recorded as unknown. + required_binaries=(BinarySpec("sibeliaz", version_args=None),), recommended_max_genomes=2000, threads_param="-t", ) diff --git a/src/tessera/cli/cmd_build_panel.py b/src/tessera/cli/cmd_build_panel.py index 3407344..d790788 100644 --- a/src/tessera/cli/cmd_build_panel.py +++ b/src/tessera/cli/cmd_build_panel.py @@ -13,7 +13,16 @@ import typer -from .main import _require_choice, _require_file, app, get_logger, stage_errors +from .main import ( + _require_choice, + _require_file, + _require_lineage_map, + _require_range, + _require_scan_windows, + app, + get_logger, + stage_errors, +) def _seed_source(candidate_pool, nextclade: bool, nextclade_dataset: str | None) -> str: @@ -84,6 +93,12 @@ def build_panel( logger = get_logger(output) with stage_errors(logger): _require_file(query, "Query file") + _require_lineage_map(lineage_map) + # Checked here, before any network or aligner work: a round count of zero + # builds nothing and still exits 0, and a bad window otherwise surfaces only + # after the panel has been recruited and aligned. + _require_range(max_rounds, "--max-rounds", lo=1) + _require_scan_windows(window_size, window_step) _require_choice(aligner, set(aligner_registry.names()), "--aligner") params = FillParams( query=query, collection=collection, output=output, diff --git a/src/tessera/cli/cmd_detect.py b/src/tessera/cli/cmd_detect.py index 97c5f7c..ad72850 100644 --- a/src/tessera/cli/cmd_detect.py +++ b/src/tessera/cli/cmd_detect.py @@ -11,7 +11,16 @@ import typer -from .main import _require_choice, _require_file, app, get_logger, stage_errors +from .main import ( + _require_choice, + _require_file, + _require_lineage_map, + _require_range, + _require_scan_windows, + app, + get_logger, + stage_errors, +) @app.command(name="detect") @@ -55,8 +64,8 @@ def detect( method: str = typer.Option( "hmm,3seq,maxchi,bootscan", "--method", help="Region caller(s): a comma-separated list of hmm/3seq/maxchi/bootscan/" - "heuristic, or 'all'. Several run as an ensemble and their regions are merged " - "(default hmm,3seq,maxchi,bootscan).", + "geneconv/barcode/heuristic, or 'all'. Several run as an ensemble and their " + "regions are merged (default hmm,3seq,maxchi,bootscan).", ), min_methods: int = typer.Option( 1, "--min-methods", @@ -113,6 +122,12 @@ def detect( logger = get_logger(output) with stage_errors(logger): _require_file(query, "Query file") + _require_lineage_map(lineage_map) + # Checked here, before any network or aligner work: a round count of zero + # builds nothing and still exits 0, and a bad window otherwise surfaces only + # after the panel has been recruited and aligned. + _require_range(max_rounds, "--max-rounds", lo=1) + _require_scan_windows(window_size, window_step) _require_choice(aligner, set(aligner_registry.names()), "--aligner") params = FillParams.for_detection( query=query, output=output, diff --git a/src/tessera/cli/cmd_fill_references.py b/src/tessera/cli/cmd_fill_references.py index 1849949..31be930 100644 --- a/src/tessera/cli/cmd_fill_references.py +++ b/src/tessera/cli/cmd_fill_references.py @@ -10,6 +10,9 @@ _require_choice, _require_directory, _require_file, + _require_lineage_map, + _require_range, + _require_scan_windows, app, get_logger, stage_errors, @@ -143,8 +146,8 @@ def fill_references( method: str = typer.Option( "hmm,3seq,maxchi,bootscan", "--method", help="Region caller(s) for the detection step: a comma-separated list of " - "hmm/3seq/maxchi/bootscan/heuristic, or 'all'. Several run as an ensemble " - "(default hmm,3seq,maxchi,bootscan).", + "hmm/3seq/maxchi/bootscan/geneconv/barcode/heuristic, or 'all'. Several run as " + "an ensemble (default hmm,3seq,maxchi,bootscan).", ), min_methods: int = typer.Option( 1, "--min-methods", @@ -203,6 +206,12 @@ def fill_references( logger = get_logger() with stage_errors(logger): _require_file(query, "Query file") + _require_lineage_map(lineage_map) + # Checked here, before any network or aligner work: a round count of zero + # builds nothing and still exits 0, and a bad window otherwise surfaces only + # after the panel has been recruited and aligned. + _require_range(max_rounds, "--max-rounds", lo=1) + _require_scan_windows(window_size, window_step) if collection is not None: _require_directory(collection, "Collection directory") _require_choice(aligner, set(aligner_registry.names()), "--aligner") diff --git a/src/tessera/cli/cmd_find_references.py b/src/tessera/cli/cmd_find_references.py index c38217a..e0e5807 100644 --- a/src/tessera/cli/cmd_find_references.py +++ b/src/tessera/cli/cmd_find_references.py @@ -6,7 +6,7 @@ import typer -from .main import app, get_logger, stage_errors +from .main import _require_file, _require_scan_windows, app, get_logger, stage_errors @app.command(name="find-references") @@ -87,6 +87,8 @@ def find_references( logger = get_logger() with stage_errors(logger): + _require_file(msa, "MSA file") + _require_scan_windows(window_size, window_step) params = FindRefParams( msa=msa, query=str(query), output=output, collection=collection, window_size=window_size, window_step=window_step, diff --git a/src/tessera/cli/cmd_reassort.py b/src/tessera/cli/cmd_reassort.py index a6dd35b..466bfc0 100644 --- a/src/tessera/cli/cmd_reassort.py +++ b/src/tessera/cli/cmd_reassort.py @@ -14,10 +14,11 @@ import typer from ..core.errors import UserInputError +from ..core.io import read_fasta from ..reassort import assign_segments from ..reassort.assign import DEFAULT_ANI_FLOOR from ..reassort.constellation import DEFAULT_MARGIN -from .main import _require_file, app, get_logger, stage_errors +from .main import _require_file, _require_range, app, get_logger, stage_errors @app.command(name="reassort") @@ -55,12 +56,26 @@ def reassort( logger = get_logger(output) with stage_errors(logger): _require_file(query, "Query file") + # ANI is a percentage. Above 100 nothing can be assigned; a negative margin + # leaves every near-best set empty, so a clonal pair reads as undetermined. + _require_range(ani_floor, "--ani-floor", lo=0.0, hi=100.0) + _require_range(margin, "--margin", lo=0.0) overrides: dict[str, str] = {} for item in dataset or []: if "=" not in item: raise UserInputError(f"--dataset must be SEGMENT=path, got {item!r}") seg, path = item.split("=", 1) overrides[seg.strip()] = path.strip() + if overrides: + # An override is looked up by segment name; one that matches no record would + # be ignored and that segment's dataset auto-detected instead. + segments = [name for name, _seq in read_fasta(query) if name] + unknown = sorted(set(overrides) - set(segments)) + if unknown: + raise UserInputError( + f"--dataset names segment(s) not in the query: {', '.join(unknown)}. " + f"Segments in {query.name}: {', '.join(segments) or '(none)'}." + ) result = assign_segments( query, dataset_overrides=overrides, @@ -102,7 +117,10 @@ def reassort( fo.write("segment\tintragenic_recombination\tn_regions\tnote\n") for sc in result.scans: flag = "yes" if sc.recombinant else ("no" if sc.scanned else "n/a") - fo.write(f"{sc.segment}\t{flag}\t{sc.n_regions}\t{sc.note}\n") + # A failure note quotes the aligner's own message, which can span + # lines and hold tabs; keep it to one cell. + note = " ".join(sc.note.split()) + fo.write(f"{sc.segment}\t{flag}\t{sc.n_regions}\t{note}\n") rollup = " | ".join(f"{sc.segment}: {sc.note}" for sc in result.scans) logger.info("Intragenic scan: %s", rollup) logger.info("Wrote %s", stsv) diff --git a/src/tessera/cli/cmd_recomb.py b/src/tessera/cli/cmd_recomb.py index 6e93c8b..f77f2c5 100644 --- a/src/tessera/cli/cmd_recomb.py +++ b/src/tessera/cli/cmd_recomb.py @@ -9,6 +9,8 @@ from .main import ( _require_choice, _require_file, + _require_lineage_map, + _require_output_directory, _require_range, app, get_logger, @@ -33,11 +35,13 @@ def recomb( "hmm,3seq,maxchi,bootscan", "--method", help="Region caller(s), comma-separated, or 'all'. Several run as an ensemble " "and their regions are merged into a consensus (agreement raises confidence); " - "the default is hmm,3seq,maxchi,bootscan (all but the legacy heuristic). Callers: " + "the default is hmm,3seq,maxchi,bootscan; geneconv, barcode and heuristic are " + "opt-in, and 'all' runs every one. Callers: " "hmm (HMM segmentation + a discordant-site " "significance test), 3seq (scan-aware triplet max-drawdown test; strong at low " "divergence), maxchi (chi-square triplet test, complementary to 3seq), bootscan " - "(distance + bootstrap support for the closest parent), barcode (clade-marker " + "(distance + bootstrap support for the closest parent), geneconv (longest " + "uninterrupted donor-match run), barcode (clade-marker " "lineage attribution; needs typed references), heuristic (legacy margin/merge). " "Pass a single name (e.g. --method hmm) for one caller.", ), @@ -148,6 +152,8 @@ def recomb( logger = get_logger() with stage_errors(logger): _require_file(msa, "MSA file") + _require_lineage_map(lineage_map) + _require_output_directory(output) _require_choice(plot_format, {"pdf", "png", "svg"}, "--plot-format") # Bound the numeric options here: out of range they reach the statistics and # surface as an internal exception under "Unexpected error". diff --git a/src/tessera/cli/cmd_type_lineages.py b/src/tessera/cli/cmd_type_lineages.py index 1821a0f..26fae7c 100644 --- a/src/tessera/cli/cmd_type_lineages.py +++ b/src/tessera/cli/cmd_type_lineages.py @@ -14,7 +14,8 @@ import typer from ..core.errors import UserInputError -from .main import _require_directory, app, get_logger, stage_errors +from ..core.io import _require_fasta, collection_genomes +from .main import _require_directory, _require_lineage_map, app, get_logger, stage_errors @app.command(name="type-lineages") @@ -51,12 +52,15 @@ def type_lineages( logger = get_logger(output) with stage_errors(logger): _require_directory(collection, "Collection directory") - genomes = sorted( - p for p in collection.iterdir() - if p.is_file() and p.suffix.lower() in (".fasta", ".fa", ".fna") - ) + _require_lineage_map(lineage_map) + # The same reading of a collection as every other command: each non-hidden + # file is a genome, whatever its extension and gzip-compressed or not, and a + # file that is not FASTA is an error rather than something to skip. + genomes = collection_genomes(collection) if not genomes: - raise UserInputError(f"No FASTA genomes found in {collection}") + raise UserInputError(f"No genome files found in {collection}") + for genome in genomes: + _require_fasta(genome) rows = assign_lineages( genomes, user_lineage_map=lineage_map, taxon=taxon, nextclade_dataset=nextclade_dataset, ref_ani_floor=ref_ani_floor, diff --git a/src/tessera/cli/main.py b/src/tessera/cli/main.py index 2239079..16f0374 100644 --- a/src/tessera/cli/main.py +++ b/src/tessera/cli/main.py @@ -78,6 +78,18 @@ def _require_file(path: Path, label: str) -> None: raise UserInputError(f"{label} is a directory, not a file: {path}") +def _require_lineage_map(path: Path | None) -> None: + """Reject a ``--lineage-map`` that does not exist. + + The readers treat a missing lineage file as "no typed names", which is right for + the automatically discovered ``lineages.tsv`` and wrong for a path the user typed: + the run would finish with an untyped report, and the barcode caller and donor + re-attribution would do nothing, without a word. + """ + if path is not None: + _require_file(path, "--lineage-map file") + + def _require_directory(path: Path, label: str) -> None: """Reject a missing or non-directory input directory. See :func:`_require_file`.""" if not Path(path).exists(): @@ -86,6 +98,22 @@ def _require_directory(path: Path, label: str) -> None: raise UserInputError(f"{label} is not a directory: {path}") +def _require_output_directory(path: Path) -> None: + """Reject an output path that exists and is not a directory. + + The writers create the directory if it is missing; given an existing file they fail + with ``[Errno 17] File exists`` from wherever the first output is written. + """ + if Path(path).exists() and not Path(path).is_dir(): + raise UserInputError(f"Output path exists and is not a directory: {path}") + + +def _require_scan_windows(window_size: int, window_step: int) -> None: + """The window options every scanning command shares; see :func:`_require_range`.""" + _require_range(window_size, "--window-size", lo=1) + _require_range(window_step, "--window-step", lo=1) + + def _parse_key_values(items: list[str], label: str) -> dict[str, str]: """Parse repeated ``key=value`` options into a dict (used for tool extras).""" out: dict[str, str] = {} diff --git a/src/tessera/core/binaries.py b/src/tessera/core/binaries.py index 445ff12..ac877c9 100644 --- a/src/tessera/core/binaries.py +++ b/src/tessera/core/binaries.py @@ -24,11 +24,13 @@ class BinarySpec: """A required external executable. ``version_args`` is the argument vector that prints a version (e.g. - ``("--version",)``). ``min_version`` is an optional dotted-string requirement. + ``("--version",)``), or ``None`` for a tool that has no version option: it is then + not executed at all and its version is recorded as ``unknown``. ``min_version`` is + an optional dotted-string requirement. """ name: str - version_args: tuple[str, ...] = field(default=("--version",)) + version_args: tuple[str, ...] | None = field(default=("--version",)) min_version: str | None = None @@ -54,6 +56,11 @@ def _query_version(name: str, version_args: tuple[str, ...]) -> str | None: parsed = _parse_version(blob) if parsed: return ".".join(map(str, parsed)) + # No dotted version in the output. A tool that exited cleanly said something about + # itself (a build date, say) and that is worth keeping; one that exited non-zero + # printed an error about the option, which is not a version. + if proc.returncode != 0: + return None return blob.strip().splitlines()[0] if blob.strip() else None @@ -70,7 +77,10 @@ def check_binaries(specs: tuple[BinarySpec, ...]) -> dict[str, str]: problems.append(f"{spec.name}: not found on PATH") continue - reported = _query_version(spec.name, spec.version_args) + reported = ( + None if spec.version_args is None + else _query_version(spec.name, spec.version_args) + ) versions[spec.name] = reported or "unknown" if spec.min_version is not None: diff --git a/src/tessera/recomb/barcode.py b/src/tessera/recomb/barcode.py index 4bd66ec..4adfd76 100644 --- a/src/tessera/recomb/barcode.py +++ b/src/tessera/recomb/barcode.py @@ -87,18 +87,21 @@ def call_regions_barcode(result: WindowSimilarity, analysis, window_size: int, p The query's backbone clade matches the most markers genome-wide; a run of windows where another clade matches its local markers better (by ``_MARGIN``) is a region for that donor clade. Returns ``(regions, major, [])`` to match ``call_regions``. + + When the caller cannot run -- the panel is untyped, or fewer than two clades carry + enough markers -- it returns ``([], None, [])``. ``major is None`` is the signal that + nothing was tested; the pipeline reports the caller as not run rather than as having + found nothing. """ from .clusters import _window_bounds from .regions import Region - labels = list(result.similarities) - default_major = labels[0] if labels else None lineage_map = getattr(params, "lineage_map", None) if not lineage_map: - return [], default_major, [] + return [], None, [] cols, alleles, rep = clade_markers(result.rows, result.query, lineage_map) if not cols: - return [], default_major, [] + return [], None, [] query = result.rows[result.query] clades = list(cols) diff --git a/src/tessera/recomb/coverage.py b/src/tessera/recomb/coverage.py index b0c4009..ee010a2 100644 --- a/src/tessera/recomb/coverage.py +++ b/src/tessera/recomb/coverage.py @@ -15,6 +15,12 @@ A gap is labelled ``divergent`` (ample comparable bases, the query really is far from all references -- a likely missing reference) or ``low_information`` (few comparable bases, so the gap is uncertain rather than informative). + +Once regions have been called, a ``divergent`` gap that sits on a called breakpoint +and is explained by that region's two parents together is relabelled ``breakpoint`` +(:func:`mark_breakpoint_gaps`): a window straddling a switch between divergent +parents matches neither of them well on its own, which is a property of the window, +not a missing reference. """ from __future__ import annotations @@ -27,7 +33,10 @@ import numpy as np from .regions import Region -from .similarity import WindowSimilarity +from .similarity import WindowSimilarity, _canonical_mask + +# The kind given to a gap that a called breakpoint explains. See mark_breakpoint_gaps. +BREAKPOINT_KIND = "breakpoint" @dataclass @@ -68,7 +77,7 @@ class CoverageGap: n_windows: int best_label: str # the closest reference across the gap (still a poor match) mean_best: float # its mean similarity over the gap - kind: str # "divergent" | "low_information" + kind: str # "divergent" | "low_information" | "breakpoint" def coverage_threshold(result: WindowSimilarity, params: CoverageParams) -> float: @@ -140,6 +149,103 @@ def span(index: int) -> tuple[int, int]: return gaps, threshold +def _mosaic_identity( + rows: dict[str, np.ndarray], query: str, parents: tuple[str, str], start: int, end: int +) -> float: + """Best identity of the query, over columns ``[start, end)``, to a single-switch + mosaic of the two ``parents`` (either order, any switch point). + + This is what a window straddling one breakpoint between the two parents can reach: + one parent up to the switch, the other after it. Identity is counted as everywhere + else -- matches over the columns where the query and the parent in force both carry + a canonical base -- so the value is on the same scale as the per-window similarity + the coverage threshold is applied to. ``nan`` when nothing is comparable. + """ + q = rows[query][start:end] + q_canon = _canonical_mask(q) + comp, match = [], [] + for label in parents: + ref = rows[label][start:end] + canon = q_canon & _canonical_mask(ref) + comp.append(np.concatenate(([0], np.cumsum(canon)))) + match.append(np.concatenate(([0], np.cumsum(canon & (q == ref))))) + best = float("nan") + for first, second in ((0, 1), (1, 0)): + # Switch after k columns: `first` explains [0, k), `second` explains [k, n). + num = match[first] + (match[second][-1] - match[second]) + den = comp[first] + (comp[second][-1] - comp[second]) + valid = den > 0 + if valid.any(): + value = float(np.max(num[valid] / den[valid])) + best = value if isnan(best) else max(best, value) + return best + + +def mark_breakpoint_gaps( + gaps: list[CoverageGap], + regions: list[Region], + result: WindowSimilarity, + window_size: int, + threshold: float, +) -> int: + """Relabel ``divergent`` gaps that a called breakpoint explains; return how many. + + A window straddling a breakpoint between two divergent parents is part one parent + and part the other, so its similarity to either alone falls below the coverage + threshold although both are in the panel. A gap is relabelled ``breakpoint`` when + + - it lies within one window width of a boundary of a called, donor-present region + (the only place a straddling window can be), and + - every under-threshold window inside it is explained by the region's two parents: + the query's identity to the best single-switch mosaic of donor and major parent + over that window reaches ``threshold``. + + The second condition asks, window by window, the question the coverage test asks, + with a mosaic of the two parents in place of a single reference. It keeps a stretch + from a source outside the panel labelled ``divergent`` even when it is short and + sits right beside a breakpoint: averaged over the whole gap such a stretch is + diluted by the well-matching flanks, but the windows that hold it stay below the + threshold whatever the switch point. ``result`` is the base-pair scan the gaps were + called on. The gaps are mutated in place; call this after region calling and before + :func:`gaps_as_regions`, which skips every kind other than ``divergent``. + """ + half = window_size // 2 + relabelled = 0 + for gap in gaps: + if gap.kind != "divergent": + continue + windows = [ + (pos - half, pos - half + window_size) + for pos, best in zip(result.positions, result.best_sim, strict=True) + if gap.msa_start <= pos - half and pos - half + window_size <= gap.msa_end + and not isnan(best) and best < threshold + ] + if not windows: + continue + for region in regions: + if region.donor_absent: + continue + parents = (region.minor_parent, region.major_parent) + if any(label not in result.rows for label in parents): + continue + near = any( + boundary - window_size <= gap.msa_start + and gap.msa_end <= boundary + window_size + for boundary in (region.msa_start, region.msa_end) + ) + if not near: + continue + explained = [ + _mosaic_identity(result.rows, result.query, parents, start, end) + for start, end in windows + ] + if all(not isnan(value) and value >= threshold for value in explained): + gap.kind = BREAKPOINT_KIND + relabelled += 1 + break + return relabelled + + def flag_undercovered_regions(regions: list[Region], threshold: float) -> None: """Mark called regions whose donor is itself a poor match (possible stand-in).""" for r in regions: diff --git a/src/tessera/recomb/diagnostics.py b/src/tessera/recomb/diagnostics.py index b88d663..6c51d34 100644 --- a/src/tessera/recomb/diagnostics.py +++ b/src/tessera/recomb/diagnostics.py @@ -54,7 +54,9 @@ class RecombinationSignal: """Parent-free recombination evidence for an alignment.""" n_informative: int # biallelic informative columns used - phi_p: float # PHI permutation p-value (one-sided; small = recombination) + # PHI permutation p-value (one-sided; small = recombination). None when the test + # could not have rejected: see recombination_signal. + phi_p: float | None phi_observed: float # the windowed mean incompatibility phi_window: int # window width, in informative-column ranks rmin: int # Hudson-Kaplan minimum number of recombination events @@ -231,6 +233,12 @@ def recombination_signal( as-is; the significance threshold is a reporting concern, applied downstream. ``query_label`` is accepted for interface symmetry; the statistics use every sequence in ``rows``. + + ``phi_p`` is ``None`` when the test is not testable at this window: with ``z`` + informative columns and ``window >= z - 1`` every pair of columns falls inside the + window, the statistic is the mean over all pairs, and no reordering of the columns + can change it -- the permutation p-value would be 1 whatever the data. Rmin, the + intervals and the profile do not depend on the permutation and are still returned. """ allele1, allele0, positions = biallelic_columns(rows) z = positions.size @@ -243,7 +251,11 @@ def recombination_signal( allele1, allele0, positions = allele1[keep], allele0[keep], positions[keep] z = positions.size incompatible = incompatibility_matrix(allele1, allele0) - p, observed = phi_pvalue(incompatible, window, seed=seed) + p: float | None + if z - 1 > window: + p, observed = phi_pvalue(incompatible, window, seed=seed) + else: + p, observed = None, phi(incompatible, window) rmin, intervals = hudson_kaplan_rmin(incompatible, positions) query_intervals = [(int(column_to_query(a)), int(column_to_query(b))) for a, b in intervals] profile = phi_profile(incompatible, positions, column_to_query, window) diff --git a/src/tessera/recomb/report.py b/src/tessera/recomb/report.py index ca9af59..fa8d924 100644 --- a/src/tessera/recomb/report.py +++ b/src/tessera/recomb/report.py @@ -16,7 +16,9 @@ - ``window_winners.tsv`` per-dataset window-win counts (ties included) - ``recombination_regions.tsv`` called recombinant regions - ``similarity_top{N}.{fmt}`` static top-N similarity plot (regions shaded) -- ``similarity_pair.{fmt}`` static major-vs-minor pairwise plot +- ``similarity_pair.{fmt}`` static pairwise plot: the major parent against the + donor of the longest called region (the two leading + window winners when no donor was called) - ``report.html`` self-contained summary (tables + interactive plot) Under informative-site windowing (near-identical panels) the windows the HMM @@ -56,10 +58,25 @@ "ReportContext", "print_summary", "print_regions", "print_coverage", "build_interactive_figure", "plot_top_n", "plot_pairwise", - "write_html_report", "write_reports", + "write_html_report", "write_reports", "pair_datasets", ] +def pair_datasets(regions: list[Region], ranked: list[str]) -> list[str]: + """The two datasets of the pairwise plot: major parent, then the leading donor. + + The leading donor is the donor of the longest called, donor-present region. Without + one there is no minor parent to show and the plot falls back to ``ranked`` (the two + leading window winners) -- which is not "major versus minor" when the panel holds a + near-duplicate of the backbone, hence the region-based choice whenever possible. + """ + present = [r for r in regions if not r.donor_absent] + if not present: + return ranked + longest = max(present, key=lambda r: r.length_bp) + return [longest.major_parent, longest.minor_parent] + + def write_reports( result: WindowSimilarity, analysis: AnalysisResult, @@ -76,36 +93,54 @@ def write_reports( output_dir.mkdir(parents=True, exist_ok=True) gaps = ctx.gaps - write_windows_tsv(result, per_window_winners, output_dir, logger) - write_stats_tsv(analysis, output_dir, logger) - write_winners_tsv(analysis, output_dir, logger) + # The files written beside the report, in the order the footer lists them. Built as + # they are written, so the footer never names a file that does not exist or a plot + # in a format that was not asked for. + companions: list[str] = [] + write_regions_tsv(regions, output_dir, logger) + companions.append("recombination_regions.tsv") + if len(ctx.methods_run) > 1 and ctx.method_breakdown is not None: + write_methods_tsv( + ctx.method_breakdown, ctx.methods_run, output_dir, logger, ctx.methods_not_run + ) + companions.append("recombination_methods.tsv") write_coverage_tsv(gaps, ctx.coverage_threshold, output_dir, logger) + companions.append("coverage_gaps.tsv") if ctx.signal is not None: write_profile_tsv(ctx.signal, output_dir, logger) - if len(ctx.methods_run) > 1 and ctx.method_breakdown is not None: - write_methods_tsv(ctx.method_breakdown, ctx.methods_run, output_dir, logger) + companions.append("recombination_profile.tsv") + write_winners_tsv(analysis, output_dir, logger) + write_stats_tsv(analysis, output_dir, logger) + write_windows_tsv(result, per_window_winners, output_dir, logger) + companions += ["window_winners.tsv", "similarity_stats.tsv", "similarity_windows.tsv"] # Rank on the windows the caller segmented: on a near-identical panel base-pair # windows tie almost everywhere, so they order the references poorly. ranking = ctx.site_analysis or analysis top_datasets = rank_datasets(ranking, top_n) logger.info("Top %d nearest datasets: %s", len(top_datasets), ", ".join(top_datasets)) - plot_top_n(result, top_datasets, regions, output_dir, plot_format, logger) + plots = [plot_top_n(result, top_datasets, regions, output_dir, plot_format, logger)] if ctx.site_result is not None: write_site_windows_tsv( ctx.site_result, winners_per_window(ctx.site_result), output_dir, logger ) - plot_top_n( + companions.append("informative_site_windows.tsv") + plots.append(plot_top_n( ctx.site_result, top_datasets, regions, output_dir, plot_format, logger, ylabel="Identity to query at informative sites", title=f"Identity at informative sites, query {result.query}", stem="informative_sites_top", - ) + )) - pair = rank_datasets(ranking, 2) - plot_pairwise(result, pair, regions, output_dir, plot_format, logger) + pair = pair_datasets(regions, rank_datasets(ranking, 2)) + plots.append(plot_pairwise(result, pair, regions, output_dir, plot_format, logger)) + # A plot function returns None when it had nothing to draw and wrote no file. + companions += [path.name for path in plots if path is not None] + if (output_dir / "run_provenance.json").exists(): + companions.append("run_provenance.json") write_html_report( - result, analysis, regions, top_datasets, provenance, output_dir, logger, ctx + result, analysis, regions, top_datasets, provenance, output_dir, logger, ctx, + companion_files=companions, ) diff --git a/src/tessera/recomb/report_assets.py b/src/tessera/recomb/report_assets.py index aa0d21c..e51bead 100644 --- a/src/tessera/recomb/report_assets.py +++ b/src/tessera/recomb/report_assets.py @@ -108,12 +108,16 @@ "The scan slides a fixed-width window along the alignment in fixed steps; each " "window is scored independently."), ("Support", - "The share of distinguishing (discordant) sites -- where the query matches one " - "candidate parent but not the other -- that favour the donor. 0.5 = no " - "preference, 1.0 = every distinguishing site favours the donor."), + "The supporting statistic of the test behind the region; its meaning depends on " + "the caller, and the statistic column of recombination_regions.tsv names it. For " + "the HMM it is the share of distinguishing (discordant) sites -- where the query " + "matches one candidate parent but not the other -- that favour the donor " + "(0.5 = no preference, 1.0 = every distinguishing site favours the donor)."), ("q-value", - "The sign-test p-value after Benjamini-Hochberg correction across all candidate " - "segments (false-discovery-rate control). A region is reported when q <= alpha."), + "The p-value of the test behind the region after Benjamini-Hochberg correction " + "across that caller's own candidates (false-discovery-rate control within one " + "caller). A region is reported when q <= alpha. For a region several callers " + "found, the most significant caller's value is shown."), ("Breakpoint", "The query position where the source switches, with a posterior-derived " "uncertainty interval from the HMM."), @@ -137,6 +141,13 @@ ("3SEQ triplet test", "Boni MF, Posada D, Feldman MW (2007). An exact nonparametric method for inferring " "mosaic structure in sequence triplets. Genetics 176(2):1035-1047."), + ("MaxChi", + "Maynard Smith J (1992). Analyzing the mosaic structure of genes. Journal of " + "Molecular Evolution 34(2):126-129."), + ("Bootscan", + "Salminen MO, Carr JK, Burke DS, McCutchan FE (1995). Identification of breakpoints " + "in intergenotypic recombinants of HIV type 1 by bootscanning. AIDS Research and " + "Human Retroviruses 11(11):1423-1425."), ("PHI test", "Bruen TC, Philippe H, Bryant D (2006). A simple and robust statistical test for " "detecting the presence of recombination. Genetics 172(4):2665-2681."), diff --git a/src/tessera/recomb/report_context.py b/src/tessera/recomb/report_context.py index 5b1bd3b..77f1e7e 100644 --- a/src/tessera/recomb/report_context.py +++ b/src/tessera/recomb/report_context.py @@ -13,7 +13,7 @@ from dataclasses import dataclass from .analyze import AnalysisResult -from .coverage import CoverageGap +from .coverage import BREAKPOINT_KIND, CoverageGap from .diagnostics import RecombinationSignal from .similarity import WindowSimilarity from .typing import LineageMap @@ -29,8 +29,16 @@ class ReportContext: lineage_map: LineageMap | None = None query_lineage: str | None = None signal: RecombinationSignal | None = None + # The run's significance level, so the report judges the PHI p-value at the same + # alpha the per-region corroboration used. + alpha: float = 0.05 organism: str | None = None methods_run: tuple[str, ...] = () + # Selected callers that could not run (barcode on an untyped panel). They stay in + # ``methods_run`` so the method table keeps a column for them, marked "not run". + methods_not_run: tuple[str, ...] = () + # Why they could not run, in the run's own words (the same text the log gives). + not_run_reason: str = "" method_breakdown: list[dict] | None = None per_major: dict[str, str] | None = None # Set only when the scan used informative-site windowing: the per-window identity @@ -44,3 +52,14 @@ class ReportContext: def gaps(self) -> list[CoverageGap]: """``coverage_gaps`` with ``None`` normalised to an empty list.""" return self.coverage_gaps or [] + + @property + def caveat_gaps(self) -> list[CoverageGap]: + """The gaps that may mean a missing reference: every kind except ``breakpoint``. + + A breakpoint gap is a window straddling a called breakpoint. It is listed in the + coverage table and ``coverage_gaps.tsv`` under its own kind, but it is not a + poorly covered stretch, so the headline caveat, the mosaic and the plots leave + it out. + """ + return [g for g in self.gaps if g.kind != BREAKPOINT_KIND] diff --git a/src/tessera/recomb/report_html.py b/src/tessera/recomb/report_html.py index 7eac156..4513cf0 100644 --- a/src/tessera/recomb/report_html.py +++ b/src/tessera/recomb/report_html.py @@ -18,7 +18,7 @@ from .regions import Region from .report_assets import _CSS, _GLOSSARY, _REFERENCES from .report_context import ReportContext -from .report_plots import GREY, _color_map, build_interactive_figure +from .report_plots import GREY, _color_map, build_interactive_figure, region_labels from .similarity import WindowSimilarity from .typing import LineageMap, typed @@ -36,6 +36,22 @@ def _swatch(color: str) -> str: return f'' +def _union_length(spans: list[tuple[int, int]]) -> int: + """Total length covered by ``spans`` (half-open intervals), overlaps counted once.""" + total = 0 + covered_to: int | None = None + for start, end in sorted(spans): + if end <= start: + continue + if covered_to is None or start > covered_to: + total += end - start + covered_to = end + elif end > covered_to: + total += end - covered_to + covered_to = end + return total + + def _summary( result: WindowSimilarity, regions: list[Region], datasets: list[str] ) -> dict: @@ -50,7 +66,9 @@ def _summary( major = "n/a" present = [r for r in regions if not r.donor_absent] absent = [r for r in regions if r.donor_absent] - recomb_bp = sum(max(0, r.query_end - r.query_start) for r in present) + # The union, not the sum: overlapping regions that name different donors are kept + # as separate rows, and adding their lengths counts the shared stretch twice. + recomb_bp = _union_length([(r.query_start, r.query_end) for r in present]) minors: list[str] = [] for r in present: if r.minor_parent not in minors: @@ -260,11 +278,20 @@ def _regions_html( def _method_comparison_html( breakdown: list[dict], methods_run: tuple[str, ...], per_major: dict[str, str], - lineage_map: LineageMap | None = None, + lineage_map: LineageMap | None = None, methods_not_run: tuple[str, ...] = (), + not_run_reason: str = "", ) -> str: """A compact region x method agreement matrix; omitted for a single-method run.""" if len(methods_run) < 2: return "" + not_run_note = "" + if methods_not_run: + names = ", ".join(html.escape(m) for m in methods_not_run) + not_run_note = ( + f'

Not run: {names}' + f'{" (" + html.escape(not_run_reason) + ")" if not_run_reason else ""}. ' + f'It tested nothing; its column below is not a negative result.

' + ) majors = ", ".join( f'{html.escape(m)} → ' f'{html.escape(typed(per_major.get(m, "n/a"), lineage_map))}' @@ -273,7 +300,7 @@ def _method_comparison_html( intro = ( f'

Each caller ran independently on the same alignment; a region ' f'found by more than one is more trustworthy (and raises the confidence above). ' - f'Backbone per method — {majors}.

' + f'Backbone per method — {majors}.

{not_run_note}' ) if not breakdown: return intro + '

No regions were called by any method.

' @@ -285,6 +312,7 @@ def _method_comparison_html( rows = "" for b in breakdown: cells = "".join( + 'not run' if m in methods_not_run else f'{"✓" if m in b["per_method_support"] else "·"}' for m in methods_run ) @@ -302,10 +330,12 @@ def _method_comparison_html( def _method_section( method_breakdown: list[dict] | None, methods_run: tuple[str, ...], per_major: dict[str, str] | None, lineage_map: LineageMap | None, + methods_not_run: tuple[str, ...] = (), not_run_reason: str = "", ) -> str: """Wrap the method-comparison table in a report section (empty for one method).""" body = _method_comparison_html( - method_breakdown or [], methods_run, per_major or {}, lineage_map + method_breakdown or [], methods_run, per_major or {}, lineage_map, methods_not_run, + not_run_reason, ) if not body: return "" @@ -361,34 +391,39 @@ def _methods_html(provenance: dict[str, str]) -> str: ) return ( '
Methods & glossary' - '

Tessera segments the query against the reference panel with an HMM ' - '(jpHMM-style) and reports a region only when its donor beats the major parent on the ' - 'sites that distinguish them (a sign test on discordant sites, immune to window ' - 'overlap; Benjamini-Hochberg FDR across segments), with a posterior breakpoint ' - 'interval. By default it runs an ensemble of callers (HMM and the 3SEQ triplet ' - 'test) and merges their regions into one consensus, so a region found by more than ' - 'one method is flagged as agreeing and treated as higher confidence (see the ' - 'caller line under Run parameters). It remains an indicative screen, not a full ' - 'phylogenetic test (e.g. GARD) -- confirm strong candidates.

' + '

Tessera runs one or more region callers on the same alignment and ' + 'merges their regions into one consensus; a region found by more than one caller ' + 'is flagged as agreeing and treated as higher confidence. The callers that ran ' + 'for this report are named in the caller line under Run parameters. The default ' + 'ensemble is four callers: an HMM segmentation of the query against the reference ' + 'panel (jpHMM-style), which reports a segment only when its donor beats the major ' + 'parent on the sites that distinguish them (a one-sided sign test on discordant ' + 'sites, with a posterior breakpoint interval); the 3SEQ triplet test; MaxChi; and ' + 'Bootscan. Each caller applies Benjamini-Hochberg correction within its own ' + 'candidates; nothing is corrected across callers or across the genome. This is an ' + 'indicative screen, not a full phylogenetic test (e.g. GARD) -- confirm strong ' + 'candidates.

' f'
{glossary}
' f'

References

    {references}
' f'

Run parameters

{params}
' ) -def _footer_html(provenance: dict[str, str]) -> str: - files = [ - "recombination_regions.tsv", "recombination_methods.tsv", "coverage_gaps.tsv", - "recombination_profile.tsv", "window_winners.tsv", "similarity_stats.tsv", - "similarity_windows.tsv", "similarity_top*.pdf", "similarity_pair.pdf", - ] - flist = ", ".join(f"{f}" for f in files) +def _footer_html(provenance: dict[str, str], companion_files: list[str] | None = None) -> str: + """The version line and the companion files that were written with this report. + + ``companion_files`` comes from the writer that produced them. With none given the + sentence is left out: a list is only worth printing if it is the list of this run. + """ ver = html.escape(provenance.get("tessera version", "")) date = html.escape(provenance.get("date (UTC)", "")) + companions = "" + if companion_files: + flist = ", ".join(f"{html.escape(f)}" for f in companion_files) + companions = f"
Companion files in this folder: {flist}.
" return ( f'
Generated by Tessera {ver} · ' - f'{date} UTC
' - f'
Companion files in this folder: {flist}.
' + f'{date} UTC{companions}' ) @@ -419,7 +454,9 @@ def _coverage_html(gaps: list[CoverageGap], threshold: float) -> str: f'{threshold:.3f} best-similarity threshold. ' f'divergent = the query is genuinely far from every reference ' f'(a likely missing reference); low information = too few comparable ' - f'bases to judge.

' + f'bases to judge; breakpoint = windows straddling a called ' + f'breakpoint, where the two parents together explain the query (not a missing ' + f'reference, and not counted in the caveat above).

' ) if not gaps: return ( @@ -468,11 +505,22 @@ def _signal_html(signal: RecombinationSignal | None, alpha: float = 0.05) -> str '

Too few informative sites in the alignment for a parent-free ' 'recombination test.

' ) - significant = signal.phi_p < alpha - verdict = ( - 'significant recombination signal' if significant - else 'no significant recombination signal' - ) + if signal.phi_p is None: + # Every pair of informative sites is inside one window, so the permutation test + # cannot reject whatever the data. Say so instead of printing p = 1. + p_cell = "not testable" + verdict = ( + f'every pair of the {signal.n_informative} informative sites lies within the ' + f'window of {signal.phi_window} ranks, so the permutation test cannot reject here; ' + f'this is not evidence against recombination. Lower ' + f'--phi-window to test' + ) + else: + p_cell = f"p = {signal.phi_p:.4g}" + verdict = ( + 'significant recombination signal' if signal.phi_p < alpha + else 'no significant recombination signal' + ) + f' (alpha {alpha:g}; {signal.n_informative} informative sites)' intervals = ", ".join( f"{_fmt_int(lo)}–{_fmt_int(hi)}" for lo, hi in signal.rmin_intervals[:8] ) @@ -487,9 +535,8 @@ def _signal_html(signal: RecombinationSignal | None, alpha: float = 0.05) -> str 'forces, with the intervals (query coordinates) as breakpoint candidates.

' ) rows = ( - f'PHI testp = {signal.phi_p:.4g}' - f'{verdict} (alpha {alpha:g}; {signal.n_informative} informative ' - f'sites)' + f'PHI test{p_cell}' + f'{verdict}' f'Min recombination events (Rmin)' f'{signal.rmin}' f'{"intervals " + intervals if intervals else "none"}' @@ -512,7 +559,7 @@ def _site_track_html( if ctx.site_result is None: return "" fig = build_interactive_figure( - ctx.site_result, datasets, regions, ctx.gaps, + ctx.site_result, datasets, regions, ctx.caveat_gaps, y_title="Identity at informative sites", value_name="identity at informative sites", ) @@ -538,15 +585,24 @@ def write_html_report( output_dir: Path, logger: logging.Logger, ctx: ReportContext, + companion_files: list[str] | None = None, ) -> Path: - """Write a single self-contained ``report.html``.""" + """Write a single self-contained ``report.html``. + + ``companion_files`` are the names of the files written beside it, listed in the + footer; ``write_reports`` supplies them. + """ gaps = ctx.gaps + # Breakpoint gaps are tabulated but are not poorly covered stretches. + caveat_gaps = ctx.caveat_gaps lineage_map = ctx.lineage_map threshold = ctx.coverage_threshold - fig = build_interactive_figure(result, datasets, regions, gaps) + fig = build_interactive_figure(result, datasets, regions, caveat_gaps) plot_div = fig.to_html(full_html=False, include_plotlyjs="inline") - colors = _color_map(datasets) + # Colour every label a region names, not only the top-N: with --top-n 1 the donor + # is outside the plotted datasets and its swatch and mosaic segment were grey. + colors = _color_map(region_labels(datasets, regions)) s = _summary(result, regions, datasets) organism_html = ( f'
{html.escape(ctx.organism)}
' if ctx.organism else "" @@ -556,6 +612,11 @@ def write_html_report( for title, body in (ctx.extra_sections or []) ) + method_section = _method_section( + ctx.method_breakdown, ctx.methods_run, ctx.per_major, lineage_map, + ctx.methods_not_run, ctx.not_run_reason, + ) + doc = ( '\n\n' '\n' @@ -566,13 +627,13 @@ def write_html_report( f'

{html.escape(result.query)}

' f"{organism_html}" f"{_verdict_html(s, result.query, colors, lineage_map, ctx.query_lineage)}" - f"{_caveat_html(gaps, threshold)}" + f"{_caveat_html(caveat_gaps, threshold)}" f"{_cards_html(s, colors, lineage_map)}" '
Query mosaic
' - f"{_mosaic_html(regions, colors, s, gaps, lineage_map)}
" + f"{_mosaic_html(regions, colors, s, caveat_gaps, lineage_map)}" '
Recombinant regions
' f'{_regions_html(regions, colors, s["query_len"], lineage_map)}
' - f"{_method_section(ctx.method_breakdown, ctx.methods_run, ctx.per_major, lineage_map)}" + f"{method_section}" '
Reference coverage
' f"{_coverage_html(gaps, threshold)}
" f"{extras}" @@ -583,7 +644,7 @@ def write_html_report( f"{plot_div}" f'{_site_track_html(ctx, datasets, regions, provenance.get("windowing", ""))}' '
Recombination signal (parent-free)
' - f"{_signal_html(ctx.signal)}
" + f"{_signal_html(ctx.signal, ctx.alpha)}" '
Window winners
' '

Windows in which each reference is the query\'s closest match ' '(ties included).

' @@ -591,7 +652,7 @@ def write_html_report( '
Per-dataset similarity statistics
' f'{_stats_html(analysis, s["major"], lineage_map)}
' f'
{_methods_html(provenance)}
' - f"{_footer_html(provenance)}" + f"{_footer_html(provenance, companion_files)}" "" ) diff --git a/src/tessera/recomb/report_plots.py b/src/tessera/recomb/report_plots.py index 46ace49..a76b9e2 100644 --- a/src/tessera/recomb/report_plots.py +++ b/src/tessera/recomb/report_plots.py @@ -29,6 +29,25 @@ def _color_map(datasets: list[str]) -> dict[str, str]: return {label: PALETTE[i % len(PALETTE)] for i, label in enumerate(datasets)} +def region_labels(datasets: list[str], regions: list[Region]) -> list[str]: + """``datasets`` followed by any parent of a called region that is not among them. + + The plots draw the top-N datasets, but a region can name a parent outside that + list (``--top-n 1`` leaves out every donor). Colours are assigned from this + extended list so a region is never drawn in the fallback grey; the plotted + datasets come first, so their colours do not depend on which regions were called. + Donor-absent regions are skipped: their label is a stand-in, not a called parent. + """ + labels = list(datasets) + for region in regions: + if region.donor_absent: + continue + for label in (region.major_parent, region.minor_parent): + if label not in labels: + labels.append(label) + return labels + + def _palette(datasets: list[str]): """Matplotlib RGBA map drawn from the shared categorical palette.""" from matplotlib.colors import to_rgba @@ -46,12 +65,22 @@ def _ylim(values) -> tuple[float, float]: return lower, 1.005 +ABSENT_LABEL = "donor absent" + + def _shade_regions(ax, regions: list[Region], colors: dict) -> None: seen: set[str] = set() for region in regions: - color = colors.get(region.minor_parent, "grey") - label = f"recombinant: {region.minor_parent}" if region.minor_parent not in seen else None - seen.add(region.minor_parent) + if region.donor_absent: + # The stand-in label is the closest reference -- often the backbone itself -- + # so naming it here would read as "recombinant from the backbone". + key, color = ABSENT_LABEL, GREY + text = f"{ABSENT_LABEL} (no close reference)" + else: + key, color = region.minor_parent, colors.get(region.minor_parent, GREY) + text = f"recombinant: {region.minor_parent}" + label = text if key not in seen else None + seen.add(key) ax.axvspan(region.msa_start, region.msa_end, color=color, alpha=0.12, label=label) @@ -80,7 +109,7 @@ def plot_top_n( logger.warning("No datasets available for the top-N plot; skipping.") return None subset = df.loc[available] - colors = _palette(available) + colors = _palette(region_labels(available, regions)) fig, ax = plt.subplots(figsize=(14, 7)) _shade_regions(ax, regions, colors) @@ -121,7 +150,7 @@ def plot_pairwise( if seq1 not in df.index or seq2 not in df.index: logger.warning("Pairwise datasets not both present; skipping pairwise plot.") return None - colors = _palette([seq1, seq2]) + colors = _palette(region_labels([seq1, seq2], regions)) fig, ax = plt.subplots(figsize=(16, 6)) _shade_regions(ax, regions, colors) @@ -157,7 +186,7 @@ def build_interactive_figure( import plotly.graph_objects as go df = result.to_dataframe() - colors = _color_map(datasets) + colors = _color_map(region_labels(datasets, regions)) fig = go.Figure() for gap in coverage_gaps or []: fig.add_vrect( @@ -168,15 +197,18 @@ def build_interactive_figure( ) seen: set[str] = set() for region in regions: - color = colors.get(region.minor_parent, GREY) + if region.donor_absent: + key, color = ABSENT_LABEL, GREY + else: + key, color = region.minor_parent, colors.get(region.minor_parent, GREY) fig.add_vrect( x0=region.msa_start, x1=region.msa_end, fillcolor=color, opacity=0.12, line_width=0, layer="below", - annotation_text=("" if region.minor_parent in seen else region.minor_parent), + annotation_text=("" if key in seen else key), annotation_position="top left", annotation_font_size=11, annotation_font_color=color, ) - seen.add(region.minor_parent) + seen.add(key) for dataset in datasets: if dataset not in df.index: continue diff --git a/src/tessera/recomb/report_text.py b/src/tessera/recomb/report_text.py index 768eb7a..52a3953 100644 --- a/src/tessera/recomb/report_text.py +++ b/src/tessera/recomb/report_text.py @@ -12,7 +12,7 @@ from pathlib import Path from .analyze import AnalysisResult, stats_sort_key, winner_label -from .coverage import CoverageGap +from .coverage import BREAKPOINT_KIND, CoverageGap from .diagnostics import RecombinationSignal from .regions import Region from .similarity import WindowSimilarity @@ -150,8 +150,12 @@ def print_coverage(gaps: list[CoverageGap], threshold: float, echo=print) -> Non for g in gaps ] print_formatted_table(rows, header=COVERAGE_HEADER, echo=echo) - echo(" ^ the closest reference here is poor; the true source may be missing. " - "Run 'tessera find-references' to search NCBI.") + if any(g.kind != BREAKPOINT_KIND for g in gaps): + echo(" ^ the closest reference here is poor; the true source may be missing. " + "Run 'tessera find-references' to search NCBI.") + if any(g.kind == BREAKPOINT_KIND for g in gaps): + echo(" breakpoint = windows straddling a called breakpoint; the two parents " + "together explain the query there (not a missing reference).") echo("") @@ -275,9 +279,14 @@ def write_regions_tsv(regions: list[Region], output_dir: Path, logger: logging.L def write_methods_tsv( breakdown: list[dict], methods_run: tuple[str, ...], output_dir: Path, - logger: logging.Logger, + logger: logging.Logger, methods_not_run: tuple[str, ...] = (), ) -> None: - """Write the per-region x per-method agreement matrix (the ensemble breakdown).""" + """Write the per-region x per-method agreement matrix (the ensemble breakdown). + + A cell is ``yes`` / ``no`` for a caller that ran, and ``not run`` for one that was + selected but could not run -- a ``no`` there would read as a caller that looked and + found nothing. + """ path = output_dir / "recombination_methods.tsv" logger.info("Writing method comparison: %s", path) with open(path, "w") as fo: @@ -285,7 +294,10 @@ def write_methods_tsv( *methods_run, "parent_free_support"]) + "\n") for b in breakdown: called = b["per_method_support"] - cells = [("yes" if m in called else "no") for m in methods_run] + cells = [ + "not run" if m in methods_not_run else ("yes" if m in called else "no") + for m in methods_run + ] fo.write("\t".join(map(str, [ b["minor_parent"], b["query_start"], b["query_end"], *cells, "yes" if b["parent_free_support"] else "no", @@ -319,8 +331,12 @@ def write_profile_tsv( path = output_dir / "recombination_profile.tsv" logger.info("Writing recombination signal profile: %s", path) with open(path, "w") as fo: + # "NA" when the PHI test could not have rejected (too few informative sites + # for the window). The Rmin stays last in the note: the harness reads it there. + phi_p = "NA" if signal.phi_p is None else f"{signal.phi_p:.4g}" + note = "" if signal.phi_p is not None else "not testable at this window; " fo.write( - f"# PHI p-value\t{signal.phi_p:.4g}\t(window {signal.phi_window} " + f"# PHI p-value\t{phi_p}\t({note}window {signal.phi_window} " f"informative sites, {signal.n_informative} sites, Rmin {signal.rmin})\n" ) fo.write("msa_pos\tquery_pos\tphi\n") diff --git a/src/tessera/recomb/run.py b/src/tessera/recomb/run.py index 1dae03e..1cffc14 100644 --- a/src/tessera/recomb/run.py +++ b/src/tessera/recomb/run.py @@ -18,6 +18,7 @@ call_coverage_gaps, flag_undercovered_regions, gaps_as_regions, + mark_breakpoint_gaps, reconcile_gaps, ) from .diagnostics import corroborating_intervals, recombination_signal @@ -110,6 +111,31 @@ class RecombParams: reattribute_margin: float = 0.03 +def caller_description(method: str, params: RecombParams) -> str: + """One caller's name and the settings that govern it, for the run provenance. + + Every caller is described under its own name. The text is what a reader uses to + tell which tests produced the regions, so a caller must never be described as a + different one. + """ + alpha = f"alpha {params.alpha:g}" + if method == "hmm": + return f"hmm (jump-rate {params.jump_rate:g}, {alpha})" + if method == "3seq": + return f"3seq (triplet max-descent test, {alpha})" + if method == "maxchi": + return f"maxchi (chi-square triplet test, scan-aware permutation, {alpha})" + if method == "bootscan": + return f"bootscan (bootstrap support, block-permutation run-length test, {alpha})" + if method == "geneconv": + return f"geneconv (longest donor-match run, permutation test, {alpha})" + if method == "barcode": + return "barcode (clade-marker attribution on a typed panel; no significance test)" + min_region = params.min_region if params.min_region is not None else params.window_size + merge_gap = params.merge_gap if params.merge_gap is not None else params.window_size + return f"{method} (min {min_region} / margin {params.margin} / merge {merge_gap})" + + def _select_windowing(bp_result, params: RecombParams, query_label: str, logger): """Choose base-pair or informative-site windowing; return ``(result, label)``. @@ -281,7 +307,14 @@ def run_recomb( bp_result.rows, query_label, bp_result.column_to_query, window=params.phi_window, ) - if signal is not None: + if signal is not None and signal.phi_p is None: + logger.info( + "Recombination signal (parent-free): PHI not testable (every pair of the " + "%d informative site(s) lies within the window of %d ranks; lower " + "--phi-window), Rmin=%d.", + signal.n_informative, signal.phi_window, signal.rmin, + ) + elif signal is not None: logger.info( "Recombination signal (parent-free): PHI p=%.4g, Rmin=%d (%d informative " "sites).", signal.phi_p, signal.rmin, signal.n_informative, @@ -320,6 +353,30 @@ def _region_params(method: str) -> RegionParams: if method == "hmm": excluded_siblings = sibs + # The barcode caller needs typed references and returns no major parent when it + # cannot run. That is "could not test", which must not be reported as "tested and + # found nothing": refuse a run that selected nothing else, and say so otherwise. + not_run = tuple(m for m in params.methods if m == "barcode" and majors[m] is None) + not_run_reason = "" + if not_run: + reason = not_run_reason = ( + "fewer than two typed clades carry enough characteristic markers" + if lineage_map else + "it needs typed references (a lineage map: --lineage-map, or a lineages.tsv " + "beside the output or the MSA) and none was found" + ) + if len(not_run) == len(params.methods): + raise UserInputError( + f"The barcode caller could not run: {reason}. No other caller was " + "selected, so this scan could not test for recombination. Supply typed " + "references or choose another --method." + ) + logger.warning( + "The barcode caller could not run (%s); it is reported as 'not run', not as " + "a negative.", reason, + ) + n_ran = len(params.methods) - len(not_run) + major_parent, per_major = reconcile_major(majors, window_wins=analysis_bp.winners_with_ties) # Parent-free corroboration needs both halves of the diagnostic: PHI to establish # that the alignment carries recombination at all, the Rmin intervals to say where. @@ -334,7 +391,17 @@ def _region_params(method: str) -> RegionParams: # it (clamped, so selecting one caller is not silently self-suppressing). This runs # *before* re-attribution: a suppressed region should not be re-attributed, and must # not announce a re-attribution in the log for a region nobody will see. - min_agree = max(1, min(params.min_methods, len(params.methods))) + min_agree = max(1, min(params.min_methods, n_ran)) + # The clamp is silent when the user simply selected fewer callers than --min-methods. + # It is not when a selected caller could not run: the user asked for corroboration + # the run cannot give, and the regions reported are then weaker than requested. + gate_lowered = bool(not_run) and min_agree < params.min_methods + if gate_lowered: + logger.warning( + "--min-methods %d cannot be met: %d caller(s) ran. The agreement gate used " + "is %d, so regions found by fewer callers than requested are reported.", + params.min_methods, n_ran, min_agree, + ) regions, method_breakdown, suppressed = filter_by_agreement( regions, method_breakdown, min_agree ) @@ -347,6 +414,12 @@ def _region_params(method: str) -> RegionParams: regions, result, lineage_map, lineage_of(major_parent, lineage_map), margin=params.reattribute_margin, logger=logger, ) + # The breakdown rows were built from the regions before re-attribution and are + # written to recombination_methods.tsv and the report's method table. Keep the + # donor they name in step with the region. reattribute_donors returns one region + # per input region in the same order, so the two lists stay parallel. + for region, row in zip(regions, method_breakdown, strict=True): + row["minor_parent"] = region.minor_parent if excluded_siblings: logger.info( "Excluded %d whole-genome sibling(s) of the query (its own lineage) from the " @@ -376,11 +449,26 @@ def _region_params(method: str) -> RegionParams: bp_result, params.window_size, coverage_params ) flag_undercovered_regions(regions, coverage_threshold) - if coverage_gaps: + # A window straddling a called breakpoint matches neither parent well on its own. + # That is not a missing reference, so such gaps are relabelled before they can + # caveat a region or be bridged to a donor-absent one. Recruitment + # (fill-references / find-references) calls call_coverage_gaps directly and is + # unaffected. + n_breakpoint = mark_breakpoint_gaps( + coverage_gaps, regions, bp_result, params.window_size, coverage_threshold, + ) + n_poor = len(coverage_gaps) - n_breakpoint + if n_poor: logger.info( "Reference coverage: %d region(s) where the closest reference is below " "%.3f -- a better reference may be missing.", - len(coverage_gaps), coverage_threshold, + n_poor, coverage_threshold, + ) + if n_breakpoint: + logger.info( + "Reference coverage: %d low-similarity stretch(es) sit on a called breakpoint " + "and are explained by the two parents there; not treated as missing references.", + n_breakpoint, ) # Bridge: a divergent coverage gap (query far from every reference) is a @@ -404,16 +492,7 @@ def _region_params(method: str) -> RegionParams: print_regions(regions, major_parent, echo=logger.info) print_coverage(coverage_gaps, coverage_threshold, echo=logger.info) - def _caller_desc(method: str) -> str: - if method == "hmm": - return f"hmm (jump-rate {params.jump_rate:g}, alpha {params.alpha:g})" - if method == "3seq": - return f"3seq (triplet max-descent test, alpha {params.alpha:g})" - min_region = params.min_region if params.min_region is not None else params.window_size - merge_gap = params.merge_gap if params.merge_gap is not None else params.window_size - return f"heuristic (min {min_region} / margin {params.margin} / merge {merge_gap})" - - caller_desc = " + ".join(_caller_desc(m) for m in params.methods) + caller_desc = " + ".join(caller_description(m, params) for m in params.methods) provenance = { "tessera version": __version__, "date (UTC)": datetime.now(UTC).strftime("%Y-%m-%d %H:%M:%S"), @@ -426,16 +505,31 @@ def _caller_desc(method: str) -> str: "metric": params.metric, "caller": f"ensemble: {caller_desc}" if len(params.methods) > 1 else caller_desc, "windowing": windowing, + # Settings that change which regions are reported. Without them a run with + # --min-methods 2 or --no-cluster-lineages is indistinguishable from a default + # run in the record. + "min methods (agreement gate)": ( + f"{min_agree} (requested {params.min_methods}; {n_ran} caller" + f"{'' if n_ran == 1 else 's'} ran)" if gate_lowered else str(min_agree) + ), + "sibling exclusion": "on" if params.exclude_siblings else "off", + "lineage clustering": "on" if params.cluster_lineages else "off", + "donor re-attribution": ( + f"on (margin {params.reattribute_margin:g})" if params.reattribute_donors else "off" + ), "major parent": major_parent or "n/a", "coverage threshold / gaps": f"{coverage_threshold:.3f} / {len(coverage_gaps)}", } + if not_run: + provenance["callers not run"] = f"{', '.join(not_run)} ({not_run_reason})" if excluded_siblings: provenance["excluded siblings (query's own lineage)"] = ", ".join( ev.label for ev in excluded_siblings ) if signal is not None: + phi_text = "not testable" if signal.phi_p is None else f"p={signal.phi_p:.4g}" provenance["recombination signal (PHI)"] = ( - f"p={signal.phi_p:.4g} ({signal.n_informative} informative sites, " + f"{phi_text} ({signal.n_informative} informative sites, " f"window {signal.phi_window})" ) provenance["min recombination events (Rmin)"] = str(signal.rmin) @@ -468,6 +562,7 @@ def _caller_desc(method: str) -> str: extra_sections=extra_sections, lineage_map=lineage_map, query_lineage=query_lineage, signal=signal, organism=params.organism, methods_run=params.methods, method_breakdown=method_breakdown, per_major=per_major, + methods_not_run=not_run, not_run_reason=not_run_reason, alpha=params.alpha, ), ) logger.info("All done.") diff --git a/tests/unit/test_barcode.py b/tests/unit/test_barcode.py index 73b3731..cc4f4a9 100644 --- a/tests/unit/test_barcode.py +++ b/tests/unit/test_barcode.py @@ -2,16 +2,22 @@ from __future__ import annotations +import csv +import json +import logging from pathlib import Path import numpy as np +import pytest +from tessera.core.errors import UserInputError from tessera.recomb.analyze import analyze from tessera.recomb.barcode import clade_markers from tessera.recomb.regions import RegionParams, call_regions +from tessera.recomb.run import RecombParams, run_recomb from tessera.recomb.similarity import compute_similarity -from ..conftest import write_fasta +from ..conftest import recombinant_msa, write_fasta L = 3000 LMAP = {"a1": "A", "a2": "A", "b1": "B", "b2": "B", "c1": "C", "c2": "C"} @@ -84,3 +90,101 @@ def test_barcode_silent_on_untyped_panel(tmp_path: Path) -> None: params = RegionParams.with_defaults(300, method="barcode") # no lineage map regions, _, _ = call_regions(result, analyze(result), 300, params) assert regions == [] + + +# --- "could not run" is not "found nothing" -------------------------------- + +def test_barcode_names_no_major_parent_when_it_cannot_run(tmp_path: Path) -> None: + """It used to return the first record of the alignment as the major parent, which a + barcode-only run then reported as the backbone.""" + result = compute_similarity(str(_panel(tmp_path)), "query", window_size=300, window_step=30) + params = RegionParams.with_defaults(300, method="barcode") # no lineage map + regions, major, _ = call_regions(result, analyze(result), 300, params) + assert regions == [] + assert major is None + + +def _run(tmp_path: Path, logger, **kwargs) -> Path: + out = tmp_path / "out" + run_recomb( + RecombParams(msa=recombinant_msa(tmp_path, recombinant=True), output=out, + query="query", plot_format="png", **kwargs), + logger, + ) + return out + + +def test_barcode_only_run_on_an_untyped_panel_is_refused(tmp_path: Path, logger) -> None: + """Nothing was tested, so there is no result to report -- not a clean negative.""" + with pytest.raises(UserInputError, match="typed references"): + _run(tmp_path, logger, methods=("barcode",)) + assert not (tmp_path / "out" / "report.html").exists() + + +def _collecting_logger() -> tuple[logging.Logger, list[str]]: + """A logger outside the ``tessera`` hierarchy that keeps its warnings in a list. + + Not ``caplog``: once a CLI test has configured the ``tessera`` logger it stops + propagating to the root logger, and whether that has happened depends on test order. + """ + messages: list[str] = [] + + class _Collect(logging.Handler): + def emit(self, record: logging.LogRecord) -> None: + messages.append(record.getMessage()) + + log = logging.getLogger("test_barcode.collect") + log.handlers = [_Collect(level=logging.WARNING)] + log.propagate = False + return log, messages + + +def test_ensemble_reports_barcode_as_not_run_on_an_untyped_panel(tmp_path: Path) -> None: + log, warnings = _collecting_logger() + out = _run(tmp_path, log, methods=("3seq", "barcode")) + assert any("barcode" in message and "could not run" in message for message in warnings) + + rows = list(csv.DictReader((out / "recombination_methods.tsv").open(), delimiter="\t")) + assert rows, "3seq should still call the insert" + assert {row["3seq"] for row in rows} == {"yes"} + assert {row["barcode"] for row in rows} == {"not run"} # was "no" + + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert "barcode" in run["callers not run"] + assert "not run" in (out / "report.html").read_text() + + +def test_agreement_gate_counts_only_the_callers_that_ran(tmp_path: Path, logger) -> None: + """--min-methods is clamped to the callers actually run. A caller that could not run + cannot agree, so it must not make the gate unreachable.""" + out = _run(tmp_path, logger, methods=("3seq", "barcode"), min_methods=2) + rows = list(csv.DictReader((out / "recombination_regions.tsv").open(), delimiter="\t")) + assert [row["methods"] for row in rows] == ["3seq"] + + +def test_a_lowered_agreement_gate_is_logged_and_recorded(tmp_path: Path) -> None: + """The user asked for two agreeing callers and only one could run. The regions are + then single-caller regions; the log and the provenance must say the gate was 1.""" + log, warnings = _collecting_logger() + out = _run(tmp_path, log, methods=("3seq", "barcode"), min_methods=2) + assert any("--min-methods 2" in m and "1 caller" in m for m in warnings) + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert run["min methods (agreement gate)"] == "1 (requested 2; 1 caller ran)" + + +def test_barcode_not_run_on_a_typed_panel_without_marked_clades(tmp_path: Path) -> None: + """Typed, but every reference is in one clade: there are no clade markers to + compete, so the caller cannot run and the warning says why.""" + log, warnings = _collecting_logger() + one_clade = {"A": "L1", "B": "L1", "other": "L1"} + out = _run(tmp_path, log, methods=("3seq", "barcode"), lineage_map=one_clade) + assert any("fewer than two typed clades" in message for message in warnings) + rows = list(csv.DictReader((out / "recombination_methods.tsv").open(), delimiter="\t")) + assert {row["barcode"] for row in rows} == {"not run"} + # The panel is typed, so "needs typed references" would be the wrong reason. + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert "fewer than two typed clades" in run["callers not run"] + assert "needs typed references" not in run["callers not run"] + report = (out / "report.html").read_text() + assert "fewer than two typed clades" in report + assert "none were available" not in report diff --git a/tests/unit/test_binaries.py b/tests/unit/test_binaries.py index 8cbeca8..8be4617 100644 --- a/tests/unit/test_binaries.py +++ b/tests/unit/test_binaries.py @@ -147,3 +147,45 @@ def test_a_binary_that_cannot_execute_is_not_a_crash(on_path) -> None: pytest.skip("file is still executable on this platform") with pytest.raises(MissingBinaryError): check_binaries((BinarySpec("broken"),)) + + +# --- a failed probe is not a version -------------------------------------- + +def failing_binary(directory: Path, name: str, output: str) -> Path: + """An executable that prints ``output`` to stderr and exits 1, like a tool given an + option it does not have.""" + path = directory / name + path.write_text(f'#!/bin/sh\nprintf %s "{output}" >&2\nexit 1\n') + path.chmod(path.stat().st_mode | stat.S_IEXEC | stat.S_IXGRP | stat.S_IXOTH) + return path + + +def test_error_output_is_not_recorded_as_a_version(on_path) -> None: + """`sibeliaz -v` prints "illegal option -- v" and exits 1; that line went into the + provenance sidecar as the aligner version.""" + failing_binary(on_path, "sibeliaz", "sibeliaz: illegal option -- v") + versions = check_binaries((BinarySpec("sibeliaz", version_args=("-v",)),)) + assert versions["sibeliaz"] == "unknown" + + +def test_failed_probe_that_still_prints_a_version_is_kept(on_path) -> None: + # Some tools print their version in a usage message and exit non-zero. + failing_binary(on_path, "usagey", "usagey 1.4.2 -- usage: usagey [options]") + assert check_binaries((BinarySpec("usagey"),))["usagey"] == "1.4.2" + + +def test_tool_without_a_version_option_is_not_probed(on_path, tmp_path) -> None: + marker = tmp_path / "was_run" + path = on_path / "noversion" + path.write_text(f'#!/bin/sh\ntouch "{marker}"\n') + path.chmod(path.stat().st_mode | stat.S_IEXEC | stat.S_IXGRP | stat.S_IXOTH) + versions = check_binaries((BinarySpec("noversion", version_args=None),)) + assert versions == {"noversion": "unknown"} + assert not marker.exists() # declared as having no version option: never executed + + +def test_sibeliaz_declares_no_version_probe() -> None: + from tessera.aligners.sibeliaz import SibeliazAligner + + (spec,) = SibeliazAligner.capabilities.required_binaries + assert spec.version_args is None diff --git a/tests/unit/test_breakpoint_gaps.py b/tests/unit/test_breakpoint_gaps.py new file mode 100644 index 0000000..f609af1 --- /dev/null +++ b/tests/unit/test_breakpoint_gaps.py @@ -0,0 +1,224 @@ +"""A coverage gap produced by a window straddling a called breakpoint is not a +missing reference. + +A window that spans a breakpoint between two divergent parents matches neither parent +well on its own, so its best similarity dips below the coverage threshold even though +both parents are in the panel. Such a gap is relabelled ``breakpoint``; it must not +caveat the region, become a donor-absent region, or reach the report headline. +""" + +from __future__ import annotations + +import csv +import re + +import numpy as np + +from tessera.recomb.coverage import ( + BREAKPOINT_KIND, + CoverageGap, + gaps_as_regions, + mark_breakpoint_gaps, +) +from tessera.recomb.regions import Region +from tessera.recomb.run import RecombParams, run_recomb + +WIDTH = 2000 +WINDOW = 400 +STEP = 40 + + +def _enc(seq: str) -> np.ndarray: + return np.frombuffer(seq.encode("ascii"), dtype=np.uint8) + + +_MINOR = "".join("C" if i % 8 == 0 else "A" for i in range(WIDTH)) + + +class _Scan: + """The part of a WindowSimilarity that gap marking reads: the aligned rows and, per + window, its centre and the best single-reference identity.""" + + def __init__(self, query: str) -> None: + self.query = "query" + self.rows = {"query": _enc(query), "major": _enc("A" * WIDTH), "minor": _enc(_MINOR)} + self.positions = [s + WINDOW // 2 for s in range(0, WIDTH - WINDOW + 1, STEP)] + q = self.rows["query"] + self.best_sim = [ + max(float(np.mean(q[p - WINDOW // 2:p + WINDOW // 2] + == self.rows[ref][p - WINDOW // 2:p + WINDOW // 2])) + for ref in ("major", "minor")) + for p in self.positions + ] + + +def _rows(query: str) -> _Scan: + """Parents that differ at every eighth column; the query is given by the caller.""" + return _Scan(query) + + +def _mosaic() -> str: + """Major up to column 1000, minor after it: a clean breakpoint at 1000.""" + return "A" * 1000 + _MINOR[1000:] + + +def _region() -> Region: + return Region( + minor_parent="minor", major_parent="major", msa_start=1000, msa_end=WIDTH, + query_start=1000, query_end=WIDTH, n_windows=10, + mean_sim_minor=1.0, mean_sim_major=0.9, margin=0.1, + ) + + +def _gap(start: int, end: int, kind: str = "divergent") -> CoverageGap: + return CoverageGap( + msa_start=start, msa_end=end, query_start=start, query_end=end, + length_bp=end - start, n_windows=2, best_label="major", mean_best=0.94, kind=kind, + ) + + +def test_gap_straddling_a_called_breakpoint_is_relabelled() -> None: + gaps = [_gap(800, 1200)] + n = mark_breakpoint_gaps(gaps, [_region()], _rows(_mosaic()), WINDOW, 0.95) + assert n == 1 + assert gaps[0].kind == BREAKPOINT_KIND + + +def test_breakpoint_gap_is_not_bridged_to_a_donor_absent_region() -> None: + gaps = [_gap(800, 1200)] + mark_breakpoint_gaps(gaps, [_region()], _rows(_mosaic()), WINDOW, 0.95) + + class _Result: + positions: list[int] = [] + similarities: dict[str, list[float]] = {} + + assert gaps_as_regions(gaps, _Result(), "major") == [] + + +def test_gap_the_two_parents_do_not_explain_stays_divergent() -> None: + # Next to the breakpoint, but the query carries a base neither parent has at one + # column in five: a third source, not a straddling window. + query = list(_mosaic()) + for i in range(800, 1200, 5): + query[i] = "G" + gaps = [_gap(800, 1200)] + n = mark_breakpoint_gaps(gaps, [_region()], _rows("".join(query)), WINDOW, 0.95) + assert n == 0 + assert gaps[0].kind == "divergent" + + +def test_gap_far_from_every_region_boundary_stays_divergent() -> None: + gaps = [_gap(100, 500)] # more than one window away from 1000 and from 2000 + n = mark_breakpoint_gaps(gaps, [_region()], _rows(_mosaic()), WINDOW, 0.95) + assert n == 0 + assert gaps[0].kind == "divergent" + + +def test_low_information_gap_is_left_alone() -> None: + gaps = [_gap(800, 1200, kind="low_information")] + mark_breakpoint_gaps(gaps, [_region()], _rows(_mosaic()), WINDOW, 0.95) + assert gaps[0].kind == "low_information" + + +def test_donor_absent_region_does_not_explain_a_gap() -> None: + region = _region() + region.donor_absent = True + gaps = [_gap(800, 1200)] + assert mark_breakpoint_gaps(gaps, [region], _rows(_mosaic()), WINDOW, 0.95) == 0 + + +def test_region_naming_a_label_outside_the_alignment_is_skipped() -> None: + # A region can carry a placeholder parent ("n/a" when no backbone was resolved). + region = _region() + region.major_parent = "n/a" + gaps = [_gap(800, 1200)] + assert mark_breakpoint_gaps(gaps, [region], _rows(_mosaic()), WINDOW, 0.95) == 0 + assert gaps[0].kind == "divergent" + + +def test_short_third_source_stretch_inside_a_wider_gap_stays_divergent() -> None: + # 120 columns right after the breakpoint carry a base neither parent has at one + # column in four (30 columns). Averaged over the whole 640-column gap the two parents + # still explain 95.3 % of it, but a 400-column window that holds the stretch is + # explained to 92.5 % at best, whatever the switch point. + query = list(_mosaic()) + for i in range(1000, 1120, 4): + query[i] = "G" + gaps = [_gap(680, 1320)] + n = mark_breakpoint_gaps(gaps, [_region()], _rows("".join(query)), WINDOW, 0.95) + assert n == 0 + assert gaps[0].kind == "divergent" + + +def test_gap_without_an_under_threshold_window_is_not_relabelled() -> None: + # Nothing in the gap is below the threshold in this scan, so there is no straddling + # window to explain; leave the label alone. + gaps = [_gap(100, 500)] + region = _region() + region.msa_start = 300 + n = mark_breakpoint_gaps(gaps, [region], _rows("A" * WIDTH), WINDOW, 0.95) + assert n == 0 + + +# --- the shipped example, end to end -------------------------------------- + +def test_clean_recombinant_is_not_reported_as_a_missing_reference( + example_data, tmp_path, logger +) -> None: + """`example_data/divergent.msa.fasta`: the donor matches the query exactly over the + insert and four callers agree. The only coverage gaps are the windows straddling + the two breakpoints.""" + out = tmp_path / "out" + run_recomb( + RecombParams(msa=example_data / "divergent.msa.fasta", output=out, query="query", + window_size=300, window_step=30, plot_format="png"), + logger, + ) + regions = list(csv.DictReader((out / "recombination_regions.tsv").open(), delimiter="\t")) + assert len(regions) == 1 + assert regions[0]["donor_undercovered"] == "no" + assert regions[0]["donor_absent"] == "no" + + lines = (out / "coverage_gaps.tsv").read_text().splitlines() + gaps = list(csv.DictReader([ln for ln in lines if not ln.startswith("#")], delimiter="\t")) + assert gaps and {g["kind"] for g in gaps} == {BREAKPOINT_KIND} + + report = (out / "report.html").read_text() + verdict = re.sub(r"<[^>]+>", "", re.search(r'

.*?

', report).group(0)) + assert "high confidence" in " ".join(verdict.split()) + assert "Possible missing reference" not in report + + +def test_unsampled_stretch_beside_the_breakpoint_keeps_its_caveat(tmp_path, logger) -> None: + """A clean insert followed by 150 bp from a lineage that is not in the panel (about + 15 % from both parents). The coverage gap there is deeper than a straddling window + explains, so it must stay `divergent` and keep caveating the region.""" + import random + + rng = random.Random(2) + + def mutate(seq: list[str], rate: float) -> list[str]: + return [rng.choice([b for b in "ACGT" if b != c]) if rng.random() < rate else c + for c in seq] + + anc = [rng.choice("ACGT") for _ in range(3000)] + a, b, og = mutate(anc, 0.055), mutate(anc, 0.055), mutate(anc, 0.15) + query = a[:1000] + b[1000:2000] + mutate(a[2000:2150], 0.15) + a[2150:] + msa = tmp_path / "third_source.fasta" + msa.write_text("".join( + f">{name}\n{''.join(seq)}\n" + for name, seq in (("query", query), ("parent_A", a), ("parent_B", b), ("outgroup", og)) + )) + out = tmp_path / "out" + run_recomb( + RecombParams(msa=msa, output=out, query="query", window_size=300, window_step=30, + plot_format="png"), + logger, + ) + lines = (out / "coverage_gaps.tsv").read_text().splitlines() + gaps = list(csv.DictReader([ln for ln in lines if not ln.startswith("#")], delimiter="\t")) + over_stretch = [g for g in gaps if int(g["msa_start"]) < 2150 and int(g["msa_end"]) > 2000] + assert over_stretch and {g["kind"] for g in over_stretch} == {"divergent"} + regions = list(csv.DictReader((out / "recombination_regions.tsv").open(), delimiter="\t")) + donor = [r for r in regions if r["minor_parent"] == "parent_B" and r["donor_absent"] == "no"] + assert donor and donor[0]["donor_undercovered"] == "yes" diff --git a/tests/unit/test_caller_provenance.py b/tests/unit/test_caller_provenance.py new file mode 100644 index 0000000..83d13f1 --- /dev/null +++ b/tests/unit/test_caller_provenance.py @@ -0,0 +1,68 @@ +"""The run provenance names the callers that ran and the settings that shape the result. + +`run_provenance.json` and the report's "Run parameters" table are what a reader uses to +reproduce a run. A caller described under another caller's name, or a setting that +changes which regions are reported but is not recorded, makes the record misleading. +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import pytest + +from tessera.recomb.regions import CALLERS +from tessera.recomb.run import RecombParams, caller_description, run_recomb + + +def _params(example_data: Path, output: Path, **kwargs) -> RecombParams: + return RecombParams( + msa=example_data / "divergent.msa.fasta", output=output, query="query", + window_size=300, window_step=30, plot_format="png", **kwargs, + ) + + +@pytest.mark.parametrize("method", CALLERS) +def test_each_caller_is_described_under_its_own_name(example_data, tmp_path, method) -> None: + description = caller_description(method, _params(example_data, tmp_path)) + assert description.startswith(method) + + +@pytest.mark.parametrize("method", [m for m in CALLERS if m != "heuristic"]) +def test_only_the_heuristic_caller_is_called_heuristic(example_data, tmp_path, method) -> None: + # maxchi, bootscan, geneconv and barcode used to fall through to the heuristic text. + assert "heuristic" not in caller_description(method, _params(example_data, tmp_path)) + + +def test_default_run_names_its_four_callers(example_data, tmp_path, logger) -> None: + out = tmp_path / "out" + run_recomb(_params(example_data, out), logger) + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert "heuristic" not in run["caller"] + for name in ("hmm", "3seq", "maxchi", "bootscan"): + assert name in run["caller"] + + +def test_provenance_records_the_settings_that_change_what_is_reported( + example_data, tmp_path, logger +) -> None: + out = tmp_path / "default" + run_recomb(_params(example_data, out), logger) + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert run["min methods (agreement gate)"] == "1" + assert run["sibling exclusion"] == "on" + assert run["lineage clustering"] == "on" + assert run["donor re-attribution"] == "off" + + out = tmp_path / "changed" + run_recomb( + _params(example_data, out, min_methods=2, exclude_siblings=False, + cluster_lineages=False, reattribute_donors=True, reattribute_margin=0.05), + logger, + ) + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert run["min methods (agreement gate)"] == "2" + assert run["sibling exclusion"] == "off" + assert run["lineage clustering"] == "off" + assert run["donor re-attribution"] == "on (margin 0.05)" diff --git a/tests/unit/test_cli_input_checks.py b/tests/unit/test_cli_input_checks.py new file mode 100644 index 0000000..4afded8 --- /dev/null +++ b/tests/unit/test_cli_input_checks.py @@ -0,0 +1,249 @@ +"""Inputs the CLI must reject before any work starts. + +Each case used to be accepted: the option was silently ignored, or the run failed +later under "Unexpected error". None of these tests may reach the network -- the +panel-building entry points are replaced with a function that fails the test if the +command gets that far. +""" + +from __future__ import annotations + +from pathlib import Path + +import pytest +from typer.testing import CliRunner + +from tessera.cli.main import app + +runner = CliRunner() + +EXAMPLE = "example_data/divergent.msa.fasta" + + +def _must_not_run(*args, **kwargs): + raise AssertionError("the command reached the pipeline; it should have been rejected") + + +@pytest.fixture +def no_pipeline(monkeypatch): + """Fail the test if a command gets past validation into panel building or typing.""" + monkeypatch.setattr("tessera.discover.iterate.fill_references", _must_not_run) + monkeypatch.setattr("tessera.discover.lineage_assign.assign_lineages", _must_not_run) + monkeypatch.setattr("tessera.discover.run.find_references", _must_not_run) + + +def _query(tmp_path: Path) -> Path: + path = tmp_path / "query.fasta" + path.write_text(">query\nACGTACGTACGT\n") + return path + + +def _collection(tmp_path: Path) -> Path: + directory = tmp_path / "collection" + directory.mkdir() + (directory / "ref.fasta").write_text(">ref\nACGTACGTACGT\n") + return directory + + +def _rejected(result, needle: str) -> None: + assert result.exit_code == 1, result.output + assert "Unexpected error" not in result.output + assert needle in " ".join(result.output.split()) + + +# --- --lineage-map must exist --------------------------------------------- + +def test_recomb_rejects_a_missing_lineage_map(tmp_path: Path) -> None: + """A mistyped path used to be ignored: the run exited 0 with an untyped report, and + the barcode caller and donor re-attribution silently did nothing.""" + result = runner.invoke(app, [ + "recomb", "--msa", EXAMPLE, "--query", "query", + "--output", str(tmp_path / "out"), "--lineage-map", str(tmp_path / "nope.tsv"), + ]) + _rejected(result, "--lineage-map file not found") + assert not (tmp_path / "out" / "report.html").exists() + + +def test_recomb_rejects_a_lineage_map_that_is_a_directory(tmp_path: Path) -> None: + result = runner.invoke(app, [ + "recomb", "--msa", EXAMPLE, "--query", "query", + "--output", str(tmp_path / "out"), "--lineage-map", str(tmp_path), + ]) + _rejected(result, "is a directory, not a file") + + +def test_type_lineages_rejects_a_missing_lineage_map(tmp_path: Path, no_pipeline) -> None: + result = runner.invoke(app, [ + "type-lineages", "--collection", str(_collection(tmp_path)), + "--output", str(tmp_path / "out"), "--lineage-map", str(tmp_path / "nope.tsv"), + ]) + _rejected(result, "--lineage-map file not found") + + +@pytest.mark.parametrize("command", ["detect", "fill-references", "build-panel"]) +def test_panel_commands_reject_a_missing_lineage_map( + tmp_path: Path, no_pipeline, command: str +) -> None: + result = runner.invoke(app, [ + command, "--query", str(_query(tmp_path)), "--output", str(tmp_path / "out"), + "--lineage-map", str(tmp_path / "nope.tsv"), + ]) + _rejected(result, "--lineage-map file not found") + + +# --- paths ------------------------------------------------------------------ + +def test_find_references_rejects_a_missing_msa(tmp_path: Path, no_pipeline) -> None: + result = runner.invoke(app, [ + "find-references", "--msa", str(tmp_path / "nope.fasta"), "--query", "query", + "--output", str(tmp_path / "out"), + ]) + _rejected(result, "MSA file not found") + + +def test_recomb_rejects_an_output_path_that_is_a_file(tmp_path: Path) -> None: + occupied = tmp_path / "out" + occupied.write_text("not a directory\n") + result = runner.invoke(app, [ + "recomb", "--msa", EXAMPLE, "--query", "query", "--output", str(occupied), + ]) + _rejected(result, "is not a directory") + assert occupied.read_text() == "not a directory\n" # left untouched + + +# --- numeric options on the panel-building commands ------------------------ + +@pytest.mark.parametrize("command", ["detect", "fill-references", "build-panel"]) +@pytest.mark.parametrize( + ("option", "value"), + [("--max-rounds", "0"), ("--window-size", "0"), ("--window-step", "0")], +) +def test_panel_commands_reject_out_of_range_options( + tmp_path: Path, no_pipeline, command: str, option: str, value: str +) -> None: + """`fill-references --max-rounds 0` exited 0 having built no alignment and no + report; a bad window reached the scan only after seeding and alignment.""" + result = runner.invoke(app, [ + command, "--query", str(_query(tmp_path)), "--output", str(tmp_path / "out"), + option, value, + ]) + _rejected(result, f"Invalid {option}") + + +@pytest.mark.parametrize(("option", "value"), [("--window-size", "0"), ("--window-step", "0")]) +def test_find_references_rejects_out_of_range_windows( + tmp_path: Path, no_pipeline, option: str, value: str +) -> None: + result = runner.invoke(app, [ + "find-references", "--msa", EXAMPLE, "--query", "query", + "--output", str(tmp_path / "out"), option, value, + ]) + _rejected(result, f"Invalid {option}") + + +# --- reassort --------------------------------------------------------------- + +def _segments(tmp_path: Path) -> Path: + path = tmp_path / "segments.fasta" + path.write_text(">HA\nACGTACGT\n>NA\nACGTACGT\n") + return path + + +@pytest.fixture +def no_reassort(monkeypatch): + monkeypatch.setattr("tessera.cli.cmd_reassort.assign_segments", _must_not_run) + + +@pytest.mark.parametrize( + ("option", "value", "needle"), + [("--ani-floor", "500", "Invalid --ani-floor"), + ("--ani-floor", "-1", "Invalid --ani-floor"), + ("--margin", "-3", "Invalid --margin")], +) +def test_reassort_rejects_out_of_range_options( + tmp_path: Path, no_reassort, option: str, value: str, needle: str +) -> None: + """A negative margin makes every near-best set empty (a clonal pair then reads as + undetermined); an ANI floor above 100 leaves every segment unassigned. Both exited 0.""" + result = runner.invoke(app, [ + "reassort", "--query", str(_segments(tmp_path)), "--output", str(tmp_path / "out"), + option, value, + ]) + _rejected(result, needle) + + +def test_reassort_rejects_a_dataset_override_for_an_unknown_segment( + tmp_path: Path, no_reassort +) -> None: + """A typo in the segment name silently fell back to dataset auto-detection.""" + result = runner.invoke(app, [ + "reassort", "--query", str(_segments(tmp_path)), "--output", str(tmp_path / "out"), + "--dataset", "HAA=nextstrain/flu/h3n2/ha", + ]) + _rejected(result, "HAA") + assert "HA, NA" in " ".join(result.output.split()) # names the segments that do exist + + +def test_segment_scan_tsv_keeps_one_row_per_segment(tmp_path: Path, monkeypatch) -> None: + """An aligner error message spans several lines and may hold tabs; written verbatim + into the note column it broke the table.""" + from tessera.reassort.assign import ReassortmentResult, SegmentAssignment + from tessera.reassort.scan import SegmentScan + + def fake_assign(query, **kwargs): + return ReassortmentResult( + segments=[SegmentAssignment("HA", "dataset", None, None, 0.0, "unassigned")], + verdict="undetermined", groups=[], pair_notes=[], + scans=[SegmentScan("HA", False, False, 0, + "scan failed: mafft failed:\nline1\tx\nline2")], + ) + + monkeypatch.setattr("tessera.cli.cmd_reassort.assign_segments", fake_assign) + out = tmp_path / "out" + result = runner.invoke(app, [ + "reassort", "--query", str(_segments(tmp_path)), "--output", str(out), + ]) + assert result.exit_code == 0, result.output + lines = (out / "segment_scan.tsv").read_text().splitlines() + assert len(lines) == 2 # header + one segment + assert all(len(line.split("\t")) == 4 for line in lines) + assert lines[1].split("\t")[3] == "scan failed: mafft failed: line1 x line2" + + +# --- type-lineages reads the same collection the other commands do --------- + +def test_type_lineages_accepts_what_a_collection_may_hold(tmp_path: Path, monkeypatch) -> None: + """Every other command treats each file in the collection as a genome, whatever its + extension and gzip-compressed or not; type-lineages filtered on three suffixes and + reported "No FASTA genomes found".""" + import gzip + + collection = tmp_path / "collection" + collection.mkdir() + with gzip.open(collection / "a.fasta.gz", "wt") as fo: + fo.write(">a\nACGTACGT\n") + (collection / "b.fas").write_text(">b\nACGTACGT\n") + (collection / ".DS_Store").write_text("not a genome") + seen: list[str] = [] + + def fake_assign(genomes, **kwargs): + seen.extend(p.name for p in genomes) + return [("a", "L1", "denovo"), ("b.fas", "L1", "denovo")] + + monkeypatch.setattr("tessera.discover.lineage_assign.assign_lineages", fake_assign) + out = tmp_path / "out" + result = runner.invoke(app, [ + "type-lineages", "--collection", str(collection), "--output", str(out), + ]) + assert result.exit_code == 0, result.output + assert seen == ["a.fasta.gz", "b.fas"] # hidden files are not genomes + assert (out / "lineages.tsv").exists() + + +def test_type_lineages_rejects_a_file_that_is_not_fasta(tmp_path: Path, no_pipeline) -> None: + collection = _collection(tmp_path) + (collection / "notes.txt").write_text("these are my notes\n") + result = runner.invoke(app, [ + "type-lineages", "--collection", str(collection), "--output", str(tmp_path / "out"), + ]) + _rejected(result, "does not look like a FASTA file") diff --git a/tests/unit/test_dependencies.py b/tests/unit/test_dependencies.py new file mode 100644 index 0000000..5daabf0 --- /dev/null +++ b/tests/unit/test_dependencies.py @@ -0,0 +1,28 @@ +"""Every declared runtime dependency is one the package imports. + +Tessera is dependency-light by design, so a declared dependency nothing imports is a +cost with no benefit: it is installed for every user and audited in CI. +""" + +from __future__ import annotations + +import re +import tomllib +from pathlib import Path + +REPO = Path(__file__).resolve().parents[2] + +# Distribution name -> import name, where they differ. +IMPORT_NAME = {"biopython": "Bio"} + + +def test_every_runtime_dependency_is_imported() -> None: + project = tomllib.loads((REPO / "pyproject.toml").read_text())["project"] + names = [re.split(r"[<>=!~\[ ;]", dep, maxsplit=1)[0] for dep in project["dependencies"]] + source = "\n".join(p.read_text() for p in (REPO / "src" / "tessera").rglob("*.py")) + unused = [ + name for name in names + if not re.search(rf"^\s*(?:import|from)\s+{IMPORT_NAME.get(name, name)}\b", + source, flags=re.M) + ] + assert unused == [] diff --git a/tests/unit/test_diagnostics.py b/tests/unit/test_diagnostics.py index c5450de..1fb0e35 100644 --- a/tests/unit/test_diagnostics.py +++ b/tests/unit/test_diagnostics.py @@ -151,3 +151,27 @@ def test_corroborating_intervals_require_a_significant_phi() -> None: def test_corroborating_intervals_without_a_signal() -> None: assert corroborating_intervals(None, alpha=0.05) == [] + + +# --- "not testable" is not "no signal" -------------------------------------- + +def test_phi_is_not_testable_when_every_site_pair_is_inside_the_window() -> None: + """With z informative sites and a window of at least z - 1 ranks, the statistic + averages over every pair of sites, so reordering the sites cannot change it and the + permutation p-value is 1 whatever the data. That was reported as "no signal".""" + rows = _block_alignment(recombinant=True) # 24 informative columns + signal = recombination_signal(rows, "s0", lambda c: c, window=100, seed=1) + assert signal is not None + assert signal.n_informative == 24 + assert signal.phi_p is None # was 1.0 + assert signal.rmin >= 1 # Rmin does not depend on the window + assert corroborating_intervals(signal, alpha=0.05) == [] + + +def test_phi_becomes_testable_one_rank_below_the_site_count() -> None: + rows = _block_alignment(recombinant=True) # z = 24, so z - 1 = 23 + at_limit = recombination_signal(rows, "s0", lambda c: c, window=23, seed=1) + below = recombination_signal(rows, "s0", lambda c: c, window=22, seed=1) + assert at_limit is not None and below is not None + assert at_limit.phi_p is None + assert below.phi_p is not None diff --git a/tests/unit/test_harness_scoring.py b/tests/unit/test_harness_scoring.py index 832db8f..9125033 100644 --- a/tests/unit/test_harness_scoring.py +++ b/tests/unit/test_harness_scoring.py @@ -210,3 +210,18 @@ def test_mask_sibling_fail_on_sibling_donor(tmp_path): si = _setup(out=tmp_path, case_type="single_insert", clade_a="A", clade_b="B.1", q_start=100, q_end=200) assert rh._score_single_insert(tmp_path, clade_of, si, 5, "tip", 1.0)["pass"] is True + + +def test_parse_signal_reads_an_untestable_phi_header(tmp_path): + """The profile header carries NA when the PHI test could not have rejected; the + harness must still find the Rmin that follows it.""" + import logging + + from tessera.recomb.diagnostics import RecombinationSignal + from tessera.recomb.report_text import write_profile_tsv + + signal = RecombinationSignal( + n_informative=54, phi_p=None, phi_observed=0.1, phi_window=100, rmin=2, + ) + write_profile_tsv(signal, tmp_path, logging.getLogger("tessera")) + assert rh.parse_signal(tmp_path / "recombination_profile.tsv") == ("NA", "2") diff --git a/tests/unit/test_reattribute.py b/tests/unit/test_reattribute.py index 7f9de7d..b18631b 100644 --- a/tests/unit/test_reattribute.py +++ b/tests/unit/test_reattribute.py @@ -2,10 +2,17 @@ from __future__ import annotations +import csv +import random +from pathlib import Path + import numpy as np from tessera.recomb.reattribute import reattribute_donors from tessera.recomb.regions import Region +from tessera.recomb.run import RecombParams, run_recomb + +from ..conftest import write_fasta class _Result: @@ -103,3 +110,50 @@ def test_keeps_an_untyped_donor_that_cannot_be_scored(): out = reattribute_donors([_region("x1", "a1", 5, 15)], _Result(rows, "q"), lm, "A", margin=0.1, min_sites=4) assert out[0].minor_parent == "x1" + + +# --- the run's other outputs follow the re-attribution --------------------- + +def _reattribution_panel(tmp_path: Path) -> tuple[Path, dict[str, str]]: + """A cowpox backbone with a variola insert at 2000-4000. The insert is copied from + ``variolaA``, but ``variolaA`` shares clade V with two distant genomes, so clade V's + consensus matches the insert worse than clade W's (``variolaB`` alone).""" + rng = random.Random(11) + + def mutate(seq: str, frac: float) -> str: + chars = list(seq) + for i in range(len(chars)): + if rng.random() < frac: + chars[i] = rng.choice("ACGT") + return "".join(chars) + + base = "".join(rng.choice("ACGT") for _ in range(6000)) + cowpox, variola = mutate(base, 0.03), mutate(base, 0.08) + variola_b = mutate(variola, 0.01) + query = list(cowpox) + query[2000:4000] = list(variola[2000:4000]) + msa = write_fasta(tmp_path / "panel.fasta", { + "query": "".join(query), "cowpox": cowpox, "variolaA": variola, + "variolaB": variola_b, "junk1": mutate(base, 0.15), "junk2": mutate(base, 0.15), + }) + lineages = {"cowpox": "CPX", "variolaA": "V", "junk1": "V", "junk2": "V", "variolaB": "W"} + return msa, lineages + + +def test_method_comparison_names_the_reattributed_donor(tmp_path: Path, logger) -> None: + """Re-attribution re-labelled the region but not the per-method breakdown, so + `recombination_methods.tsv` and the report's method table kept the old donor.""" + msa, lineages = _reattribution_panel(tmp_path) + out = tmp_path / "out" + run_recomb( + RecombParams(msa=msa, output=out, query="query", plot_format="png", + lineage_map=lineages, reattribute_donors=True, cluster_lineages=False), + logger, + ) + regions = list(csv.DictReader((out / "recombination_regions.tsv").open(), delimiter="\t")) + methods = list(csv.DictReader((out / "recombination_methods.tsv").open(), delimiter="\t")) + assert [r["minor_parent"] for r in regions] == ["variolaB"] # re-attributed from variolaA + assert [m["minor_parent"] for m in methods] == ["variolaB"] + assert [(m["query_start"], m["query_end"]) for m in methods] == [ + (r["query_start"], r["query_end"]) for r in regions + ] diff --git a/tests/unit/test_report_plots.py b/tests/unit/test_report_plots.py new file mode 100644 index 0000000..d50412c --- /dev/null +++ b/tests/unit/test_report_plots.py @@ -0,0 +1,91 @@ +"""The plots label what the regions table says, in the same colours.""" + +from __future__ import annotations + +from pathlib import Path + +from matplotlib.figure import Figure + +from tessera.recomb.regions import Region +from tessera.recomb.report import pair_datasets +from tessera.recomb.report_plots import ( + GREY, + _color_map, + _palette, + _shade_regions, + build_interactive_figure, + region_labels, +) +from tessera.recomb.similarity import compute_similarity + +from ..conftest import recombinant_msa + + +def _region(minor: str, major: str = "A", *, start: int = 2000, end: int = 4000, + donor_absent: bool = False) -> Region: + return Region( + minor_parent=minor, major_parent=major, msa_start=start, msa_end=end, + query_start=start, query_end=end, n_windows=10, + mean_sim_minor=0.99, mean_sim_major=0.95, margin=0.04, donor_absent=donor_absent, + ) + + +# --- colours cover every label a region names ------------------------------ + +def test_region_labels_appends_region_parents_to_the_plotted_datasets() -> None: + labels = region_labels(["A"], [_region("B"), _region("other", donor_absent=True)]) + assert labels == ["A", "B"] # the donor is added; a donor-absent stand-in is not + assert region_labels(["A", "B"], [_region("B")]) == ["A", "B"] # no duplicates + + +def test_donor_outside_the_top_n_keeps_a_colour(tmp_path: Path) -> None: + """With --top-n 1 only the backbone was in the colour map, so the donor band and + its label were drawn grey.""" + result = compute_similarity(str(recombinant_msa(tmp_path, recombinant=True)), "query") + fig = build_interactive_figure(result, ["A"], [_region("B")]) + donor_colour = _color_map(["A", "B"])["B"] + band = next(a for a in fig.layout.annotations if a.text == "B") + assert band.font.color == donor_colour + assert donor_colour != GREY + + +# --- a donor-absent region is not a recombinant from the backbone ---------- + +def test_static_plot_does_not_call_a_donor_absent_region_recombinant() -> None: + """The stand-in label of a donor-absent region is the closest reference, often the + backbone itself: the legend read "recombinant: " in the backbone's colour.""" + ax = Figure().subplots() # no pyplot state, no display backend needed + _shade_regions(ax, [_region("A", donor_absent=True)], _palette(["A", "B"])) + labels = ax.get_legend_handles_labels()[1] + assert labels == ["donor absent (no close reference)"] + + +def test_static_plot_still_labels_a_called_donor() -> None: + ax = Figure().subplots() + _shade_regions(ax, [_region("B")], _palette(["A", "B"])) + labels = ax.get_legend_handles_labels()[1] + assert labels == ["recombinant: B"] + + +def test_interactive_plot_marks_a_donor_absent_region(tmp_path: Path) -> None: + result = compute_similarity(str(recombinant_msa(tmp_path, recombinant=True)), "query") + fig = build_interactive_figure(result, ["A", "B"], [_region("A", donor_absent=True)]) + texts = [a.text for a in fig.layout.annotations] + assert texts == ["donor absent"] + assert fig.layout.annotations[0].font.color == GREY + + +# --- the pair plot is major versus leading minor --------------------------- + +def test_pair_plot_shows_the_backbone_and_the_leading_donor() -> None: + """It showed the top two window winners, which are the backbone and a + near-duplicate of it when the panel holds one.""" + regions = [_region("variola", "cowpox", start=2000, end=2600), + _region("camelpox", "cowpox", start=4000, end=4100)] + assert pair_datasets(regions, ["cowpox", "cowpox2"]) == ["cowpox", "variola"] + + +def test_pair_plot_falls_back_to_the_window_ranking_without_a_called_donor() -> None: + assert pair_datasets([], ["cowpox", "cowpox2"]) == ["cowpox", "cowpox2"] + absent = [_region("cowpox", "cowpox", donor_absent=True)] + assert pair_datasets(absent, ["cowpox", "cowpox2"]) == ["cowpox", "cowpox2"] diff --git a/tests/unit/test_report_signal.py b/tests/unit/test_report_signal.py new file mode 100644 index 0000000..dd690db --- /dev/null +++ b/tests/unit/test_report_signal.py @@ -0,0 +1,67 @@ +"""The parent-free signal section states the alpha the run used, and says "not +testable" when the PHI test could not have rejected.""" + +from __future__ import annotations + +import json +import logging + +from tessera.recomb.diagnostics import RecombinationSignal +from tessera.recomb.report_html import _signal_html +from tessera.recomb.report_text import write_profile_tsv +from tessera.recomb.run import RecombParams, run_recomb + + +def _signal(phi_p: float | None, n_informative: int = 500) -> RecombinationSignal: + return RecombinationSignal( + n_informative=n_informative, phi_p=phi_p, phi_observed=0.1, phi_window=100, + rmin=2, rmin_intervals=[(100, 200), (900, 1000)], profile=[(150, 150, 0.2)], + ) + + +def test_signal_section_says_not_testable() -> None: + html = _signal_html(_signal(None, n_informative=54), alpha=0.05) + assert "not testable" in html + assert "p =" not in html + assert "no significant recombination signal" not in html + assert "54 informative sites" in html + # 54 sites and a window of 53 ranks is also untestable (every pair is inside it), so + # "do not exceed the window" would be false there; describe the actual condition. + assert "do not exceed" not in html + assert "every pair" in html + assert "--phi-window" in html # tells the reader what to change + + +def test_profile_tsv_header_marks_an_untestable_phi(tmp_path) -> None: + write_profile_tsv(_signal(None, n_informative=54), tmp_path, logging.getLogger("tessera")) + header = (tmp_path / "recombination_profile.tsv").read_text().splitlines()[0] + assert header.split("\t")[1] == "NA" + assert "not testable" in header + + +def test_run_reports_an_untestable_phi(example_data, tmp_path, logger) -> None: + # cryptic_insert has 54 informative sites: fewer than the default window of 100. + out = tmp_path / "cryptic" + run_recomb( + RecombParams(msa=example_data / "cryptic_insert.msa.fasta", output=out, + query="query", plot_format="png"), + logger, + ) + run = json.loads((out / "run_provenance.json").read_text())["run"] + assert run["recombination signal (PHI)"].startswith("not testable") + assert "not testable" in (out / "report.html").read_text() + + +def test_report_states_the_alpha_the_run_used(example_data, tmp_path, logger) -> None: + """The section always printed "alpha 0.05" and judged significance at 0.05, while + the per-region PHI flag used --alpha.""" + out = tmp_path / "divergent" + run_recomb( + RecombParams(msa=example_data / "divergent.msa.fasta", output=out, query="query", + window_size=300, window_step=30, plot_format="png", alpha=0.2, + phi_window=20), + logger, + ) + report = (out / "report.html").read_text() + assert "(alpha 0.2;" in report + assert "(alpha 0.05;" not in report diff --git a/tests/unit/test_report_stale_text.py b/tests/unit/test_report_stale_text.py new file mode 100644 index 0000000..8ce2e07 --- /dev/null +++ b/tests/unit/test_report_stale_text.py @@ -0,0 +1,78 @@ +"""Report and help text that describes what the run did, not an earlier version of it.""" + +from __future__ import annotations + +import re +from pathlib import Path + +from typer.testing import CliRunner + +from tessera.cli.main import app +from tessera.recomb.report_assets import _REFERENCES +from tessera.recomb.run import RecombParams, run_recomb + + +def _footer(out: Path) -> str: + report = (out / "report.html").read_text() + return re.search(r"
.*?
", report, flags=re.S).group(0) + + +def _run(example_data: Path, out: Path, logger, **kwargs) -> Path: + run_recomb( + RecombParams(msa=example_data / "divergent.msa.fasta", output=out, query="query", + window_size=300, window_step=30, **kwargs), + logger, + ) + return out + + +def test_footer_lists_the_plots_in_the_format_they_were_written( + example_data, tmp_path, logger +) -> None: + out = tmp_path / "out" + out.mkdir() + (out / "similarity_pair.pdf").write_text("left over from an earlier run") + footer = _footer(_run(example_data, out, logger, plot_format="png")) + assert "similarity_pair.png" in footer + assert "similarity_top3.png" in footer + # Was hard-coded; and a file this run did not write is not listed even if present. + assert ".pdf" not in footer + + +def test_footer_lists_only_files_that_exist(example_data, tmp_path, logger) -> None: + out = _run(example_data, tmp_path / "out", logger, plot_format="png", + methods=("hmm",), phi=False) + footer = _footer(out) + listed = re.findall(r"(.*?)", footer) + assert listed, "the footer should list the companion files" + for name in listed: + assert (out / name).exists(), f"{name} is listed but was not written" + # A single-caller run writes no method comparison; --no-phi writes no profile. + assert "recombination_methods.tsv" not in listed + assert "recombination_profile.tsv" not in listed + assert "run_provenance.json" in listed + + +def test_methods_paragraph_describes_the_default_ensemble( + example_data, tmp_path, logger +) -> None: + report = (_run(example_data, tmp_path / "out", logger, plot_format="png") + / "report.html").read_text() + methods = re.search(r'
.*?
', report, flags=re.S).group(0) + assert "(HMM and the 3SEQ triplet test)" not in methods # the two-caller description + for caller in ("3SEQ", "MaxChi", "Bootscan"): + assert caller in methods + + +def test_references_cite_every_default_caller() -> None: + titles = " ".join(title for title, _ in _REFERENCES) + assert "MaxChi" in titles + assert "Bootscan" in titles + + +def test_method_help_lists_every_caller_and_the_real_default() -> None: + result = CliRunner().invoke(app, ["recomb", "--help"], env={"COLUMNS": "200"}) + assert result.exit_code == 0 + text = " ".join(result.output.replace("│", " ").split()) + assert "geneconv" in text # run by `all`, and was missing from the list + assert "all but the legacy heuristic" not in text # geneconv and barcode are opt-in too diff --git a/tests/unit/test_report_summary.py b/tests/unit/test_report_summary.py new file mode 100644 index 0000000..62cbfd9 --- /dev/null +++ b/tests/unit/test_report_summary.py @@ -0,0 +1,68 @@ +"""The report's headline numbers agree with the regions they summarise.""" + +from __future__ import annotations + +import numpy as np + +from tessera.recomb.regions import Region +from tessera.recomb.report_html import _summary, _verdict_html + + +class _Result: + """The one attribute ``_summary`` reads: a 6000-base, gap-free query.""" + + query = "query" + query_cumulative = np.arange(6001) + + +def _region(minor: str, start: int, end: int, *, donor_absent: bool = False) -> Region: + return Region( + minor_parent=minor, major_parent="backbone", msa_start=start, msa_end=end, + query_start=start, query_end=end, n_windows=5, + mean_sim_minor=0.99, mean_sim_major=0.92, margin=0.07, + qvalue=1e-9, support=0.95, methods=("hmm", "3seq"), donor_absent=donor_absent, + ) + + +def test_overlapping_regions_are_counted_once() -> None: + """The ensemble keeps overlapping regions that name different donors separate. + Adding their lengths reported a 2.2 kb union as 3.3 kb (54.6 % of the query).""" + regions = [_region("donorB", 1900, 3100), _region("donorA", 2025, 4100)] + s = _summary(_Result(), regions, ["backbone", "donorA", "donorB"]) + assert s["recomb_bp"] == 2200 # the union 1900-4100, not 1200 + 2075 + assert round(s["pct"], 1) == 36.7 + assert s["n_regions"] == 2 # both regions are still listed + + +def test_recombinant_fraction_never_exceeds_the_query() -> None: + regions = [_region(f"donor{i}", 0, 6000) for i in range(3)] + s = _summary(_Result(), regions, ["backbone"]) + assert s["recomb_bp"] == 6000 + assert s["pct"] == 100.0 + + +def test_touching_and_nested_regions() -> None: + touching = [_region("donorA", 1000, 2000), _region("donorB", 2000, 3000)] + assert _summary(_Result(), touching, ["backbone"])["recomb_bp"] == 2000 + nested = [_region("donorA", 1000, 4000), _region("donorB", 2000, 2500)] + assert _summary(_Result(), nested, ["backbone"])["recomb_bp"] == 3000 + + +def test_disjoint_regions_still_add_up() -> None: + regions = [_region("donorA", 500, 1500), _region("donorB", 3000, 3500)] + assert _summary(_Result(), regions, ["backbone"])["recomb_bp"] == 1500 + + +def test_donor_absent_regions_stay_out_of_the_recombinant_span() -> None: + regions = [_region("donorA", 500, 1500), _region("x", 3000, 5000, donor_absent=True)] + s = _summary(_Result(), regions, ["backbone"]) + assert s["recomb_bp"] == 1000 + assert s["n_absent"] == 1 + + +def test_verdict_states_the_union() -> None: + regions = [_region("donorB", 1900, 3100), _region("donorA", 2025, 4100)] + s = _summary(_Result(), regions, ["backbone", "donorA", "donorB"]) + verdict = _verdict_html(s, "query", {}) + assert "2.2 kb" in verdict + assert "36.7%" in verdict diff --git a/validation/README.md b/validation/README.md index f052270..3996976 100644 --- a/validation/README.md +++ b/validation/README.md @@ -61,10 +61,19 @@ python validation/run_specificity.py --reps 3 # quick look python validation/run_specificity.py --min-methods 1 # without the agreement gate ``` -It is deliberately sensitive to the failure mode it was built for. With the default -agreement gate the scan is clean; dropping the gate with `--min-methods 1` surfaces the -single-caller regions again (measured at 4 replicates: 9/16 runs, 10 false regions, -almost all from one caller). Treat a non-zero total as a regression to explain. +It is deliberately sensitive to the failure mode it was built for. **The harness's own +default is `--min-methods 2`; the `tessera` CLI default is `--min-methods 1`**, so a clean +default harness run does not describe the shipped default. Measured at 3 replicates on +2026-10-01 (commit `457bdfb`): + +| gate | runs with a false region | false regions | source | +|---|---|---|---| +| `--min-methods 2` (harness default) | 0/12 (CI 0-24 %) | 0 | -- | +| `--min-methods 1` (CLI default) | 7/12 (58 %, CI 32-81 %) | 8 | hmm = 8 | + +The positive control was detected 3/3 with the correct donor at both gates (median +breakpoint error 55 bp). Treat a non-zero total at `--min-methods 2` as a regression to +explain, and quote the `--min-methods 1` row when describing what a default run reports. **Caveat.** JC69 on a fixed topology is simpler than real viral evolution, and modest replicate counts carry real sampling error -- hence the intervals. These numbers @@ -105,11 +114,13 @@ have not been fetched, so a partial setup still reports cleanly. | `enterovirus_e11` | enterovirus ~7.3 kb | Echovirus-11 x Coxsackievirus-B1, breakpoint in P2 | mafft | recombination detected; both parents named (checks parents-present, not the ambiguous backbone direction) | | `hiv_crf02ag` | HIV-1 ~9.2 kb | CRF02_AG (IbNG): A backbone + subtype-G segments | mafft | major A; G donor region(s) over the pol/vif and vpu-env inserts | | `hcv_2k1b` | HCV ~9.4 kb | RF1_2k/1b: genotype-2k 5' + 1b 3', breakpoint in NS2/NS3 | mafft | major 1b; genotype-2 donor over the 5'; **precise breakpoint ~nt 3187 recovered** | -| `hcv_clonal_1b` | HCV ~9.4 kb | pure genotype-1b (non-recombinant control) | mafft | resolves to 1b throughout; **0 regions** (real-data specificity) | +| `hcv_clonal_1b` | HCV ~9.4 kb | pure genotype-1b (non-recombinant control) | mafft | resolves to 1b throughout, but **currently FAILS**: one 12 bp MaxChi-only region (GT2a donor, q = 0.045) is reported where none is expected (real-data specificity) | -Each reproduces its published event (or, for the clonal control, its *absence* of -recombination) end-to-end; the current run is **7 PASS, 0 FAIL** -(`orthopox_example` SKIPs until its 7-genome collection is built). `hcv_2k1b` is the +The six recombinant datasets that ran reproduce their published events end-to-end +(`orthopox_example` SKIPs until its 7-genome collection is built). The clonal control does +not currently pass: the run on 2026-10-01 was **6 PASS, 1 FAIL** (`hcv_clonal_1b`), **1 +SKIP**. The failing region rests on a single caller at the CLI default `--min-methods 1`; +it is recorded here as an open specificity item rather than hidden. `hcv_2k1b` is the first real **precise-breakpoint** check and `hcv_clonal_1b` the first real **false-positive** (specificity) check. Accessions are listed per dataset in `datasets.json` (`provenance` field) and were confirmed