diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 666a65a..de5a510 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -16,7 +16,7 @@ jobs: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: '3.12' + python-version: '3.11' - name: Install dependencies run: | python -m pip install --upgrade pip diff --git a/.gitignore b/.gitignore index d548845..9f32e87 100755 --- a/.gitignore +++ b/.gitignore @@ -119,3 +119,6 @@ vartracker/.snakemake/ !vartracker/test_data/ !vartracker/test_data/** vartracker/test_data/**/.DS_Store + +# Claude Code local config +.claude/ diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..6938a42 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,69 @@ +# Changelog + +All notable changes to vartracker are documented in this file, in addition to +(not replacing) the project's versioned GitHub releases. + +Format loosely follows [Keep a Changelog](https://keepachangelog.com/). + +## [2.3.0] - 2026-07-31 + +### Added + +- `--local-csq` option (`vartracker vcf`/`bam`/`end-to-end`) to switch `bcftools csq` to + independent, per-variant consequence calling instead of the default joint/compound calling - + useful for gene-dense, high-variant-density data where unphased, sub-consensus variants can + cluster in the same gene without genotype evidence they co-occur. +- `--out`/`--outdir` for `vartracker plot heatmap`, to control the output filename and location. +- `--multiallelic-overflow` (`error`/`drop-lowest-af`/`skip-site`) to control how sites with more + than two surviving ALT alleles in a single sample are handled before `bcftools csq`. +- LoFreq primer-overlap rescue for amplicon data (`--primer-bed`, `--lofreq-primer-rescue`, + `--lofreq-rescue-*` thresholds), with raw, rescued, and filtered-out LoFreq calls written to + per-sample audit TSVs. +- Consensus genome generation in the `bam`/`end-to-end` workflows (`--consensus-snp-min-af`, + `--consensus-snp-thresh`, `--consensus-indel-thresh`), writing both a simple and an + IUPAC-aware consensus FASTA per sample. +- Gene-name disambiguation across contigs/replicons: `ambiguous_gene_names()` and updates to + `gene_lengths_from_gff3`/`generate_gene_table` so genes that legitimately reuse a name on + different contigs (e.g. `repA` on a chromosome and a plasmid) are tracked and reported separately. +- Gene-wise summary plot (`mutations_per_gene.pdf`) is now capped to the top genes by newly + emerged variants on large (e.g. bacterial) genomes, via `--max-plot-genes`/`--plot-genes`; the + tabular output is unaffected. +- New standalone plotting subcommands (`plot heatmap`, `plot genome`, `plot trajectory`, + `plot turnover`, `plot lifespan`) and heatmap filtering options (`--aa-exclude`, `--aa-include`, + `--only-persistent`, `--only-new`, `--gene-include`/`--gene-exclude`, `--variant-type`, `--qc`, + `--min-prop-passing-qc`, `--min-persistence`, `--min-max-af`, `--min-sample-af`, + `--sample-subset`, `--hide-singletons`, `--min-depth`, `--x-labels`, `--title`, + `--literature-csv`). +- Structured logging for Snakemake pipeline steps. +- Documentation: a guide for building a custom `--literature-csv` for non-SARS-CoV-2 pathogens, a + QC-column interpretation guide, and a bacterial-genome validation/limitations section. + +### Changed + +- `vartracker plot heatmap` CLI flags were de-prefixed for standalone use (e.g. + `--heatmap-aa-exclude` -> `--aa-exclude`); the legacy prefixed forms are kept as aliases. +- The default heatmap now always shows a variant's canonical row, even if that row is + joint/compound, instead of unconditionally hiding all joint rows; `--include-joint` reveals + any additional joint/compound annotation-group rows a variant has. +- Per-sample variant QC is now reported directly (`per_sample_variant_qc`, `P`/`F`) rather than + only as an overall pass/fail, and the default heatmap visually flags failing cells. + +### Fixed + +- Duplicate `results.csv` rows for the same variant (one per `bcftools csq` annotation group) were + previously plotted as independent data points, inflating turnover counts; plotting now collapses + to one row per (variant, sample), taking the union of presence and the maximum allele frequency. +- `joint_variant` could miss rows with a genuinely compound `bcsq_nt_notation` when a variant had + several sibling rows at the same position; it now flags every such row. +- Heatmap label truncation could make two distinct variants render an identical label, silently + dropping one of them; colliding labels are now disambiguated with a positional suffix. +- Gene names that collide across contigs/replicons were previously merged into a single row in + gene-length and gene-wise summary tables; they are now tracked and reported per-contig. +- Percent-encoded gene names read from `bcftools`' `BCSQ` field (e.g. `ercS%27`) were not + URL-decoded on the variant-calling side, so they never matched the (correctly decoded) + annotation side and were silently dropped from gene-wise summaries. +- Corrected several low-frequency multiallelic-site handling bugs, so distinct ALT alleles at the + same position are kept separate through preprocessing and correctly rejoined before + `bcftools csq`. +- Fixed potential results-row duplication for some samples during heatmap generation. +- Fixed longitudinal variant-rescue edge cases and LoFreq strand-bias rescue thresholds. diff --git a/CITATION.cff b/CITATION.cff index 55d2fd2..1381fd3 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,8 +1,8 @@ cff-version: 1.2.0 message: "If you use this software, please cite it as below." title: "vartracker" -version: "2.2.1" -date-released: "2026-05-06" +version: "2.3.0" +date-released: "2026-07-31" license: MIT repository-code: "https://github.com/charlesfoster/vartracker" url: "https://github.com/charlesfoster/vartracker" @@ -22,7 +22,7 @@ preferred-citation: - family-names: Foster given-names: Charles title: "vartracker" - version: "2.2.1" + version: "2.3.0" doi: "10.5281/zenodo.18452274" url: "https://github.com/charlesfoster/vartracker" - date-released: "2026-05-06" + date-released: "2026-07-31" diff --git a/README.md b/README.md index 30a3e45..bf41af9 100755 --- a/README.md +++ b/README.md @@ -23,6 +23,7 @@ A bioinformatics pipeline to summarise variants called against a reference in a - [Quick Start](#quick-start) - [Output](#output) - [What does vartracker do?](#what-does-vartracker-do) +- [Limitations](#limitations) - [Citation](#citation) - [License](#license) - [Contributing](#contributing) @@ -167,7 +168,7 @@ Docker is a self-contained reproducible option. If you publish the image, record set it when running to include it in the run manifest: ```bash -export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.2.1 +export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.3.0 export VARTRACKER_CONTAINER_DIGEST=sha256:... ``` @@ -309,8 +310,11 @@ for both `.depth.txt` and `_depth.txt` patterns when preparing its internal test ### Mode-specific options - `vartracker vcf` – accepts core analysis options such as `--min-snv-freq`, `--min-indel-freq`, - `--allele-frequency-tag`, `--multiallelic-overflow`, `--name`, `--outdir`, `--sample-cap`, `--manifest-level`, and literature controls - (`--search-pokay`, `--literature-csv`). Use `--test` to run the bundled smoke test. + `--allele-frequency-tag`, `--multiallelic-overflow`, `--local-csq`, `--name`, `--outdir`, + `--sample-cap`, `--manifest-level`, and literature controls (`--search-pokay`, + `--literature-csv`). Use `--test` to run the bundled smoke test. + `--max-plot-genes` and `--plot-genes` control the gene-wise summary figure only (see + [Limitations](#limitations)); the tabular/TSV output always includes every annotated gene. - `vartracker bam` – everything from `vcf`, plus Snakemake options: `--snakemake-outdir`, `--cores`, `--snakemake-dryrun`, `--verbose`, `--redo`, `--rulegraph`, `--primer-bed`, `--lofreq-primer-rescue`, `--consensus-snp-min-af`, @@ -318,24 +322,38 @@ for both `.depth.txt` and `_depth.txt` patterns when preparing its internal test - `vartracker end-to-end` – similar to `bam`, with optional amplicon clipping controls: `--primer-bed` and `--ampliconclip-tolerance` (default: `1`). Supplying `--primer-bed` also enables LoFreq primer-overlap rescue by default. -- `vartracker plot heatmap` (`hm`) – regenerate the heatmap from an existing vartracker results CSV, including all heatmap customization filters. +- `vartracker plot heatmap` (`hm`) – regenerate the heatmap from an existing vartracker results CSV, including all heatmap customisation filters. - `vartracker plot genome` – plot SNP positions along the genome or a selected gene region using all observed allele-frequency values for each variant. - `vartracker plot trajectory` – plot allele-frequency trajectories for a selected or auto-ranked subset of variants, optionally in takeover mode using threshold lines and threshold-based filtering. - `vartracker plot turnover` – plot new-versus-lost longitudinal turnover from the filtered result set. - `vartracker plot lifespan` – plot first-to-last detection spans for a selected or auto-ranked subset of variants. +QC threshold note: +- `--min-snv-freq`, `--min-indel-freq`, and `--min-depth` are configurable allele-frequency and + read-depth thresholds applied when summarising and visualising longitudinal variant calls. Their + defaults reflect our own genomic surveillance and longitudinal sequencing workflows and should be + treated as starting points, not universally applicable QC recommendations. The appropriate + thresholds for a given study depend on its objective and on the sequencing protocol, depth, + variant caller, and empirically established error profile of the upstream workflow: use more + stringent thresholds when specificity is prioritised or the input data have higher error rates, + and only lower thresholds for low-frequency variant analysis when this is supported by a suitably + validated upstream workflow. + Consequence-calling note: - Vartracker keeps distinct ALT alleles at the same position separate during preprocessing, then rejoins them immediately before `bcftools csq` so codon-level consequences can still be inferred correctly. -- If more than two ALT alleles remain present in a single sample at one genomic position after frequency filtering, vartracker defaults to stopping with an informative error before `bcftools csq`. This is the safest behavior and the default `--multiallelic-overflow error` mode. +- If more than two ALT alleles remain present in a single sample at one genomic position after frequency filtering, vartracker defaults to stopping with an informative error before `bcftools csq`. This is the safest behaviour and the default `--multiallelic-overflow error` mode. - `--multiallelic-overflow drop-lowest-af` continues by removing the lowest-frequency retained ALT allele(s) for the affected sample before `bcftools csq`, and prints a warning describing the site and the dropped allele(s). - `--multiallelic-overflow skip-site` continues by skipping consequence calling for the affected site entirely, leaving those variants in the results as unannotated rows and printing a warning describing the site. Heatmap filtering: -- `vcf`, `bam`, and `end-to-end` always write the default heatmap. To customize heatmap content after a run, use `vartracker plot heatmap results.csv [options]`. -- By default, all consequence classes are included except joint variants. Use `--include-joint` to show joint variants. +- `vcf`, `bam`, and `end-to-end` always write the default heatmap. To customise heatmap content after a run, use `vartracker plot heatmap results.csv [options]`. +- By default, each variant is shown once, using its canonical row (whether that row is joint or + not - see [Limitations](#limitations)). Use `--include-joint` to additionally reveal extra + joint/compound annotation-group rows for variants that have more than one. - `--aa-exclude`: comma-separated `type_of_change` patterns to exclude. Wildcards are supported. - `--aa-include`: comma-separated `type_of_change` patterns to include. -- `--only-persistent`: only include `new_persistent` variants. +- `--only-persistent`: only include new variants present at the final timepoint (`new_persistent` or + `new_intermittent`; see [Persistence labels](#persistence-labels)). - `--only-new`: only include variants with `variant_status == new`. - `--gene-include` and `--gene-exclude`: comma-separated gene patterns. - `--variant-type`: comma-separated variant-type patterns such as `snp` or `indel`. @@ -350,13 +368,19 @@ Heatmap filtering: - `--x-labels sample-number`: label heatmap x-axis columns by `sample_number` instead of sample name. - `--title`: set the heatmap plot title. The default is `Variant allele frequencies`. - `--literature-csv`: include literature links in the interactive HTML heatmap using a literature hits CSV. +- `--out` (`vartracker plot heatmap` only): write the heatmap using this path as the base name, + e.g. `--out plots/myheatmap` writes `plots/myheatmap.pdf` and `plots/myheatmap.html`. +- `--outdir` (`vartracker plot heatmap` only): output directory for heatmap files (default: + beside `results.csv`). - Example: `--aa-exclude "synonymous,*frameshift*,stop_gained"` Standalone plot filtering: - `--gene`, `--effect`, `--min-af`, `--max-af`: restrict the plotted result set before ranking/selection. - `--variants` or `--variant-file`: explicitly choose variants and preserve that order. - `--sample-min`, `--sample-max`: restrict the passage/sample-number window. -- `--persistent-only` and `--new-only`: keep only persistent new variants or only variants with `variant_status == new`. +- `--persistent-only` and `--new-only`: keep only new variants present at the final timepoint + (`new_persistent` or `new_intermittent`; see [Persistence labels](#persistence-labels)) or only + variants with `variant_status == new`. - `trajectory` and `lifespan` auto-select a limited subset by default (`--top-n`) to stay readable. - `turnover` uses all filtered variants by default and is also written automatically during the main `vcf`/`bam`/`end-to-end` workflows as `variant_turnover_plot.pdf`. - `genome` uses SNPs only by default, keeps all observed allele-frequency values for each plotted variant, and writes `variant_genome_plot.pdf` during the main workflows. @@ -371,7 +395,7 @@ Genome plot options: - `--gene`: zoom to a single gene region. - `--aa-scale`: with `--gene`, use amino-acid coordinates on the x-axis. - `--cds-scale`: with `--gene`, use CDS-relative nucleotide coordinates on the x-axis. -- `--focus-coords`: highlight nucleotide or amino-acid coordinate ranges, depending on the current x-axis mode. Separate color groups with `;`, ranges within a group with `,`, and optionally prefix a group with `Name:`. +- `--focus-coords`: highlight nucleotide or amino-acid coordinate ranges, depending on the current x-axis mode. Separate colour groups with `;`, ranges within a group with `,`, and optionally prefix a group with `Name:`. - `--focus-region-file`: read named focus region groups from a `.json`, `.csv`, or `.tsv` file for an inset legend. - `--show-intersections`: add a compact `Region | Variant` table below the genome plot for highlighted-region hits. - In the genome plot, undetected samples are rendered at the detection threshold rather than zero; by default this floor is `0.03`, or `--min-af` if supplied, and the dashed guide line follows that same threshold. @@ -430,6 +454,41 @@ vartracker [mode] input_data.csv --literature-csv pokay_database.csv -o results/ Alternatively, pass `--search-pokay` to automatically download and search against the Pokay SARS-CoV-2 literature database. +#### Building a custom literature database for other pathogens + +`--search-pokay` only covers SARS-CoV-2. For other pathogens, supply your own CSV via +`--literature-csv`. Variant lookup during vartracker analysis is based on a CSV-format file that +is either generated automatically (`--search-pokay`) or supplied by the user (`--literature-csv +`). The expected structure of the file is described using the `vartracker schema literature` +command. In brief, after deriving appropriate information from the scientific literature, users +can create their own lookup table by creating a new CSV file whereby each row corresponds to a +variant of interest, with `gene` and `mutation` required and `category`, `information`, and +`reference` recommended: + +- **`gene`**: must exactly match (case-sensitive) the gene/product name assigned to that variant + by `bcftools csq` using the GFF3/GenBank annotation supplied via `--gff3`. For non-SARS-CoV-2 + pathogens this is simply the gene name as it appears in your annotation file — vartracker's + SARS-CoV-2-specific remapping of `ORF1ab` into individual `nsp1`–`nsp16` names does not apply + outside SARS-CoV-2, so for other pathogens use the gene names exactly as they appear in your + GFF3. +- **`mutation`**: the amino acid consequence in short-hand notation *without* a gene prefix (e.g. + `D614G`, not `S:D614G`). Each row describes a single mutation; if you have information on + several mutations in the same gene, add one row per mutation. Note that matching is done via + substring containment on this column, so avoid overly short or ambiguous notations that could + unintentionally match unrelated variants (e.g. a bare position number). +- **`category`**: a free-text label used to group/colour variants in output tables and the + heatmap. There's no fixed vocabulary — choose categories meaningful for your pathogen (e.g. + "resistance", "immune_escape", "homoplasy"). +- **`information`**: free-text description of the mutation's putative effect, drawn from the + literature. +- **`reference`**: one or more supporting DOIs or URLs, semicolon-delimited if there are multiple. + +A minimal template with this exact structure is provided at +`test_data/mock_literature/mock_literature.csv`; the SARS-CoV-2-specific `pokay_database.csv` +generated by `--search-pokay` follows the same schema and can also be used as a real-world +formatting reference, bearing in mind its `ORF1ab`/`nsp` gene naming is SARS-CoV-2-specific and +shouldn't be copied for other pathogens. + ### Command Line Reference ``` @@ -516,7 +575,8 @@ vartracker produces several output files: - **`_variants.rescued.tsv`** (`bam`/`end-to-end`): LoFreq primer-overlap rescue audit table, empty when rescue is disabled or no variants are rescued - **`_variants.filtered_out.tsv`** (`bam`/`end-to-end`): Raw LoFreq calls excluded from the final VCF, including filter reason and call metrics - **new_mutations.csv**: Mutations not present in the first sample -- **persistent_new_mutations.csv**: New mutations that persist to the final sample +- **persistent_new_mutations.csv**: New mutations present at the final sample (`new_persistent` or + `new_intermittent`; see [Persistence labels](#persistence-labels)) - **cumulative_mutations.pdf**: Plot showing mutation accumulation over time - **mutations_per_gene.pdf**: Gene-wise mutation statistics - **variant_allele_frequency_heatmap.html**: Interactive heatmap with optional literature annotations @@ -527,6 +587,64 @@ vartracker produces several output files: By default the manifest is lightweight. Use `--manifest-level deep` to checksum all referenced input files (FASTQ/BAM/VCF/coverage) and include file sizes. +### Persistence labels + +The `persistence_status` column classifies each variant from `variant_status` +(`original`: present in the first sample; `new`: absent in the first sample) plus its presence +pattern across the rest of the samples: + +- `original_retained`: an `original` variant continuously present through the final sample. +- `original_intermittent`: an `original` variant present in the final sample, but absent from at + least one sample in between (i.e. lost and regained). +- `original_lost`: an `original` variant absent by the final sample. +- `new_persistent`: a `new` variant continuously present from its first appearance through the + final sample. +- `new_intermittent`: a `new` variant present in the final sample, but absent from at least one + sample between its first appearance and the final sample (i.e. it appeared, disappeared in a + later sample, then reappeared). +- `new_transient`: a `new` variant absent by the final sample. + +These labels are driven by presence/absence, not allele frequency, and depend only on the first, +last, and intervening samples - they say nothing on their own about whether an intervening absence +reflects genuine loss or a QC dropout (see the `per_sample_variant_qc` column in +[Output schema](#output-schema)). `--only-persistent` / `--persistent-only` filters (heatmap and +standalone plots) and `persistent_new_mutations.csv` include both `new_persistent` and +`new_intermittent` variants, since both reached the final timepoint; the label only distinguishes +the path taken to get there. + +### Interpreting the QC columns + +`results.csv` records presence/absence per sample (`presence_absence`, `Y`/`N`), but an `N` does +not always mean the variant was confidently confirmed absent. At low sequencing depth, a variant +can go undetected simply because there was insufficient coverage to call it either way - this is +indistinguishable, from the VCF alone, from genuine absence. The QC columns exist to flag this: + +- `per_sample_variant_qc`: a per-sample `P`/`F` flag. `F` means that sample had no + variant-supporting read *and* site coverage below `--min-depth` (default: 10) - i.e. absence + could not be confidently distinguished from dropout/non-detection at that sample. `P` means the + call (presence or absence) was made with confidence. +- `all_samples_pass_qc`: `true` only if every sample is `P`. +- `proportion_samples_passing_qc`: the fraction of samples that are `P`. + +**Practical guidance:** if `all_samples_pass_qc` is `false` for a variant, inspect +`per_sample_variant_qc` to see exactly which sample(s) it failed at - e.g. `P / P / F / P / P / P` +identifies the third sample as the QC failure. Before treating an `N` in `presence_absence` as +evidence a variant was truly lost or never present, check the corresponding position in +`per_sample_variant_qc`: an `N` paired with `F` should be read as "not detected", not "confirmed +absent" - this is especially relevant for low-frequency variants near the allele-frequency or depth +thresholds (`--min-snv-freq`, `--min-indel-freq`, `--min-depth`), where dropout is more likely than +at high-confidence, high-depth sites. This ambiguity also propagates into `persistence_status` (see +[Persistence labels](#persistence-labels)): an apparent loss-then-reappearance (`*_intermittent`) +may reflect genuine intermittent presence, or simply a low-coverage sample in between. + +**QC in the heatmap.** The default heatmap marks `F` cells visually rather than just via colour: the +static PDF draws an unfilled black-bordered rectangle over any cell whose sample failed QC for that +variant; the interactive HTML version uses a dark inset ring plus a hover tooltip reading +`QC=FAIL`. To exclude variants that don't pass QC from a plot entirely (rather than just flagging +the cells), use `--qc` and `--min-prop-passing-qc` (see +[Mode-specific options](#mode-specific-options)), or inspect `per_sample_variant_qc` directly for +the samples of interest. + ### Output schema The results table schema is documented in `docs/OUTPUT_SCHEMA.md`. You can also print it from the CLI: @@ -552,7 +670,7 @@ vartracker schema literature The pipeline performs the following analysis: -1. **VCF Standardization**: Normalizes and standardizes input VCF files, preserving distinct ALT alleles at the same genomic position +1. **VCF Standardisation**: Normalises and standardises input VCF files, preserving distinct ALT alleles at the same genomic position 2. **Variant Merging**: Combines all longitudinal samples 3. **Annotation**: Adds amino acid consequences using `bcftools csq` on the merged VCF so sample-specific joint consequences are inferred from each sample's surviving ALT combination 4. **Comprehensive Analysis**: For each variant, determines: @@ -566,12 +684,83 @@ The pipeline performs the following analysis: 5. **Visualization**: Generates plots for mutation accumulation and gene-wise statistics 6. **Functional Annotation**: (optional) Searches against literature databases for known functional impacts +## Limitations + +vartracker was designed for viral pathogens with small, compact genomes (SARS-CoV-2: ~30 kb, +12 genes). The underlying analysis - VCF standardisation, merging, annotation, and +original/new/persistent/transient classification - scales to larger genomes without modification. +The practical constraint on larger genomes (e.g. bacterial pathogens, which can carry thousands of +annotated genes) is **visualisation**, not computation: + +- The gene-wise summary figure (`mutations_per_gene.pdf`) plots one bar per gene per panel. On a + genome with thousands of annotated genes this becomes unreadable as a static image regardless of how many variants + are actually present, because the plot iterates over every annotated gene, not just genes that + carry a variant. +- By default, the figure is capped to the top 30 genes, ranked by number of newly emerged variants + (ties broken by total variant count), via `--max-plot-genes`. Use `--plot-genes` to instead name + an explicit set of genes to plot. **This cap applies to the figure only** - the tabular/TSV output + always contains every annotated gene, so no data is discarded by this option. +- When the figure is truncated, this is stated directly on the figure itself (e.g. "top 30 of 412 + genes with variants"); if nothing was truncated, no such note is shown. + +### Bacterial genomes + +vartracker works well and efficiently at bacterial genome scale. It has been validated in `vcf` +mode (i.e. from pre-called VCFs and coverage files, not the `bam`/`end-to-end` read-mapping +workflow) against simulated *Pseudomonas aeruginosa* PAO1 data (NC_002516.2, 6.26 Mb, 5,573 CDS +features) across two scenarios - 80 and 1,000 simulated variants, each across 6 timepoints. The +smaller, 80-variant scenario completed in approximately 11 seconds of wall-clock time with +approximately 0.9 GiB peak memory; the larger, 1,000-variant scenario completed in approximately +22 seconds with approximately 2.8 GiB peak memory. At this scale, the practical caveat is not +runtime or memory but the **interpretability of joint/compound amino-acid consequences in +gene-dense hotspots**, discussed below. + +**Joint vs local `bcftools csq` calling.** By default, vartracker calls consequences jointly (the +`bcftools csq` default), so that variants close enough together to plausibly affect the same +codon(s) are described together as a single, compound amino-acid change. This is the correct +behaviour for genuinely linked variants, but on gene-dense, high-variant-density data - common in +bacterial within-host or experimental-evolution datasets, and rare in vartracker's original viral +use case - many unphased, sub-consensus variants can cluster in the same gene without genotype +evidence that they actually co-occur on the same haplotype. Joint calling then produces long, +compound descriptions that are technically correct but hard to read, and can fragment a single +variant's presence/absence trajectory across samples. The `--local-csq` option (see +`vartracker --help`) switches to independent, SnpEff-like per-variant consequence calling, at the +cost of no longer detecting genuinely combined effects between physically linked variants. Whichever +mode is used, rows describing a joint/compound consequence are flagged in the `joint_variant` +column of `results.csv`, which can be used to identify or filter these rows after the fact. This +column is now fully reliable: on the PAO1 validation dataset, 100% of genuinely compound +`bcftools csq` rows are correctly flagged. + +**Heatmap `--include-joint` semantics.** By default, the heatmap shows one row per variant - its +canonical row, whether that row happens to be joint or not. `--include-joint` additionally reveals +extra joint/compound annotation-group rows for variants that have more than one. This is a different +kind of control from the gene-wise figure's `--max-plot-genes` cap: there is no row cap on the +heatmap. + +**Heatmap legibility at bacterial scale.** Unlike the gene-wise figure, the heatmap has no built-in +row cap, and is effectively illegible as a static image. However, you can still open it and scroll to read the rows. Alternatively, for large numbers of variants, narrow the heatmap using the +"Heatmap filtering" options described under [Mode-specific options](#mode-specific-options) - for +example `--gene-include`, `--hide-singletons`, `--min-max-af`, and `--only-persistent` - or generate +multiple heatmaps over subsets of genes/samples rather than relying on a single, unfiltered plot. + +**Coverage-file disk and memory footprint.** Disk and memory usage for coverage/depth files scale +with genome length multiplied by timepoint count. For reference, the PAO1 validation used 6 depth +files at approximately 141 MB each (846 MB total) for one 6.3 Mb genome across 6 timepoints. Users +planning many-timepoint experimental-evolution designs (often dozens of timepoints) on genomes +larger than PAO1 should budget disk and memory accordingly. + +Separately, the bundled `pokay` functional-annotation database +(see [Using Literature Database](#using-literature-database)) is specific to SARS-CoV-2 mutations +and is not applied to, or meaningful for, other pathogens. A custom literature CSV following the +same schema can be supplied via `--literature-csv` for other organisms; see +[Building a custom literature database for other pathogens](#building-a-custom-literature-database-for-other-pathogens). + ## Citation When using vartracker, please cite the software release you used. Citation metadata is provided in `CITATION.cff`, and GitHub releases are archived on Zenodo. -- Foster, C. (2026). *vartracker* (Version 2.2.1). Zenodo. https://doi.org/10.5281/zenodo.18452274 +- Foster, C. (2026). *vartracker* (Version 2.3.0). Zenodo. https://doi.org/10.5281/zenodo.18452274 Note: the DOI above is the Zenodo concept DOI for all versions; a version-specific DOI is minted by Zenodo after each GitHub release. diff --git a/docs/OUTPUT_SCHEMA.md b/docs/OUTPUT_SCHEMA.md index 578c9d1..ff824f1 100644 --- a/docs/OUTPUT_SCHEMA.md +++ b/docs/OUTPUT_SCHEMA.md @@ -24,13 +24,13 @@ Columns that encode per-sample values are slash-separated and ordered by the inp | type_of_variant | string | Variant type derived from the VCF entry. | | snp, indel | | type_of_change | string | Functional change classification from bcftools csq. | | synonymous, missense, frameshift, ... | | variant_status | string | Whether the variant is present in the first sample. | | original, new | -| persistence_status | string | Persistence class based on first and last sample presence. | | original_retained, original_lost, new_persistent, new_transient | +| persistence_status | string | Persistence class based on presence at the first and last sample, and whether presence was continuous in between. The '_intermittent' classes mean the variant was present at both the first/last (original) or reappeared by the last sample (new), but was absent from at least one sample in between. | | original_retained, original_intermittent, original_lost, new_persistent, new_intermittent, new_transient | | presence_absence | string (slash-separated) | Per-sample presence (Y) or absence (N), ordered by input. | | Y/N | | first_appearance | string | Sample name where the variant first appears. | | | | last_appearance | string | Sample name where the variant last appears. | | | -| all_samples_pass_qc | boolean | True if every sample passes per-sample variant QC. | | true, false | +| all_samples_pass_qc | boolean | True only if every sample is 'P' in per_sample_variant_qc. False means at least one sample's presence/absence call could not be confidently distinguished from dropout/non-detection - inspect per_sample_variant_qc to see which sample(s) are affected before treating this variant's presence/absence pattern as reliable. | | true, false | | proportion_samples_passing_qc | number | Proportion of samples passing per-sample variant QC. | fraction | 0-1 | -| per_sample_variant_qc | string (slash-separated) | Per-sample QC flags (P/F) ordered by input. | | P, F | +| per_sample_variant_qc | string (slash-separated) | Per-sample QC flags (P/F) ordered by input. A sample fails ('F') when there is no variant-supporting read and site coverage is below --min-depth, i.e. when genuine absence of the variant cannot be distinguished from dropout/non-detection; 'P' means absence (or presence) was confidently called at that sample. E.g. 'P / P / F / P / P / P' identifies the third sample as the one where QC failed. The default heatmap marks 'F' cells visually: an unfilled black-bordered rectangle in the static PDF, and a dark inset ring plus a 'QC=FAIL' hover tooltip in the interactive HTML version. | | P, F | | aa1_total_properties | string | Physicochemical properties for the reference amino acid. | | semicolon-separated properties | | aa2_total_properties | string | Physicochemical properties for the alternate amino acid. | | semicolon-separated properties | | aa1_unique_properties | string | Properties unique to the reference amino acid. | | semicolon-separated properties | diff --git a/pyproject.toml b/pyproject.toml index 790e3cc..99696a7 100755 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "vartracker" -version = "2.2.1" +version = "2.3.0" authors = [ {name = "Dr Charles Foster"}, ] diff --git a/scripts/validation/pao1/README.md b/scripts/validation/pao1/README.md new file mode 100644 index 0000000..ef3415e --- /dev/null +++ b/scripts/validation/pao1/README.md @@ -0,0 +1,73 @@ +# PAO1 bacterial-genome validation + +Scripts to reproduce the *Pseudomonas aeruginosa* PAO1 (NC_002516.2) bacterial-scale validation +referenced in the main README's [Bacterial genomes](../../../README.md#bacterial-genomes) section +and in the accompanying manuscript. This regenerates two synthetic longitudinal datasets from +scratch, then runs `vartracker vcf` on each: + +- **scenario_a**: 80 variants, uniform gene weighting, 6 timepoints. +- **scenario_b**: 1,000 variants, Zipf-biased gene weighting (concentrating variants in a smaller + set of genes, to stress-test gene-dense hotspots), 6 timepoints. + +## What's here + +- `simulate_pao1.py`: builds both scenarios against a supplied PAO1 reference FASTA/GFF3. For each + scenario it assigns simulated variants a `persistence_status` pattern (`original_retained`, + `original_lost`, `new_persistent`, `new_transient`, `original_intermittent`, + `new_intermittent`), synthesises per-timepoint VCFs and matching depth/coverage files, and writes + a `_ground_truth.csv` (intended pattern, presence/absence, allele frequency per + timepoint) alongside the `_input.csv` manifest `vartracker vcf` expects. All randomness + is seeded (`scenario_a`: seed 42; `scenario_b`: seed 43; coverage: seed 123), so re-running + produces identical output. +- `run_validation.sh`: end-to-end driver - builds the reference bundle, runs the simulation, then + runs `vartracker vcf` on both scenarios (benchmarked with GNU time if available). + +## Requirements + +- `vartracker` installed and on `PATH` (see the main [README](../../../README.md#installation)). +- `bcftools` and `bgzip` (external dependencies of `simulate_pao1.py` and of vartracker itself). +- `python3` with `numpy` and `pandas`. +- Optional, to reproduce reported timing/peak-memory figures: GNU time - `gtime` on macOS + (`brew install gnu-time`) or `/usr/bin/time -v` on Linux. Without it the scenarios still run, just + unbenchmarked. + +## Usage + +```bash +scripts/validation/pao1/run_validation.sh [outdir] # outdir defaults to ./pao1_validation +``` + +This is equivalent to running the steps individually: + +```bash +vartracker prepare reference --accessions NC_002516.2 \ + --outdir pao1_validation/refs/pao1 --prefix pao1_ref + +python3 scripts/validation/pao1/simulate_pao1.py \ + --fasta pao1_validation/refs/pao1/pao1_ref.fa \ + --gff3 pao1_validation/refs/pao1/pao1_ref.gff3 \ + --base-outdir pao1_validation + +vartracker vcf pao1_validation/scenario_a/scenario_a_input.csv \ + --reference pao1_validation/refs/pao1/pao1_ref.fa \ + --gff3 pao1_validation/refs/pao1/pao1_ref.gff3 \ + --outdir pao1_validation/results/pao1_80_variants + +vartracker vcf pao1_validation/scenario_b/scenario_b_input.csv \ + --reference pao1_validation/refs/pao1/pao1_ref.fa \ + --gff3 pao1_validation/refs/pao1/pao1_ref.gff3 \ + --outdir pao1_validation/results/pao1_1000_variants +``` + +## Checking the results + +- Compare `results/pao1_*_variants/results.csv` (specifically `variant`, `variant_status`, and + `persistence_status`) against the corresponding `scenario_{a,b}_ground_truth.csv` to confirm every + simulated variant and its intended temporal pattern was recovered. +- Timing and peak memory are hardware-dependent, so exact figures will vary by machine; the README + and manuscript report approximately 11 s / 0.9 GiB peak memory for scenario_a and approximately + 22 s / 2.8 GiB peak memory for scenario_b, measured via `gtime -v` on the reference machine used + for validation. +- Neither the generated coverage files (large - tens to hundreds of MB per timepoint) nor the + reference bundle are committed to this repository; both are deterministically regenerated by the + steps above. diff --git a/scripts/validation/pao1/run_validation.sh b/scripts/validation/pao1/run_validation.sh new file mode 100644 index 0000000..c649fda --- /dev/null +++ b/scripts/validation/pao1/run_validation.sh @@ -0,0 +1,61 @@ +#!/usr/bin/env bash +set -euo pipefail + +# Reproduces the PAO1 bacterial-genome validation described in the vartracker +# README ("Bacterial genomes" section under Limitations) and manuscript: a +# synthetic 80-variant scenario and a synthetic 1000-variant scenario, each +# spanning 6 longitudinal timepoints, simulated against Pseudomonas aeruginosa +# PAO1 (NC_002516.2) and analysed with `vartracker vcf`. +# +# Usage: +# scripts/validation/pao1/run_validation.sh [outdir] +# +# Requires on PATH: vartracker, bcftools, bgzip, python3 (with numpy and +# pandas). Optional, to reproduce the reported timing/peak-memory figures: +# GNU time - `gtime` on macOS (brew install gnu-time) or `/usr/bin/time -v` +# on Linux. Without either, the analysis still runs, just unbenchmarked. + +OUTDIR="${1:-pao1_validation}" +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" + +mkdir -p "$OUTDIR" + +echo "== Step 1: build the PAO1 reference bundle (NC_002516.2) ==" +vartracker prepare reference \ + --accessions NC_002516.2 \ + --outdir "$OUTDIR/refs/pao1" \ + --prefix pao1_ref + +echo "== Step 2: simulate scenario_a (80 variants, uniform gene weighting) and scenario_b (1000 variants, Zipf-biased gene weighting) ==" +python3 "$SCRIPT_DIR/simulate_pao1.py" \ + --fasta "$OUTDIR/refs/pao1/pao1_ref.fa" \ + --gff3 "$OUTDIR/refs/pao1/pao1_ref.gff3" \ + --base-outdir "$OUTDIR" + +if command -v gtime >/dev/null 2>&1; then + TIME_CMD="gtime -v" +elif /usr/bin/time -v true >/dev/null 2>&1; then + TIME_CMD="/usr/bin/time -v" +else + echo "No GNU time found (gtime / time -v); continuing without wall-clock/memory benchmarking." >&2 + TIME_CMD="" +fi + +echo "== Step 3: run vartracker vcf on scenario_a ==" +$TIME_CMD vartracker vcf "$OUTDIR/scenario_a/scenario_a_input.csv" \ + --reference "$OUTDIR/refs/pao1/pao1_ref.fa" \ + --gff3 "$OUTDIR/refs/pao1/pao1_ref.gff3" \ + --outdir "$OUTDIR/results/pao1_80_variants" + +echo "== Step 4: run vartracker vcf on scenario_b ==" +$TIME_CMD vartracker vcf "$OUTDIR/scenario_b/scenario_b_input.csv" \ + --reference "$OUTDIR/refs/pao1/pao1_ref.fa" \ + --gff3 "$OUTDIR/refs/pao1/pao1_ref.gff3" \ + --outdir "$OUTDIR/results/pao1_1000_variants" + +echo +echo "Done." +echo "Cross-check $OUTDIR/results/pao1_*_variants/results.csv against" +echo "$OUTDIR/scenario_{a,b}/scenario_{a,b}_ground_truth.csv, and (if GNU time" +echo "ran above) compare Elapsed / Maximum resident set size against the" +echo "figures quoted in the README's Bacterial genomes section." diff --git a/scripts/validation/pao1/simulate_pao1.py b/scripts/validation/pao1/simulate_pao1.py new file mode 100644 index 0000000..54c8c27 --- /dev/null +++ b/scripts/validation/pao1/simulate_pao1.py @@ -0,0 +1,445 @@ +#!/usr/bin/env python3 +"""Simulate longitudinal VCF + coverage data against the PAO1 reference bundle +for vartracker bacterial-scale validation (Scenario A ~80 variants, Scenario B ~1000). +""" + +import argparse +import subprocess +import sys +from pathlib import Path + +import numpy as np +import pandas as pd + +N_TIMEPOINTS = 6 +BASES = ["A", "C", "G", "T"] + +PATTERNS = [ + "original_retained", + "original_lost", + "new_persistent", + "new_transient", + "original_intermittent", + "new_intermittent", +] +PATTERN_WEIGHTS = [0.30, 0.15, 0.25, 0.15, 0.08, 0.07] +MIN_INTERMITTENT = 3 + + +def parse_fasta(path): + header = None + seq_chunks = [] + with open(path) as fh: + for line in fh: + line = line.rstrip("\n") + if line.startswith(">"): + header = line[1:].split()[0] + else: + seq_chunks.append(line) + return header, "".join(seq_chunks).upper() + + +def parse_cds(gff3_path, contig): + cds = [] + with open(gff3_path) as fh: + for line in fh: + if line.startswith("#"): + continue + parts = line.rstrip("\n").split("\t") + if len(parts) < 9 or parts[0] != contig or parts[2] != "CDS": + continue + start, end, strand, attrs = int(parts[3]), int(parts[4]), parts[6], parts[8] + gene = None + for kv in attrs.split(";"): + if kv.startswith("gene="): + gene = kv.split("=", 1)[1] + break + if gene is None or end - start < 20: + continue + cds.append({"gene": gene, "start": start, "end": end, "strand": strand}) + return cds + + +def zipf_weights(n, rng, exponent=1.1): + order = rng.permutation(n) + ranks = np.empty(n, dtype=float) + ranks[order] = np.arange(1, n + 1) + w = 1.0 / (ranks**exponent) + return w / w.sum() + + +def choose_positions( + cds_list, n_variants, ref_seq, rng, indel_fraction=0.08, biased=False +): + """Pick n_variants distinct (pos, ref, alt, gene, var_type) tuples from CDS regions.""" + genes = np.array([c["gene"] for c in cds_list]) + unique_genes, gene_index = np.unique(genes, return_inverse=True) + # map each unique gene name -> list of cds_list indices (handles duplicate-name genes e.g. speA) + gene_to_cds_idx = {g: [] for g in unique_genes} + for i, g in enumerate(genes): + gene_to_cds_idx[g].append(i) + + if biased: + weights = zipf_weights(len(unique_genes), rng) + else: + weights = np.full(len(unique_genes), 1.0 / len(unique_genes)) + + used_positions = set() + variants = [] + attempts = 0 + max_attempts = n_variants * 50 + + while len(variants) < n_variants and attempts < max_attempts: + attempts += 1 + gene_choice = rng.choice(unique_genes, p=weights) + cds_idx = rng.choice(gene_to_cds_idx[gene_choice]) + cds = cds_list[cds_idx] + is_indel = rng.random() < indel_fraction + span = rng.integers(1, 3) if is_indel else 0 + lo, hi = cds["start"] + 3, cds["end"] - 3 - span + if hi <= lo: + continue + pos = int(rng.integers(lo, hi)) + span_positions = set(range(pos, pos + span + 2)) + if span_positions & used_positions: + continue + + ref_base = ref_seq[pos - 1] + if ref_base not in BASES: + continue + + if not is_indel: + alt_base = rng.choice([b for b in BASES if b != ref_base]) + ref, alt, var_type = ref_base, alt_base, "snp" + else: + if rng.random() < 0.5: + # insertion + inserted = "".join(rng.choice(BASES, size=span)) + ref, alt, var_type = ref_base, ref_base + inserted, "indel" + else: + # deletion: anchor base + (span) deleted ref bases + del_seq = ref_seq[pos - 1 : pos - 1 + span + 1] + if len(del_seq) < span + 1 or any(b not in BASES for b in del_seq): + continue + ref, alt, var_type = del_seq, del_seq[0], "indel" + + used_positions |= span_positions + variants.append( + { + "chrom": None, # filled by caller + "pos": pos, + "ref": ref, + "alt": alt, + "gene": gene_choice, + "var_type": var_type, + } + ) + + if len(variants) < n_variants: + print( + f"WARNING: only placed {len(variants)}/{n_variants} variants after {attempts} attempts", + file=sys.stderr, + ) + return variants + + +def gen_presence_af(pattern, var_type, rng): + af_floor = 0.12 if var_type == "indel" else 0.05 + af_ceiling = 0.95 + present = [False] * N_TIMEPOINTS + + if pattern == "original_retained": + present = [True] * N_TIMEPOINTS + base = rng.uniform(0.2, 0.5) + af_vals = [ + np.clip(base + rng.uniform(-0.05, 0.05), af_floor, af_ceiling) + for _ in range(N_TIMEPOINTS) + ] + + elif pattern == "original_lost": + last_present = int(rng.integers(1, 5)) # 1..4 inclusive -> absent by tp5 + present = [i <= last_present for i in range(N_TIMEPOINTS)] + n_present = last_present + 1 + af_vals = list( + np.clip( + np.linspace(0.55, 0.15, n_present) + + rng.uniform(-0.03, 0.03, n_present), + af_floor, + af_ceiling, + ) + ) + + elif pattern == "new_persistent": + first_present = int(rng.integers(1, 5)) # 1..4 + present = [i >= first_present for i in range(N_TIMEPOINTS)] + n_present = N_TIMEPOINTS - first_present + af_vals = list( + np.clip( + np.linspace(0.12, 0.85, n_present) + + rng.uniform(-0.03, 0.03, n_present), + af_floor, + af_ceiling, + ) + ) + + elif pattern == "new_transient": + first_present = int(rng.integers(1, 4)) # 1..3 + last_present = int( + rng.integers(first_present, min(first_present + 3, 5)) + ) # < 5 + present = [first_present <= i <= last_present for i in range(N_TIMEPOINTS)] + n_present = last_present - first_present + 1 + base = rng.uniform(0.15, 0.4) + af_vals = list( + np.clip( + [base + rng.uniform(-0.05, 0.05) for _ in range(n_present)], + af_floor, + af_ceiling, + ) + ) + + elif pattern == "original_intermittent": + gap_start = int(rng.integers(1, 4)) # 1..3 + gap_len = int( + rng.integers(1, 5 - gap_start) + ) # keep gap inside 1..4, tp5 stays present + present = [True] * N_TIMEPOINTS + for i in range(gap_start, min(gap_start + gap_len, 5)): + present[i] = False + present[0] = True + present[5] = True + n_present = sum(present) + base = rng.uniform(0.2, 0.45) + af_vals = list( + np.clip( + [base + rng.uniform(-0.05, 0.05) for _ in range(n_present)], + af_floor, + af_ceiling, + ) + ) + + elif pattern == "new_intermittent": + first_present = int(rng.integers(1, 3)) # 1..2 + gap_start = int(rng.integers(first_present + 1, 5)) # inside 1..4 + gap_len = int(rng.integers(1, max(2, 5 - gap_start))) + present = [False] * N_TIMEPOINTS + for i in range(first_present, N_TIMEPOINTS): + present[i] = True + for i in range(gap_start, min(gap_start + gap_len, 5)): + present[i] = False + present[0] = False + present[5] = True + n_present = sum(present) + base = rng.uniform(0.15, 0.4) + af_vals = list( + np.clip( + [base + rng.uniform(-0.05, 0.05) for _ in range(n_present)], + af_floor, + af_ceiling, + ) + ) + + else: + raise ValueError(pattern) + + # sanity checks on invariants + if pattern.startswith("original"): + assert present[0] is True + else: + assert present[0] is False + if pattern in ( + "original_retained", + "new_persistent", + "original_intermittent", + "new_intermittent", + ): + assert present[5] is True + else: + assert present[5] is False + + af_iter = iter(af_vals) + af_full = [round(float(next(af_iter)), 4) if p else None for p in present] + dp_full = [int(rng.integers(80, 200)) if p else None for p in present] + return present, af_full, dp_full + + +def assign_patterns(n, rng): + patterns = list(rng.choice(PATTERNS, size=n, p=PATTERN_WEIGHTS)) + for target in ("original_intermittent", "new_intermittent"): + count = patterns.count(target) + if count < MIN_INTERMITTENT: + donor_pool = [ + i + for i, p in enumerate(patterns) + if p not in ("original_intermittent", "new_intermittent") + ] + n_needed = MIN_INTERMITTENT - count + idx_to_convert = rng.choice(donor_pool, size=n_needed, replace=False) + for idx in idx_to_convert: + patterns[idx] = target + return patterns + + +def simulate_scenario( + name, n_variants, cds_list, contig, ref_seq, seed, biased, outdir: Path +): + rng = np.random.default_rng(seed) + variants = choose_positions(cds_list, n_variants, ref_seq, rng, biased=biased) + for v in variants: + v["chrom"] = contig + + patterns = assign_patterns(len(variants), rng) + + records = [] + for v, pattern in zip(variants, patterns): + present, af, dp = gen_presence_af(pattern, v["var_type"], rng) + variant_status = "original" if pattern.startswith("original") else "new" + records.append( + { + **v, + "pattern": pattern, + "variant_status": variant_status, + "present": present, + "af": af, + "dp": dp, + } + ) + + records.sort(key=lambda r: r["pos"]) + + # sanity: no duplicate positions + positions = [r["pos"] for r in records] + assert len(positions) == len(set(positions)), "duplicate positions generated" + + outdir.mkdir(parents=True, exist_ok=True) + vcf_dir = outdir / "vcfs" + vcf_dir.mkdir(exist_ok=True) + + for tp in range(N_TIMEPOINTS): + rows = [r for r in records if r["present"][tp]] + rows.sort(key=lambda r: r["pos"]) + raw_path = vcf_dir / f"pao1_tp{tp}.raw.vcf" + with open(raw_path, "w") as fh: + fh.write("##fileformat=VCFv4.2\n") + fh.write(f"##contig=\n") + fh.write('##INFO=\n') + fh.write( + '##INFO=\n' + ) + fh.write("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n") + for r in rows: + af_val = r["af"][tp] + dp_val = r["dp"][tp] + qual = 60 + fh.write( + f"{r['chrom']}\t{r['pos']}\t.\t{r['ref']}\t{r['alt']}\t{qual}\tPASS\tDP={dp_val};AF={af_val}\n" + ) + + final_gz = vcf_dir / f"pao1_tp{tp}.vcf.gz" + sort_cmd = f"bcftools sort -Ov {raw_path} 2>/dev/null | bgzip > {final_gz}" + subprocess.run(sort_cmd, shell=True, check=True) + subprocess.run(["bcftools", "index", "-f", str(final_gz)], check=True) + raw_path.unlink() + + # ground truth table for cross-checking step 4 + gt_rows = [] + for r in records: + gt_rows.append( + { + "chrom": r["chrom"], + "pos": r["pos"], + "ref": r["ref"], + "alt": r["alt"], + "gene": r["gene"], + "var_type": r["var_type"], + "variant_status": r["variant_status"], + "intended_pattern": r["pattern"], + "presence": "/".join("Y" if p else "N" for p in r["present"]), + "af": "/".join("" if a is None else str(a) for a in r["af"]), + } + ) + gt_df = pd.DataFrame(gt_rows).sort_values(["gene", "pos"]) + gt_df.to_csv(outdir / f"{name}_ground_truth.csv", index=False) + + # input CSV for vartracker vcf mode + csv_rows = [] + for tp in range(N_TIMEPOINTS): + csv_rows.append( + { + "sample_name": f"pao1_tp{tp}", + "sample_number": tp, + "reads1": "", + "reads2": "", + "bam": "", + "vcf": str((vcf_dir / f"pao1_tp{tp}.vcf.gz").resolve()), + "coverage": str( + (outdir.parent / "coverage" / f"pao1_tp{tp}_depth.txt").resolve() + ), + } + ) + pd.DataFrame(csv_rows).to_csv(outdir / f"{name}_input.csv", index=False) + + # summary counts + pattern_counts = pd.Series([r["pattern"] for r in records]).value_counts() + gene_counts = pd.Series([r["gene"] for r in records]).value_counts() + print( + f"[{name}] {len(records)} distinct variants placed across {gene_counts.shape[0]} genes" + ) + print(f"[{name}] pattern mix:\n{pattern_counts}") + print(f"[{name}] top genes by variant count:\n{gene_counts.head(10)}") + return records + + +def generate_coverage(contig, genome_length, coverage_dir: Path, seed=123): + coverage_dir.mkdir(parents=True, exist_ok=True) + pos = np.arange(1, genome_length + 1, dtype=np.int64) + for tp in range(N_TIMEPOINTS): + out_path = coverage_dir / f"pao1_tp{tp}_depth.txt" + if out_path.exists(): + continue + rng = np.random.default_rng(seed + tp) + depth = rng.integers(80, 151, size=genome_length, dtype=np.int32) + df = pd.DataFrame({"chrom": contig, "pos": pos, "depth": depth}) + df.to_csv(out_path, sep="\t", header=False, index=False) + print(f"wrote coverage file {out_path} ({genome_length} lines)") + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--fasta", required=True) + ap.add_argument("--gff3", required=True) + ap.add_argument("--base-outdir", required=True) + args = ap.parse_args() + + base = Path(args.base_outdir) + contig, ref_seq = parse_fasta(args.fasta) + print(f"Loaded reference contig {contig}, length {len(ref_seq)}") + cds_list = parse_cds(args.gff3, contig) + print(f"Parsed {len(cds_list)} CDS features on {contig}") + + generate_coverage(contig, len(ref_seq), base / "coverage") + + simulate_scenario( + "scenario_a", + 80, + cds_list, + contig, + ref_seq, + seed=42, + biased=False, + outdir=base / "scenario_a", + ) + simulate_scenario( + "scenario_b", + 1000, + cds_list, + contig, + ref_seq, + seed=43, + biased=True, + outdir=base / "scenario_b", + ) + + +if __name__ == "__main__": + main() diff --git a/tests/test_analysis.py b/tests/test_analysis.py index a97c007..3030c57 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -2,18 +2,329 @@ from __future__ import annotations +import os +import string + +import matplotlib.pyplot as plt import pandas as pd import pytest +import seaborn as sns from vartracker.analysis import ( + _build_variant_key, _heatmap_figure_size, search_literature, _prepare_variant_heatmap_matrix, process_joint_variants, generate_variant_heatmap, + generate_gene_table, + plot_gene_table, + parse_plot_genes_arg, + select_canonical_row_positions, + select_genes_for_plot, ) +def _gene_table_rows(gene_table, gene): + rows = gene_table[gene_table["gene"] == gene] + return dict(zip(rows["type"], rows["number"])) + + +def test_generate_gene_table_keeps_single_contig_genes_unqualified(): + """Baseline: single-contig data (e.g. viral references) is unaffected.""" + table = pd.DataFrame( + { + "gene": ["S", "S", "N"], + "chrom": ["NC_045512.2", "NC_045512.2", "NC_045512.2"], + "type_of_change": ["missense", "missense", "synonymous"], + "presence_absence": ["NY", "NY", "YY"], + } + ) + + gene_table = generate_gene_table(table, gene_lengths={"S": 3822, "N": 1260}) + + assert set(gene_table["gene"]) == {"S", "N"} + assert _gene_table_rows(gene_table, "S")["total"] == 2 + + +def test_generate_gene_table_splits_colliding_gene_names_across_contigs(): + """Same gene name on chromosome + plasmid must not be summed together.""" + table = pd.DataFrame( + { + "gene": ["repA", "repA", "dnaA"], + "chrom": ["chrom1", "plasmid1", "chrom1"], + "type_of_change": ["missense", "missense", "missense"], + "presence_absence": ["NY", "NY", "NY"], + } + ) + + gene_table = generate_gene_table( + table, + gene_lengths={ + "dnaA": 300, + "repA (chrom1)": 201, + "repA (plasmid1)": 150, + }, + ) + + genes_present = set(gene_table["gene"]) + assert "repA (chrom1)" in genes_present + assert "repA (plasmid1)" in genes_present + assert "repA" not in genes_present + + # Each contig's copy of repA should retain its own count, not a merged total. + assert _gene_table_rows(gene_table, "repA (chrom1)")["total"] == 1 + assert _gene_table_rows(gene_table, "repA (plasmid1)")["total"] == 1 + + +def test_generate_gene_table_ambiguous_genes_param_prevents_data_loss(): + """A gene flagged ambiguous by the annotation, but with variants on only + one contig in this particular table, must still surface under its + qualified label so it matches the gene-length scaffold and isn't dropped + by the scaffold merge.""" + table = pd.DataFrame( + { + "gene": ["repA"], + "chrom": ["chrom1"], + "type_of_change": ["missense"], + "presence_absence": ["NY"], + } + ) + + gene_table = generate_gene_table( + table, + gene_lengths={ + "repA (chrom1)": 201, + "repA (plasmid1)": 150, + }, + ambiguous_genes={"repA"}, + ) + + genes_present = set(gene_table["gene"]) + assert "repA (chrom1)" in genes_present + assert "repA" not in genes_present + assert _gene_table_rows(gene_table, "repA (chrom1)")["total"] == 1 + # The scaffold still surfaces the other contig's copy with a zero count. + assert _gene_table_rows(gene_table, "repA (plasmid1)")["total"] == 0 + + +def _make_variant_row(gene, type_of_change, presence_absence, chrom="chrom1"): + return { + "gene": gene, + "chrom": chrom, + "type_of_change": type_of_change, + "presence_absence": presence_absence, + } + + +def test_select_genes_for_plot_keeps_everything_when_annotation_under_cap(): + """When the whole annotated gene set already fits under the cap (e.g. a + viral reference), nothing is filtered - this is what keeps SARS-CoV-2 + figures byte-identical regardless of the new default cap.""" + table = pd.DataFrame( + [ + _make_variant_row("S", "missense", "NY"), + _make_variant_row("N", "synonymous", "YY"), + ] + ) + gene_table = generate_gene_table( + table, gene_lengths={"S": 3822, "N": 1260, "E": 228} + ) + + plot_table, subtitle = select_genes_for_plot(gene_table, max_plot_genes=30) + + assert subtitle is None + pd.testing.assert_frame_equal( + plot_table.reset_index(drop=True), gene_table.reset_index(drop=True) + ) + + +def test_select_genes_for_plot_ranks_by_new_mutations_with_total_tiebreak(): + """Ranking must use newly emerged variants first, falling back to total + variants only to break ties - not the other way around.""" + table = pd.DataFrame( + [ + # Gene A: 3 new mutations, 5 total. + _make_variant_row("A", "missense", "NY"), + _make_variant_row("A", "missense", "NY"), + _make_variant_row("A", "missense", "NY"), + _make_variant_row("A", "synonymous", "YY"), + _make_variant_row("A", "synonymous", "YY"), + # Gene B: 3 new mutations (tied with A), only 2 total. + _make_variant_row("B", "missense", "NY"), + _make_variant_row("B", "missense", "NY"), + _make_variant_row("B", "missense", "NY"), + # Gene C: fewer new mutations than A/B, but more total variants. + _make_variant_row("C", "missense", "NY"), + *[_make_variant_row("C", "synonymous", "YY") for _ in range(9)], + # Gene D: no new mutations at all. + _make_variant_row("D", "synonymous", "YY"), + ] + ) + gene_table = generate_gene_table( + table, gene_lengths={"A": 100, "B": 100, "C": 100, "D": 100, "E": 100} + ) + + plot_table, subtitle = select_genes_for_plot(gene_table, max_plot_genes=2) + + assert subtitle == "top 2 of 4 genes with variants" + assert set(plot_table["gene"]) == {"A", "B"} + + +def test_select_genes_for_plot_omits_subtitle_when_variant_genes_fit_cap(): + """Even if the annotated genome is large, if the genes that actually + carry variants already fit under the cap, nothing was truncated and the + subtitle must be omitted (e.g. never print "top 30 of 12").""" + table = pd.DataFrame( + [ + _make_variant_row("A", "missense", "NY"), + _make_variant_row("B", "missense", "NY"), + ] + ) + gene_lengths = {"A": 100, "B": 100} + gene_lengths.update({f"unused{i}": 100 for i in range(20)}) + gene_table = generate_gene_table(table, gene_lengths=gene_lengths) + + plot_table, subtitle = select_genes_for_plot(gene_table, max_plot_genes=10) + + assert subtitle is None + assert set(plot_table["gene"]) == {"A", "B"} + + +def test_select_genes_for_plot_cap_zero_or_negative_raises(): + table = pd.DataFrame([_make_variant_row("A", "missense", "NY")]) + gene_table = generate_gene_table(table, gene_lengths={"A": 100}) + + with pytest.raises(ValueError): + select_genes_for_plot(gene_table, max_plot_genes=0) + + with pytest.raises(ValueError): + select_genes_for_plot(gene_table, max_plot_genes=-5) + + +def test_select_genes_for_plot_explicit_list_warns_and_drops_invalid(capsys): + table = pd.DataFrame( + [ + _make_variant_row("X", "missense", "NY"), + _make_variant_row("Y", "missense", "NY"), + ] + ) + gene_lengths = {"X": 100, "Y": 100, "Z": 100} # Z is annotated but has no variants + gene_table = generate_gene_table(table, gene_lengths=gene_lengths) + + plot_table, subtitle = select_genes_for_plot( + gene_table, plot_genes=["X", "Z", "FAKE"] + ) + + assert subtitle is None + assert set(plot_table["gene"]) == {"X"} + captured = capsys.readouterr() + assert "Z" in captured.out and "no variants" in captured.out + assert "FAKE" in captured.out and "reference annotation" in captured.out + + +def test_select_genes_for_plot_explicit_list_all_invalid_raises(): + table = pd.DataFrame([_make_variant_row("X", "missense", "NY")]) + gene_table = generate_gene_table(table, gene_lengths={"X": 100}) + + with pytest.raises(ValueError): + select_genes_for_plot(gene_table, plot_genes=["FAKE1", "FAKE2"]) + + +def test_select_genes_for_plot_explicit_list_overrides_max_plot_genes(): + table = pd.DataFrame( + [ + _make_variant_row("A", "missense", "NY"), + _make_variant_row("B", "missense", "NY"), + _make_variant_row("C", "missense", "NY"), + ] + ) + gene_table = generate_gene_table(table, gene_lengths={"A": 100, "B": 100, "C": 100}) + + plot_table, subtitle = select_genes_for_plot( + gene_table, max_plot_genes=1, plot_genes=["A", "B"] + ) + + assert subtitle is None + assert set(plot_table["gene"]) == {"A", "B"} + + +def test_parse_plot_genes_arg_comma_separated_deduplicates_and_preserves_order(): + genes = parse_plot_genes_arg("geneA, geneB,geneA, geneC") + assert genes == ["geneA", "geneB", "geneC"] + + +def test_parse_plot_genes_arg_reads_from_file(tmp_path): + gene_file = tmp_path / "genes.txt" + gene_file.write_text("geneA\n\n# a comment\ngeneB\ngeneA\n") + + genes = parse_plot_genes_arg(str(gene_file)) + + assert genes == ["geneA", "geneB"] + + +def _golden_plot_gene_table(gene_table, pname, outdir): + """Frozen copy of `plot_gene_table` as it existed before the + `--max-plot-genes`/`--plot-genes` options were added. Used only as a + ground truth for the byte-identical regression test below.""" + g = sns.catplot( + x="gene", + y="number", + col="type", + col_wrap=3, + data=gene_table, + kind="bar", + height=4, + aspect=1.2, + ) + g.set_axis_labels("", "Number of Mutations") + g.fig.subplots_adjust(top=0.9) + g.fig.suptitle(f"{pname}", weight="bold") + + for ax in g.axes.flat: + for label in ax.get_xticklabels(): + label.set_rotation(90) + + for ax, title in zip(g.fig.axes, list(gene_table["type"].unique())): + if pd.isna(title): + title_str = "None" + else: + title_str = str(title) if not isinstance(title, str) else title + ax.set_title(string.capwords(title_str.replace("_", " "))) + + plt.savefig( + os.path.join(outdir, "mutations_per_gene.pdf"), dpi=300, bbox_inches="tight" + ) + plt.close() + + +def test_plot_gene_table_default_matches_pre_cap_baseline(tmp_path, monkeypatch): + """With the default cap (30) and a SARS-CoV-2-sized (12 gene) reference, + output must be byte-identical to the pre-existing plotting behaviour, so + published manuscript figures don't change.""" + monkeypatch.setenv("SOURCE_DATE_EPOCH", "0") + + gene_names = [f"gene{i}" for i in range(12)] + rows = [] + for i, gene in enumerate(gene_names[:10]): + rows.append(_make_variant_row(gene, "missense", "NY" if i % 2 == 0 else "YY")) + table = pd.DataFrame(rows) + gene_lengths = {gene: 1000 for gene in gene_names} + gene_table = generate_gene_table(table, gene_lengths=gene_lengths) + + golden_dir = tmp_path / "golden" + new_dir = tmp_path / "new" + golden_dir.mkdir() + new_dir.mkdir() + + _golden_plot_gene_table(gene_table.copy(deep=True), "Baseline", str(golden_dir)) + plot_gene_table(gene_table.copy(deep=True), "Baseline", str(new_dir)) + + golden_bytes = (golden_dir / "mutations_per_gene.pdf").read_bytes() + new_bytes = (new_dir / "mutations_per_gene.pdf").read_bytes() + assert golden_bytes == new_bytes + + def test_search_literature_handles_nullable_boolean_masks(tmp_path): table = pd.DataFrame( { @@ -291,6 +602,74 @@ def test_prepare_variant_heatmap_matrix_excludes_wildcard_consequence_types(): assert list(matrix.index) == ["S:N501Y\n(A23063T)"] +def test_prepare_variant_heatmap_matrix_disambiguates_colliding_base_labels(): + """Two genuinely different variants that render to the same (truncated) + base label - a real bcftools-csq artefact on joint-called bacterial data + - must not collapse into a single heatmap row. A third, non-colliding + variant must come out byte-identical to plain `_resolve_variant_labels` + output, i.e. no unwanted disambiguation suffix.""" + long_aa = "CLLLDEFGHIJKLMNOPQRPA" # 21 chars: long enough to truncate + table = pd.DataFrame( + [ + { + "gene": "PA3025", + "amino_acid_consequence": long_aa, + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "alt_freq": "0.6", + "samples": "P0", + "variant": "C100T", + "start": 100, + }, + { + "gene": "PA3025", + "amino_acid_consequence": long_aa, + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "alt_freq": "0.7", + "samples": "P0", + "variant": "C200A", + "start": 200, + }, + { + "gene": "PA1000", + "amino_acid_consequence": "D50G", + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "alt_freq": "0.8", + "samples": "P0", + "variant": "A300G", + "start": 300, + }, + ] + ) + + matrix = _prepare_variant_heatmap_matrix(table, ["P0"], 0.0, 0.0) + + head = long_aa[:10] + tail = long_aa[-10:] + middle = len(long_aa) - 20 + truncated = f"{head}+{middle}{tail}" if middle > 0 else long_aa + + labels = list(matrix.index) + colliding_labels = [label for label in labels if f"PA3025:{truncated}" in label] + + # Both distinct variants survive (previously the second was silently + # dropped by `drop_duplicates(subset=["base_label"])`), and their + # displayed labels differ from one another. + assert len(colliding_labels) == 2 + assert colliding_labels[0] != colliding_labels[1] + assert any("@100" in label for label in colliding_labels) + assert any("@200" in label for label in colliding_labels) + + # The non-colliding variant is untouched: identical to what + # `_resolve_variant_labels` alone would produce. + assert "PA1000:D50G\n(A300G)" in labels + + def test_prepare_variant_heatmap_matrix_applies_extended_filters(): table = pd.DataFrame( [ @@ -368,6 +747,126 @@ def test_prepare_variant_heatmap_matrix_applies_extended_filters(): assert list(matrix.index) == ["S:D215G\n(A22206G)"] +def test_prepare_variant_heatmap_matrix_only_persistent_includes_new_intermittent(): + """--only-persistent must keep new_intermittent alongside new_persistent, + since both reached the final timepoint - only the path there differs.""" + table = pd.DataFrame( + [ + { + "gene": "S", + "amino_acid_consequence": "D215G", + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "persistence_status": "new_persistent", + "alt_freq": "0.0 / 0.6 / 0.7", + "samples": "P0 / P1 / P2", + "variant": "A22206G", + "start": 22206, + }, + { + "gene": "S", + "amino_acid_consequence": "E484K", + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "persistence_status": "new_intermittent", + "alt_freq": "0.0 / 0.6 / 0.7", + "samples": "P0 / P1 / P2", + "variant": "G23012A", + "start": 23012, + }, + { + "gene": "S", + "amino_acid_consequence": "N501Y", + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "persistence_status": "new_transient", + "alt_freq": "0.0 / 0.6 / 0.0", + "samples": "P0 / P1 / P2", + "variant": "A23063T", + "start": 23063, + }, + ] + ) + + matrix = _prepare_variant_heatmap_matrix( + table, + ["P0", "P1", "P2"], + 0.2, + 0.2, + only_persistent=True, + ) + + assert set(matrix.index) == {"S:D215G\n(A22206G)", "S:E484K\n(G23012A)"} + + +def test_prepare_variant_heatmap_matrix_default_includes_canonical_joint_only_variant(): + """When every row for a variant is joint, the default (`include_joint=False`) + heatmap must not silently drop it entirely - its (only, hence canonical) + row is exempted from the implicit joint exclusion.""" + table = pd.DataFrame( + [ + { + "gene": "PA1199", + "amino_acid_consequence": "D50G", + "nsp_aa_change": "", + "type_of_change": "joint_missense", + "type_of_variant": "snp", + "alt_freq": "0.6", + "samples": "P0", + "variant": "A100G", + "start": 100, + }, + ] + ) + + matrix = _prepare_variant_heatmap_matrix(table, ["P0"], 0.0, 0.0) + + assert list(matrix.index) == ["PA1199:D50G\n(A100G)"] + + +def test_prepare_variant_heatmap_matrix_default_hides_noncanonical_joint_fragment(): + """A variant with both a canonical (non-joint) row and a joint fragment + row must still hide the joint fragment by default - the canonical-row + exemption only rescues variants whose ONLY row is joint, it must not + let every joint row through.""" + table = pd.DataFrame( + [ + { + "gene": "S", + "amino_acid_consequence": "D50G", + "nsp_aa_change": "", + "type_of_change": "missense", + "type_of_variant": "snp", + "presence_absence": "Y / Y / Y", + "alt_freq": "0.5 / 0.6 / 0.7", + "samples": "P0 / P1 / P2", + "variant": "A100G", + "start": 100, + }, + { + "gene": "S", + "amino_acid_consequence": "D50G+X60Y", + "nsp_aa_change": "", + "type_of_change": "joint_missense", + "type_of_variant": "snp", + "presence_absence": "N / N / Y", + "alt_freq": "0.0 / 0.0 / 0.9", + "samples": "P0 / P1 / P2", + "variant": "A100G", + "start": 100, + }, + ] + ) + + matrix = _prepare_variant_heatmap_matrix(table, ["P0", "P1", "P2"], 0.0, 0.0) + + assert list(matrix.index) == ["S:D50G\n(A100G)"] + assert matrix.loc["S:D50G\n(A100G)", "P2"] == 0.7 + + def test_process_joint_variants_only_adds_single_joint_prefix(tmp_path): csv_path = tmp_path / "results.csv" pd.DataFrame( @@ -477,6 +976,40 @@ def test_process_joint_variants_matches_main_row_by_presence_pattern(tmp_path): assert result.loc[2, "type_of_change"] == "joint_missense" +def test_process_joint_variants_flags_compound_row_without_at_pointer(tmp_path): + """A row whose own `bcsq_nt_notation` is genuinely compound (multiple + '+'-joined `pos ref>alt` terms) must be flagged joint even when no + `@`-pointer row in the table points at it - the one-to-one `@`-pointer + matching only ever touches one sibling row per position, so a compound + sibling that nothing points to would otherwise be missed entirely.""" + csv_path = tmp_path / "results.csv" + pd.DataFrame( + [ + { + "start": 100, + "gene": "PA3025", + "amino_acid_consequence": "CLLL", + "nsp_aa_change": "", + "bcsq_nt_notation": "100A>T+150C>G", + "bcsq_aa_notation": "p.CLLL", + "aa1_total_properties": "", + "aa2_total_properties": "", + "aa1_unique_properties": "", + "aa2_unique_properties": "", + "aa1_weight": "", + "aa2_weight": "", + "weight_difference": "", + "type_of_change": "missense", + }, + ] + ).to_csv(csv_path, index=False) + + result = process_joint_variants(str(csv_path)) + + assert bool(result.loc[0, "joint_variant"]) is True + assert result.loc[0, "type_of_change"] == "joint_missense" + + def test_process_joint_variants_is_order_invariant_for_overlapping_gene_rows(tmp_path): shared_rows = [ { @@ -649,3 +1182,93 @@ def test_heatmap_figure_size_enforces_minimum_row_height(): assert width == 4.0 assert height == pytest.approx(5.2) + + +def test_select_canonical_row_positions_prefers_union_over_non_joint(): + """When every row for a variant is joint, the row carrying the true + (union) trajectory must win even though every candidate is joint - this + is the PA1199-shaped case where "prefer non-joint" alone would pick a + fragment row with the wrong trajectory.""" + table = pd.DataFrame( + [ + { + "chrom": "NC_002516.2", + "start": 1300712, + "ref": "AG", + "alt": "A", + "presence_absence": "N / N / Y / N / N / N", + "joint_variant": True, + "amino_acid_consequence": "short", + "type_of_change": "joint_stop_gained&frameshift", + }, + { + "chrom": "NC_002516.2", + "start": 1300712, + "ref": "AG", + "alt": "A", + "presence_absence": "N / N / N / Y / N / N", + "joint_variant": True, + "amino_acid_consequence": "short", + "type_of_change": "joint_missense&inframe_altering", + }, + { + "chrom": "NC_002516.2", + "start": 1300712, + "ref": "AG", + "alt": "A", + "presence_absence": "N / N / Y / Y / Y / Y", + "joint_variant": True, + "amino_acid_consequence": "short", + "type_of_change": "joint_stop_gained&frameshift", + }, + ] + ) + keys = [ + _build_variant_key(row, "fallback") for row in table.itertuples(index=False) + ] + + canonical = select_canonical_row_positions(table, keys) + + assert len(canonical) == 1 + (position,) = canonical.values() + assert table.iloc[position]["presence_absence"] == "N / N / Y / Y / Y / Y" + + +def test_select_canonical_row_positions_prefers_standalone_when_it_covers_the_union(): + """The pre-existing viral case: a standalone (non-joint) row already + carries the full trajectory alongside a joint row masked to a subset of + samples - the standalone row must still win.""" + table = pd.DataFrame( + [ + { + "chrom": "segA", + "start": 100, + "ref": "A", + "alt": "T", + "presence_absence": "Y / Y / Y", + "joint_variant": False, + "amino_acid_consequence": "D614G", + "type_of_change": "missense", + }, + { + "chrom": "segA", + "start": 100, + "ref": "A", + "alt": "T", + "presence_absence": "N / N / Y", + "joint_variant": True, + "amino_acid_consequence": "D614G+N501Y", + "type_of_change": "joint_missense", + }, + ] + ) + keys = [ + _build_variant_key(row, "fallback") for row in table.itertuples(index=False) + ] + + canonical = select_canonical_row_positions(table, keys) + + assert len(canonical) == 1 + (position,) = canonical.values() + assert table.iloc[position]["presence_absence"] == "Y / Y / Y" + assert bool(table.iloc[position]["joint_variant"]) is False diff --git a/tests/test_annotation_processing.py b/tests/test_annotation_processing.py new file mode 100644 index 0000000..f2bdd74 --- /dev/null +++ b/tests/test_annotation_processing.py @@ -0,0 +1,53 @@ +"""Tests for annotation_processing helpers, especially multi-contig handling.""" + +from __future__ import annotations + +from vartracker.annotation_processing import ( + ambiguous_gene_names, + gene_lengths_from_gff3, +) + +SINGLE_CONTIG_GFF3 = """\ +##gff-version 3 +chrom1\tcustom\tCDS\t1\t300\t.\t+\t0\tID=cds-dnaA;gene=dnaA +chrom1\tcustom\tCDS\t400\t600\t.\t+\t0\tID=cds-repA;gene=repA +""" + +MULTI_CONTIG_COLLISION_GFF3 = """\ +##gff-version 3 +chrom1\tcustom\tCDS\t1\t300\t.\t+\t0\tID=cds-dnaA;gene=dnaA +chrom1\tcustom\tCDS\t400\t600\t.\t+\t0\tID=cds-repA-chrom;gene=repA +plasmid1\tcustom\tCDS\t1\t150\t.\t+\t0\tID=cds-repA-plasmid;gene=repA +plasmid1\tcustom\tCDS\t200\t350\t.\t+\t0\tID=cds-mobA;gene=mobA +""" + + +def test_gene_lengths_single_contig_unqualified(tmp_path): + """Single-contig references (e.g. viral genomes) keep plain gene names.""" + gff_path = tmp_path / "single.gff3" + gff_path.write_text(SINGLE_CONTIG_GFF3) + + lengths = gene_lengths_from_gff3(gff_path) + + assert lengths["dnaA"] == 300 + assert lengths["repA"] == 201 + assert ambiguous_gene_names(gff_path) == set() + + +def test_gene_lengths_disambiguates_colliding_names_across_contigs(tmp_path): + """Identically-named genes on different contigs/replicons must not be summed.""" + gff_path = tmp_path / "multi.gff3" + gff_path.write_text(MULTI_CONTIG_COLLISION_GFF3) + + lengths = gene_lengths_from_gff3(gff_path) + + # Colliding gene name is split per-contig, each keeping its own length. + assert lengths["repA (chrom1)"] == 201 + assert lengths["repA (plasmid1)"] == 150 + assert "repA" not in lengths + + # Non-colliding gene names on either contig stay unqualified. + assert lengths["dnaA"] == 300 + assert lengths["mobA"] == 151 + + assert ambiguous_gene_names(gff_path) == {"repA"} diff --git a/tests/test_cli.py b/tests/test_cli.py index 3f29f0c..40db029 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -34,7 +34,8 @@ def test_plot_heatmap_help_uses_standalone_option_names(): assert "--x-labels" in heatmap_options assert "--heatmap-aa-exclude" not in out assert "--heatmap-include-joint" not in out - assert "--outdir" not in out + assert "--out" in out + assert "--outdir" in out assert "--name" not in out assert "--min-snv-freq" not in out assert "--min-indel-freq" not in out diff --git a/tests/test_main.py b/tests/test_main.py index 3885c39..5d65799 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -1080,6 +1080,79 @@ def fake_vcf(args): assert modes_checked == ["e2e"] +def _prepare_e2e_run(monkeypatch, tmp_path): + updated_csv = tmp_path / "samples_updated.csv" + vcf_out = tmp_path / "vcf.gz" + cov_out = tmp_path / "coverage.txt" + vcf_out.write_text("", encoding="utf-8") + cov_out.write_text("", encoding="utf-8") + updated_csv.write_text( + "sample_name,sample_number,reads1,reads2,bam,vcf,coverage\n" + f"Sample1,0,,,,{vcf_out},{cov_out}\n", + encoding="utf-8", + ) + + monkeypatch.setattr(main_module, "validate_dependencies", lambda mode="vcf": None) + monkeypatch.setattr( + main_module, "run_e2e_workflow", lambda **kwargs: str(updated_csv) + ) + monkeypatch.setattr(main_module, "_run_vcf_command", lambda args: 0) + + reads1 = tmp_path / "reads1.fastq" + reads2 = tmp_path / "reads2.fastq" + reads1.write_text("", encoding="utf-8") + reads2.write_text("", encoding="utf-8") + samples_csv = tmp_path / "reads.csv" + samples_csv.write_text( + "sample_name,sample_number,reads1,reads2,bam,vcf,coverage\n" + f"Sample1,0,{reads1},{reads2},,,\n", + encoding="utf-8", + ) + return samples_csv + + +def test_e2e_warns_when_primer_bed_missing(monkeypatch, tmp_path, capsys): + samples_csv = _prepare_e2e_run(monkeypatch, tmp_path) + + exit_code = main_module.main( + [ + "end-to-end", + str(samples_csv), + "--reference", + "ref.fasta", + "--outdir", + str(tmp_path / "results"), + ] + ) + + assert exit_code == 0 + captured = capsys.readouterr() + assert "no --primer-bed supplied" in captured.out + + +def test_e2e_no_warning_when_primer_bed_supplied(monkeypatch, tmp_path, capsys): + samples_csv = _prepare_e2e_run(monkeypatch, tmp_path) + primer_bed = tmp_path / "primers.bed" + primer_bed.write_text("", encoding="utf-8") + + exit_code = main_module.main( + [ + "end-to-end", + str(samples_csv), + "--reference", + "ref.fasta", + "--primer-bed", + str(primer_bed), + "--outdir", + str(tmp_path / "results"), + ] + ) + + assert exit_code == 0 + captured = capsys.readouterr() + assert "no --primer-bed supplied" not in captured.out + + def test_e2e_dryrun_skips_vcf(monkeypatch, tmp_path_factory): reads1 = tmp_path_factory.mktemp("reads") / "reads1.fastq" reads2 = reads1.parent / "reads2.fastq" diff --git a/tests/test_plotting.py b/tests/test_plotting.py index fb3b200..8b9190f 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -142,6 +142,81 @@ def test_prepare_plot_inputs_and_filters(): assert filtered_long["sample_number"].max() == 3 +def test_apply_shared_plot_filters_persistent_only_includes_new_intermittent(): + """persistent_only must keep new_intermittent alongside new_persistent - + both are new variants present at the final timepoint, just differing in + whether presence was continuous along the way.""" + table = pd.DataFrame( + [ + { + "chrom": "segA", + "start": 23403, + "end": 23403, + "gene": "S", + "variant": "A23403G", + "amino_acid_consequence": "S:D614G", + "nsp_aa_change": "", + "type_of_variant": "snp", + "type_of_change": "missense", + "variant_status": "new", + "persistence_status": "new_persistent", + "presence_absence": "N / Y / Y / Y", + "alt_freq": "0.0 / 0.20 / 0.45 / 0.70", + "samples": "P0 / P1 / P2 / P3", + "sample_number": "0 / 1 / 2 / 3", + "per_sample_variant_qc": "P / P / P / P", + "reference": "", + }, + { + "chrom": "segA", + "start": 23012, + "end": 23012, + "gene": "S", + "variant": "G23012A", + "amino_acid_consequence": "S:E484K", + "nsp_aa_change": "", + "type_of_variant": "snp", + "type_of_change": "missense", + "variant_status": "new", + "persistence_status": "new_intermittent", + "presence_absence": "N / Y / N / Y", + "alt_freq": "0.0 / 0.20 / 0.0 / 0.45", + "samples": "P0 / P1 / P2 / P3", + "sample_number": "0 / 1 / 2 / 3", + "per_sample_variant_qc": "P / P / P / P", + "reference": "", + }, + { + "chrom": "segA", + "start": 23063, + "end": 23063, + "gene": "S", + "variant": "A23063T", + "amino_acid_consequence": "S:N501Y", + "nsp_aa_change": "", + "type_of_variant": "snp", + "type_of_change": "missense", + "variant_status": "new", + "persistence_status": "new_transient", + "presence_absence": "N / Y / N / N", + "alt_freq": "0.0 / 0.20 / 0.0 / 0.0", + "samples": "P0 / P1 / P2 / P3", + "sample_number": "0 / 1 / 2 / 3", + "per_sample_variant_qc": "P / P / P / P", + "reference": "", + }, + ] + ) + + summary, long_df, _, _ = prepare_plot_inputs(table) + filtered_summary, _ = apply_shared_plot_filters( + summary, long_df, persistent_only=True + ) + + labels = {label.split(" (", 1)[0] for label in filtered_summary["variant_label"]} + assert labels == {"S:D614G", "S:E484K"} + + def test_auto_select_variants_prefers_literature_and_persistent(): summary, _, _, _ = prepare_plot_inputs(_results_table()) @@ -290,17 +365,241 @@ def test_prepare_plot_inputs_uses_unique_internal_variant_keys(): summary, long_df, _, _ = prepare_plot_inputs(table) assert summary["variant_id"].nunique() == 2 + # Both variants would otherwise render to the same "INTERGENIC:None" + # display label. variant_id (the internal key) was always unique for + # this case; now the *displayed* label is also disambiguated (via + # find_colliding_base_labels/disambiguate_base_label) with a positional + # + ref>alt suffix, so the two variants are distinguishable on plots too. assert {label.split(" (", 1)[0] for label in summary["variant_label"]} == { - "INTERGENIC:None" + "INTERGENIC:None @100A>G", + "INTERGENIC:None @200C>T", } assert ( long_df.loc[ - long_df["variant_label"].str.startswith("INTERGENIC:None"), "variant_id" + long_df["variant_label"].str.startswith("INTERGENIC:None @"), + "variant_id", ].nunique() == 2 ) +def test_prepare_plot_inputs_collapses_duplicate_rows_for_trajectory(): + # Models the real PA1199 case from PAO1 validation + # (NC_002516.2:1300712 AG>A): one physical variant described by + # bcftools-csq under three distinct joint annotation groups, each row + # masked to a different single sample, collectively present in all 3 + # samples even though no single row shows that on its own. + base = { + "chrom": "NC_002516.2", + "start": 1300712, + "end": 1300713, + "ref": "AG", + "alt": "A", + "gene": "PA1199", + "variant": "1300712_AGdelG", + "nsp_aa_change": "", + "type_of_variant": "indel", + "variant_status": "new", + "persistence_status": "new_persistent", + "samples": "P1 / P2 / P3", + "sample_number": "1 / 2 / 3", + "per_sample_variant_qc": "P / P / P", + "joint_variant": True, + } + table = pd.DataFrame( + [ + { + **base, + "amino_acid_consequence": "PA1199:frameshift_group_A", + "type_of_change": "joint_frameshift", + "presence_absence": "Y / N / N", + "alt_freq": "0.50 / 0.03 / 0.02", + }, + { + **base, + "amino_acid_consequence": "PA1199:frameshift_group_B", + "type_of_change": "joint_frameshift", + "presence_absence": "N / Y / N", + "alt_freq": "0.03 / 0.60 / 0.04", + }, + { + **base, + "amino_acid_consequence": "PA1199:frameshift_group_C", + "type_of_change": "joint_frameshift", + "presence_absence": "N / N / Y", + "alt_freq": "0.02 / 0.05 / 0.70", + }, + ] + ) + + summary, long_df, _, _ = prepare_plot_inputs(table) + + assert long_df["variant_id"].nunique() == 1 + variant_id = long_df["variant_id"].iloc[0] + variant_rows = long_df[long_df["variant_id"] == variant_id].sort_values( + "sample_number" + ) + + assert len(variant_rows) == 3 + assert list(variant_rows["sample_number"]) == [1, 2, 3] + assert variant_rows["present"].all() + # Each sample's AF must be the max across the variant's 3 source rows for + # that sample, not zero and not an arbitrary single row's value. + assert list(variant_rows["allele_frequency"].round(2)) == [0.50, 0.60, 0.70] + + +def test_prepare_plot_inputs_collapses_duplicate_rows_for_turnover_counts(): + # Same fragmented-row scenario as above, but shaped for a clear + # new/lost/new transition sequence across 4 samples so turnover counting + # can be checked against a known-correct answer. + base = { + "chrom": "NC_002516.2", + "start": 1300712, + "end": 1300713, + "ref": "AG", + "alt": "A", + "gene": "PA1199", + "variant": "1300712_AGdelG", + "amino_acid_consequence": "PA1199:frameshift", + "nsp_aa_change": "", + "type_of_variant": "indel", + "type_of_change": "joint_frameshift", + "variant_status": "new", + "persistence_status": "new_transient", + "samples": "P1 / P2 / P3 / P4", + "sample_number": "1 / 2 / 3 / 4", + "per_sample_variant_qc": "P / P / P / P", + "joint_variant": True, + } + table = pd.DataFrame( + [ + { + **base, + "presence_absence": "Y / N / N / N", + "alt_freq": "0.40 / 0.0 / 0.0 / 0.0", + }, + { + **base, + "presence_absence": "N / N / Y / Y", + "alt_freq": "0.0 / 0.0 / 0.30 / 0.50", + }, + ] + ) + + summary, long_df, _, _ = prepare_plot_inputs(table) + + # No duplicate (variant_id, sample_number) pairs: exactly one row per + # sample for this variant, as plot_variant_turnover assumes. + assert not long_df.duplicated(subset=["variant_id", "sample_number"]).any() + + variant_id = long_df["variant_id"].iloc[0] + trajectory = long_df[long_df["variant_id"] == variant_id].sort_values( + "sample_number" + ) + presence_sequence = list(trajectory["present"]) + # Union of the two rows' presence vectors: present, absent, present, present. + assert presence_sequence == [True, False, True, True] + + # Naive new/lost transition scan mirroring what plot_variant_turnover + # does over a collapsed long_df (one row per variant per sample). + transitions = 0 + previous_present = False + for present in presence_sequence: + if present != previous_present: + transitions += 1 + previous_present = present + # new at sample 1, lost at sample 2, new (re-emergence) at sample 3. + assert transitions == 3 + + +def test_prepare_plot_inputs_summary_uses_canonical_row_metadata(): + variant_a_base = { + "chrom": "NC_002516.2", + "start": 2000000, + "end": 2000000, + "ref": "C", + "alt": "T", + "gene": "PA2000", + "variant": "2000000C>T", + "nsp_aa_change": "", + "type_of_variant": "snp", + "variant_status": "new", + "persistence_status": "new_persistent", + "samples": "P1 / P2 / P3", + "sample_number": "1 / 2 / 3", + "per_sample_variant_qc": "P / P / P", + } + variant_b_base = { + "chrom": "NC_002516.2", + "start": 3000000, + "end": 3000000, + "ref": "G", + "alt": "A", + "gene": "PA3000", + "variant": "3000000G>A", + "nsp_aa_change": "", + "type_of_variant": "snp", + "variant_status": "new", + "persistence_status": "new_persistent", + "samples": "P1 / P2 / P3", + "sample_number": "1 / 2 / 3", + "per_sample_variant_qc": "P / P / P", + } + + table = pd.DataFrame( + [ + # Variant A: one joint fragment row (masked, NOT the union) and + # one standalone row that IS the union -> standalone wins + # because it's both the union match AND non-joint. + { + **variant_a_base, + "amino_acid_consequence": "PA2000:joint_frag", + "type_of_change": "joint_missense", + "presence_absence": "N / N / Y", + "alt_freq": "0.0 / 0.0 / 0.30", + "joint_variant": True, + }, + { + **variant_a_base, + "amino_acid_consequence": "PA2000:D10E", + "type_of_change": "missense", + "presence_absence": "Y / Y / Y", + "alt_freq": "0.20 / 0.25 / 0.30", + "joint_variant": False, + }, + # Variant B: both rows are joint; one is a masked fragment, the + # other is the union -> the union row wins even though it is + # also joint (no non-joint candidate exists to prefer instead). + { + **variant_b_base, + "amino_acid_consequence": "PA3000:joint_frag_x", + "type_of_change": "joint_missense_x", + "presence_absence": "N / N / Y", + "alt_freq": "0.0 / 0.0 / 0.40", + "joint_variant": True, + }, + { + **variant_b_base, + "amino_acid_consequence": "PA3000:joint_frag_y", + "type_of_change": "joint_missense_y", + "presence_absence": "Y / Y / Y", + "alt_freq": "0.10 / 0.20 / 0.40", + "joint_variant": True, + }, + ] + ) + + summary, long_df, _, _ = prepare_plot_inputs(table) + + row_a = summary[summary["gene"] == "PA2000"] + assert len(row_a) == 1 + assert row_a["type_of_change"].iloc[0] == "missense" + + row_b = summary[summary["gene"] == "PA3000"] + assert len(row_b) == 1 + assert row_b["type_of_change"].iloc[0] == "joint_missense_y" + + def test_genome_plot_defaults_to_snps_only(capsys): table = pd.DataFrame( [ diff --git a/tests/test_vcf_processing.py b/tests/test_vcf_processing.py index 6816f40..d6468ed 100644 --- a/tests/test_vcf_processing.py +++ b/tests/test_vcf_processing.py @@ -14,6 +14,7 @@ from vartracker.analysis import process_joint_variants from vartracker.vcf_processing import ( _derive_vcf_output_paths, + _summarise_sample_trajectory, annotate_vcf, format_vcf, merge_consequences, @@ -32,6 +33,71 @@ def test_derive_vcf_output_paths_handles_gz(tmp_path): assert os.path.dirname(out) == str(tmp_path) +def test_summarise_sample_trajectory_original_retained_no_gap(): + result = _summarise_sample_trajectory(["0.1", "0.1", "0.1"], ["s0", "s1", "s2"]) + assert result["variant_status"] == "original" + assert result["persistent_status"] == "original_retained" + assert result["presence_absence"] == ["Y", "Y", "Y"] + + +def test_summarise_sample_trajectory_original_intermittent_has_gap(): + """Present at first and last timepoints, but absent in between - this is + the "original" side of the reviewer's ambiguity, distinct from + continuous presence.""" + result = _summarise_sample_trajectory(["0.1", ".", "0.1"], ["s0", "s1", "s2"]) + assert result["variant_status"] == "original" + assert result["persistent_status"] == "original_intermittent" + assert result["presence_absence"] == ["Y", "N", "Y"] + assert result["first_appearance"] == "s0" + assert result["last_appearance"] == "s2" + + +def test_summarise_sample_trajectory_original_lost(): + result = _summarise_sample_trajectory(["0.1", "0.1", "."], ["s0", "s1", "s2"]) + assert result["persistent_status"] == "original_lost" + + +def test_summarise_sample_trajectory_new_persistent_no_gap(): + result = _summarise_sample_trajectory([".", "0.1", "0.1"], ["s0", "s1", "s2"]) + assert result["variant_status"] == "new" + assert result["persistent_status"] == "new_persistent" + + +def test_summarise_sample_trajectory_new_intermittent_has_gap(): + """The exact scenario the reviewer flagged: a new mutation that appears, + disappears in a subsequent passage, then reappears by the final + timepoint. Previously indistinguishable from continuous persistence.""" + result = _summarise_sample_trajectory( + [".", "0.1", ".", "0.1"], ["s0", "s1", "s2", "s3"] + ) + assert result["variant_status"] == "new" + assert result["persistent_status"] == "new_intermittent" + assert result["first_appearance"] == "s1" + assert result["last_appearance"] == "s3" + + +def test_summarise_sample_trajectory_new_transient(): + result = _summarise_sample_trajectory([".", "0.1", "."], ["s0", "s1", "s2"]) + assert result["persistent_status"] == "new_transient" + + +def test_summarise_sample_trajectory_absent_throughout_is_new_transient(): + result = _summarise_sample_trajectory([".", ".", "."], ["s0", "s1", "s2"]) + assert result["persistent_status"] == "new_transient" + assert result["first_appearance"] == "None" + assert result["last_appearance"] == "None" + + +def test_summarise_sample_trajectory_leading_absence_is_not_a_gap(): + """A run of absence *before* the first appearance of a new variant is + expected (that's just what "new" means) and must not itself be treated + as a gap.""" + result = _summarise_sample_trajectory( + [".", ".", "0.1", "0.1"], ["s0", "s1", "s2", "s3"] + ) + assert result["persistent_status"] == "new_persistent" + + def test_derive_vcf_output_paths_handles_plain_vcf(tmp_path): out, csq, log = _derive_vcf_output_paths("/data/alpha.vcf", tmp_path, "alpha") @@ -447,6 +513,34 @@ def test_process_vcf_merges_starred_and_unstarred_equivalent_bcsq(tmp_path): assert row["persistence_status"] == "original_retained" +def test_process_vcf_decodes_percent_encoded_gene_name_in_bcsq(tmp_path): + vcf_path = tmp_path / "annotated.vcf" + vcf_path.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\n" + "chr1\t4\t.\tA\tG\t.\tPASS\tDP=100;AF=0.1;" + "BCSQ=missense|ercS%27|tx|protein_coding|+|10K>10A|100A>G" + "\tGT:DP:AF:BCSQ\t1:100:0.1:1\n", + encoding="utf-8", + ) + + cov1 = tmp_path / "s1.depth.txt" + _write_depth_file(cov1) + + table = process_vcf(str(vcf_path), [str(cov1)], 10, ["s1"]) + + assert len(table) == 1 + assert table.iloc[0]["gene"] == "ercS'" + + def test_process_vcf_keeps_single_site_variant_present_when_sample_has_joint_csq( tmp_path, ): @@ -489,6 +583,74 @@ def test_process_vcf_keeps_single_site_variant_present_when_sample_has_joint_csq assert joint["alt_freq"] == ". / 0.200" +def test_process_vcf_recovers_presence_when_no_sample_gets_a_standalone_bcsq( + tmp_path, +): + """When a variant sits in a gene dense enough that every sample's own + joint-consequence call pairs it with a *different* co-occurring + neighbour, bcftools never reports it standalone for any sample, so each + annotation-group row is masked to a different, non-overlapping subset of + samples - fragmenting a variant that is genuinely present in all three + samples across two partial rows, with neither showing its true, full + trajectory. This is the bacterial-scale hotspot-gene scenario (see the + PAO1 validation): the variant's real presence must still surface in one + row so persistence_status isn't computed from a fragment. The two joint + rows themselves are kept unconditionally (see the module-level comment + in process_vcf on why row count is never reduced by deleting them) - + a single-sample joint call can be genuine, novel biology, not just an + artefact, and there is no reliable way to tell the two apart from + presence data alone.""" + vcf_path = tmp_path / "annotated.vcf" + vcf_path.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2\ts3\n" + "chr1\t4\t.\tA\tG\t.\tPASS\tDP=100;AF=0.2;" + "BCSQ=missense|GENE1|tx|protein_coding|+|2K>2A|4A>G+5A>C," + "frameshift|GENE1|tx|protein_coding|+|2KAAAAAAAAAAAA>2K|4A>G+6A>T" + "\tGT:DP:AF:BCSQ\t1:100:0.1:1\t1:100:0.2:4\t1:100:0.3:0\n", + encoding="utf-8", + ) + + cov1 = tmp_path / "s1.depth.txt" + cov2 = tmp_path / "s2.depth.txt" + cov3 = tmp_path / "s3.depth.txt" + _write_depth_file(cov1) + _write_depth_file(cov2) + _write_depth_file(cov3) + + table = process_vcf( + str(vcf_path), [str(cov1), str(cov2), str(cov3)], 10, ["s1", "s2", "s3"] + ) + + # Both joint-detail rows survive, each masked to the one sample whose + # own bitmask matched that specific compound description. + joint_with_5 = table[table["bcsq_nt_notation"].eq("4A>G+5A>C")].iloc[0] + assert joint_with_5["presence_absence"] == "Y / N / N" + + joint_with_6 = table[table["bcsq_nt_notation"].eq("4A>G+6A>T")].iloc[0] + assert joint_with_6["presence_absence"] == "N / Y / N" + + # No group's bitmask matched s3, yet s3 genuinely carries the ALT allele + # (GT=1, AF=0.3) - the fix must add one extra row from the true, + # unmasked per-sample presence so this isn't silently dropped. + recovered = table[table["presence_absence"].eq("Y / Y / Y")] + assert len(recovered) == 1 + recovered_row = recovered.iloc[0] + assert recovered_row["alt_freq"] == "0.100 / 0.200 / 0.300" + assert recovered_row["variant_status"] == "original" + assert recovered_row["persistence_status"] == "original_retained" + + assert len(table) == 3 + + @pytest.mark.skipif(shutil.which("bcftools") is None, reason="bcftools not available") def test_merge_then_annotate_preserves_joint_annotations_across_samples(tmp_path): ref, gff = _write_minimal_reference_bundle(tmp_path) diff --git a/vartracker/_version.py b/vartracker/_version.py index 11e633a..e463cb7 100644 --- a/vartracker/_version.py +++ b/vartracker/_version.py @@ -2,18 +2,40 @@ from __future__ import annotations +from pathlib import Path + try: # Python 3.8+ from importlib import metadata except ImportError: # pragma: no cover -- fallback for very old Pythons import importlib_metadata as metadata # type: ignore -# NOTE: When bumping the project version remember to update this fallback value -# alongside the version declared in pyproject.toml. -_FALLBACK_VERSION = "2.2.1" +_LAST_RESORT_VERSION = "0.0.0+unknown" + + +def _read_pyproject_version() -> str | None: + """Read `project.version` straight from pyproject.toml. + + Only reachable when vartracker isn't pip-installed (no dist-info to read + metadata from), e.g. a git checkout run in place - so pyproject.toml is + still sitting one directory above this file. Keeps the version fallback + in sync with pyproject.toml automatically, rather than a second + hand-maintained copy of the version string. + """ + pyproject_path = Path(__file__).resolve().parent.parent / "pyproject.toml" + if not pyproject_path.exists(): + return None + try: + import tomllib + except ImportError: # pragma: no cover -- Python <3.11 + return None + with pyproject_path.open("rb") as handle: + data = tomllib.load(handle) + return data.get("project", {}).get("version") + try: __version__ = metadata.version("vartracker") except metadata.PackageNotFoundError: # pragma: no cover - during local dev - __version__ = _FALLBACK_VERSION + __version__ = _read_pyproject_version() or _LAST_RESORT_VERSION __all__ = ["__version__"] diff --git a/vartracker/analysis.py b/vartracker/analysis.py index e895d93..5a32cac 100644 --- a/vartracker/analysis.py +++ b/vartracker/analysis.py @@ -19,6 +19,7 @@ import seaborn as sns from .constants import REF_GENE_LENGTHS, NSP_LENGTHS, NSPS from .core import get_logo +from .vcf_processing import NEW_ENDS_PRESENT_STATUSES # Global configuration for plotting plt.rcdefaults() @@ -100,6 +101,130 @@ def _presence_vector(value: object) -> tuple[str, ...]: return tuple(_parse_slash_separated_tokens(value)) +def _build_variant_key(row: object, fallback_label: str) -> str: + chrom = str(getattr(row, "chrom", "")).strip() + start = str(getattr(row, "start", "")).strip() + ref = str(getattr(row, "ref", "")).strip() + alt = str(getattr(row, "alt", "")).strip() + variant_name = str(getattr(row, "variant", "")).strip() + + if chrom and start and ref and alt: + return f"{chrom}:{start}:{ref}>{alt}" + if chrom and start and variant_name: + return f"{chrom}:{start}:{variant_name}" + if variant_name: + return variant_name + return fallback_label + + +def _row_is_joint(row: object) -> bool: + """Best-effort joint/compound flag for a results-table row. + + Prefers the explicit `joint_variant` column (set by + `process_joint_variants`); falls back to the `type_of_change` prefix for + tables that predate that column (e.g. synthetic test fixtures). + """ + joint_variant = getattr(row, "joint_variant", None) + if joint_variant is not None and not ( + isinstance(joint_variant, float) and pd.isna(joint_variant) + ): + return bool(joint_variant) + change_type = str(getattr(row, "type_of_change", "")).strip().lstrip("*") + return re.sub(r"^(joint_)+", "", change_type) != change_type + + +def select_canonical_row_positions( + table: pd.DataFrame, variant_keys: Sequence[str] +) -> dict[str, int]: + """Map each variant key to the positional index of its canonical row. + + `variant_keys[i]` is the caller's key for the i-th row of `table` (as + produced by `_build_variant_key`), so callers that already build keys in + their own pass cannot drift from this one. + + A single physical variant can have more than one row when bcftools csq + describes it under several distinct joint/compound annotation groups + (see the row-canonicalisation logic in `vcf_processing.py::process_vcf`). + The canonical row is chosen, in order: + + 1. Among rows whose `presence_absence` equals the element-wise union of + all of the variant's rows' presence vectors (i.e. the row - or rows - + that alone already carry the variant's true, full trajectory). If no + row matches the union exactly (not expected, but cheap to guard), + fall back to the row(s) with the most "Y" tokens. + 2. Among those, prefer a non-joint row (see `_row_is_joint`) - this + reproduces today's behaviour exactly for the standalone-vs-joint-detail + case that already existed before row-canonicalisation. + 3. Tie-break deterministically: shortest `amino_acid_consequence` (least + compound description), then lowest positional index. + + Returns a dict of variant_key -> positional index into `table` (i.e. an + index usable with `table.iloc[...]`, not `table.loc[...]`). + """ + positions_by_key: dict[str, list[int]] = {} + for position, key in enumerate(variant_keys): + positions_by_key.setdefault(key, []).append(position) + + canonical: dict[str, int] = {} + for key, positions in positions_by_key.items(): + if len(positions) == 1: + canonical[key] = positions[0] + continue + + presence_vectors = { + position: _presence_vector( + getattr(table.iloc[position], "presence_absence", "") + ) + for position in positions + } + length = max((len(vec) for vec in presence_vectors.values()), default=0) + union = tuple( + ( + "Y" + if any( + vec[idx] == "Y" + for vec in presence_vectors.values() + if idx < len(vec) + ) + else "N" + ) + for idx in range(length) + ) + + candidates = [ + position for position in positions if presence_vectors[position] == union + ] + if not candidates: + max_present = max( + sum(token == "Y" for token in vec) for vec in presence_vectors.values() + ) + candidates = [ + position + for position in positions + if sum(token == "Y" for token in presence_vectors[position]) + == max_present + ] + + rows_by_position = {position: table.iloc[position] for position in candidates} + non_joint_candidates = [ + position + for position in candidates + if not _row_is_joint(rows_by_position[position]) + ] + if non_joint_candidates: + candidates = non_joint_candidates + + def _tie_break_key(position: int) -> tuple[int, int]: + aa_consequence = str( + getattr(rows_by_position[position], "amino_acid_consequence", "") + ) + return (len(aa_consequence), position) + + canonical[key] = min(candidates, key=_tie_break_key) + + return canonical + + def _find_joint_main_index(tab: pd.DataFrame, joint_index: int, main_pos: int) -> int: candidates = tab[ (tab["start"] == main_pos) @@ -163,10 +288,6 @@ def process_joint_variants(path): # Get indices where bcsq_aa_notation starts with "@" idx = tab[tab["bcsq_aa_notation"].str.startswith("@", na=False)].index - # If no joint variants, return the table as is - if idx.empty: - return tab - for i in idx: try: # Find the main pos and main index number @@ -206,6 +327,21 @@ def process_joint_variants(path): print(f"Warning: Could not process joint variant at index {i}: {str(e)}") continue + # Second, additive pass: some rows carry a genuinely compound + # `bcsq_nt_notation` (e.g. "190263TCC>T+190295T>G+190354G>A") - more than + # one "pos ref>alt" term joined by "+" - but were never flagged joint by + # the `@`-pointer matching above. This happens when a variant has several + # sibling rows at the same position (routine now that rows are + # canonicalised), so the one-to-one `@`-pointer match only touches one + # sibling even though every sibling's own notation is plainly compound. + # Flag any such row directly from its own notation, so `joint_variant` + # and the `joint_` `type_of_change` prefix never disagree with what the + # row itself encodes. + compound_nt_mask = tab["bcsq_nt_notation"].astype(str).str.count(">") > 1 + for i in tab.index[compound_nt_mask]: + tab.at[i, "joint_variant"] = True + tab.at[i, "type_of_change"] = _ensure_joint_prefix(tab.at[i, "type_of_change"]) + tab.to_csv(path, index=None) return tab @@ -282,13 +418,23 @@ def generate_cumulative_lineplot(table, pname, sample_number_list, outname): def generate_gene_table( - table: pd.DataFrame, gene_lengths: Dict[str, int] | None = None + table: pd.DataFrame, + gene_lengths: Dict[str, int] | None = None, + ambiguous_genes: set | None = None, ): """ Generate gene-wise mutation statistics table. Args: table (pd.DataFrame): Variants table + gene_lengths (dict, optional): Per-gene CDS lengths, e.g. from + `gene_lengths_from_gff3`. + ambiguous_genes (set, optional): Gene names known (from the reference + annotation, e.g. via `ambiguous_gene_names`) to occur on more than + one contig/replicon. These are always split per-contig even if a + given results table only contains variants for one of the + colliding copies, so real data for the other copy is never + silently dropped when the gene-length scaffold is merged in. Returns: pd.DataFrame: Gene statistics table @@ -353,15 +499,28 @@ def make_gene_rows(gene, df): return pd.DataFrame(gene_result) - # Process original genes + # Process original genes. Genes sharing a name across more than one contig + # (e.g. bacterial plasmids reusing names like 'repA') are disambiguated by + # contig so their variant counts are never silently summed together; genes + # confined to a single contig keep their plain name, matching prior + # behaviour for single-contig/viral references and multi-segment genomes + # with unique per-segment gene names. result = [] genes_in_table = table["gene"].unique() + gene_chrom_counts = ( + table.groupby("gene")["chrom"].nunique() if "chrom" in table.columns else {} + ) + known_ambiguous = ambiguous_genes or set() for gene in genes_in_table: if gene == "INTERGENIC": continue - df = table[table["gene"] == gene] - result.append(make_gene_rows(gene, df)) + gene_df = table[table["gene"] == gene] + if gene in known_ambiguous or gene_chrom_counts.get(gene, 1) > 1: + for chrom, chrom_df in gene_df.groupby("chrom"): + result.append(make_gene_rows(f"{gene} ({chrom})", chrom_df)) + else: + result.append(make_gene_rows(gene, gene_df)) if include_nsps: table = table.copy() @@ -400,29 +559,180 @@ def assign_nsp(row): return gene_table -def plot_gene_table(gene_table, pname, outdir): +def parse_plot_genes_arg(value: str) -> List[str]: + """ + Parse a `--plot-genes` CLI value into a gene name list. + + `value` is either a comma-separated list of gene names, or a path to a + file with one gene name per line (blank lines and lines starting with + `#` are ignored). Order is preserved and duplicates are dropped. + """ + if os.path.isfile(value): + with open(value, encoding="utf-8") as handle: + raw_values = [ + line.strip() + for line in handle + if line.strip() and not line.strip().startswith("#") + ] + else: + raw_values = [token.strip() for token in value.split(",") if token.strip()] + + seen: set = set() + genes: List[str] = [] + for gene in raw_values: + if gene not in seen: + seen.add(gene) + genes.append(gene) + return genes + + +def select_genes_for_plot( + gene_table: pd.DataFrame, + max_plot_genes: Optional[int] = 30, + plot_genes: Optional[List[str]] = None, +) -> Tuple[pd.DataFrame, Optional[str]]: + """ + Select which genes appear on the gene-wise summary FIGURE. + + This is presentation-layer only: it never modifies `gene_table` itself + (the tabular/TSV output built by `generate_gene_table` always contains + every annotated gene), it only decides which subset is handed to + `plot_gene_table` for rendering. + + When `plot_genes` is not given, and the number of annotated genes + already fits within `max_plot_genes`, every gene is kept (this is what + makes SARS-CoV-2-scale runs, with their handful of genes, identical to + plotting with no cap at all). Otherwise the figure is limited to genes + that carry at least one variant, ranked by number of newly emerged + variants (total variants as tiebreak) - ranking by total variants alone + would surface long genes or reference/sample divergence rather than + anything that changed over the longitudinal series. + + Args: + gene_table: Output of `generate_gene_table`. + max_plot_genes: Cap on number of genes shown, or None for no cap. + plot_genes: Explicit gene names to show; overrides `max_plot_genes`. + Names absent from the reference annotation or with no variants + in this dataset are dropped with a warning. + + Returns: + (plot_table, subtitle): `plot_table` is the subset of `gene_table` + to plot (original gene order preserved); `subtitle` describes the + truncation applied (e.g. "top 30 of 412 genes with variants"), or + None if nothing was truncated. + """ + if gene_table.empty: + return gene_table, None + + genes_in_order = list(dict.fromkeys(gene_table["gene"])) + totals = gene_table.loc[gene_table["type"] == "total"].set_index("gene")["number"] + variant_genes = [g for g in genes_in_order if totals.get(g, 0) > 0] + + if plot_genes: + annotation_genes = set(genes_in_order) + variant_gene_set = set(variant_genes) + valid: List[str] = [] + for gene in plot_genes: + if gene not in annotation_genes: + print( + f"Warning: --plot-genes gene '{gene}' was not found in the " + "reference annotation; skipping." + ) + continue + if gene not in variant_gene_set: + print( + f"Warning: --plot-genes gene '{gene}' has no variants in " + "this dataset; skipping." + ) + continue + valid.append(gene) + + if not valid: + raise ValueError( + "None of the genes named in --plot-genes are valid: each gene " + "must be present in the reference annotation and carry at " + "least one variant in this dataset." + ) + + selected = set(valid) + plot_table = gene_table[gene_table["gene"].isin(selected)] + return plot_table, None + + if max_plot_genes is None: + return gene_table, None + if max_plot_genes <= 0: + raise ValueError( + f"--max-plot-genes must be a positive integer, got {max_plot_genes}." + ) + + if len(genes_in_order) <= max_plot_genes: + return gene_table, None + + n_with_variants = len(variant_genes) + if n_with_variants <= max_plot_genes: + plot_table = gene_table[gene_table["gene"].isin(variant_genes)] + return plot_table, None + + new_counts = gene_table.loc[gene_table["type"] == "new_mutations"].set_index( + "gene" + )["number"] + ranked = sorted( + variant_genes, key=lambda g: (-new_counts.get(g, 0), -totals.get(g, 0)) + ) + selected = set(ranked[:max_plot_genes]) + plot_table = gene_table[gene_table["gene"].isin(selected)] + subtitle = f"top {max_plot_genes} of {n_with_variants} genes with variants" + return plot_table, subtitle + + +def plot_gene_table( + gene_table, + pname, + outdir, + max_plot_genes: Optional[int] = 30, + plot_genes: Optional[List[str]] = None, +): """ Plot gene-wise mutation statistics. Args: - gene_table (pd.DataFrame): Gene statistics table + gene_table (pd.DataFrame): Gene statistics table (all annotated + genes; the tabular/TSV output is unaffected by this function). pname (str): Project name for plot title outdir (str): Output directory + max_plot_genes: Cap on number of genes shown in the FIGURE, ranked + by newly emerged variants (default: 30). Does not affect any + other output. + plot_genes: Explicit gene names to show in the FIGURE; overrides + `max_plot_genes`. """ try: + plot_table, subtitle = select_genes_for_plot( + gene_table, max_plot_genes=max_plot_genes, plot_genes=plot_genes + ) + g = sns.catplot( x="gene", y="number", col="type", col_wrap=3, - data=gene_table, + data=plot_table, kind="bar", height=4, aspect=1.2, ) - g.set_axis_labels("", "Number of Mutations") + x_label = ( + "Gene (top genes ranked by newly emerged variants; " + "total variants as tiebreak)" + if subtitle + else "" + ) + g.set_axis_labels(x_label, "Number of Mutations") g.fig.subplots_adjust(top=0.9) - g.fig.suptitle(f"{pname}", weight="bold") + title = f"{pname}" + if subtitle: + title = f"{title}\n{subtitle}" if title else subtitle + g.fig.suptitle(title, weight="bold") # Rotate labels by 90 degrees for ax in g.axes.flat: @@ -430,7 +740,7 @@ def plot_gene_table(gene_table, pname, outdir): label.set_rotation(90) # Make subplot titles nicer - for ax, title in zip(g.fig.axes, list(gene_table["type"].unique())): + for ax, title in zip(g.fig.axes, list(plot_table["type"].unique())): # Convert title to string to handle numeric values and NaN if pd.isna(title): title_str = "None" @@ -629,6 +939,99 @@ def _resolve_variant_labels(row) -> Tuple[str, str, str]: return gene_label, display_label, base_label +def find_colliding_base_labels( + base_labels: Sequence[str], variant_keys: Sequence[str] +) -> set[str]: + """Base labels that map to more than one distinct variant key. + + `base_labels` and `variant_keys` are parallel sequences, one entry per + results-table row (see `_build_variant_key` for how a variant key is + derived). The same physical variant can legitimately appear on more than + one row (e.g. once standalone and once as a joint/compound fragment) and + still share a `base_label` - that is not a collision. A collision is a + `base_label` shared by rows whose variant keys differ, i.e. genuinely + different variants that `_resolve_variant_labels` happens to render + identically. This includes, but is not limited to, truncation of long + compound amino-acid-consequence strings to the same head+N+tail form - + on real bacterial (joint-csq) data, bcftools csq can also give two + different nucleotide variants an identical, untruncated amino-acid + description. + """ + keys_by_label: dict[str, set[str]] = {} + for label, key in zip(base_labels, variant_keys): + keys_by_label.setdefault(label, set()).add(key) + return {label for label, keys in keys_by_label.items() if len(keys) > 1} + + +def disambiguate_base_label(base_label: str, row: object) -> str: + """Append a positional suffix to a base label already known to collide. + + Only call this for labels already returned by `find_colliding_base_labels` + - i.e. labels genuinely shared by more than one distinct variant. The + row's `start` is usually enough to disambiguate two otherwise + identically-rendered variants (e.g. `"PA3025:CLLL+255RPA @3388987"`); + `ref`/`alt` are folded in too whenever available, so that two distinct + variants at the exact same position (a multiallelic site) still end up + with different labels rather than requiring a second "is this still + colliding?" pass. + + Note: a disambiguated label no longer matches the literature-anchor + `_canonical_label` lookup. That's acceptable - collisions only arise for + long compound joint labels, which never matched the literature database + anyway (see `_canonical_label`). + """ + start = str(getattr(row, "start", "")).strip() + candidate = f"{base_label} @{start}" + ref = str(getattr(row, "ref", "")).strip() + alt = str(getattr(row, "alt", "")).strip() + if ref or alt: + candidate = f"{candidate}{ref}>{alt}" + return candidate + + +def compute_variant_labeling( + table: pd.DataFrame, +) -> tuple[list[str], set[str], dict[str, int]]: + """Table-wide bookkeeping shared by every plot/heatmap caller: a variant + key per row, the set of base labels that collide across distinct + variants, and each variant's canonical row position. + + Returns (variant_keys, colliding_base_labels, canonical_positions), + where variant_keys[i] is the key for table's i-th row (see + `_build_variant_key`). + """ + base_labels: list[str] = [] + variant_keys: list[str] = [] + for row in table.itertuples(index=False): + _, _, row_base_label = _resolve_variant_labels(row) + base_labels.append(row_base_label) + variant_keys.append(_build_variant_key(row, row_base_label)) + + colliding_base_labels = find_colliding_base_labels(base_labels, variant_keys) + canonical_positions = select_canonical_row_positions(table, variant_keys) + return variant_keys, colliding_base_labels, canonical_positions + + +def disambiguate_labels_if_colliding( + base_label: str, + display_label: str, + row: object, + colliding_base_labels: set[str], +) -> tuple[str, str]: + """Apply `disambiguate_base_label` to (base_label, display_label) when + base_label is a known collision; otherwise return them unchanged. + """ + if base_label not in colliding_base_labels: + return base_label, display_label + base_label = disambiguate_base_label(base_label, row) + nuc_change = str(getattr(row, "variant", "")).strip() + if nuc_change and nuc_change not in {"", "None"}: + display_label = f"{base_label}\n({nuc_change})" + else: + display_label = base_label + return base_label, display_label + + def _prepare_variant_heatmap_matrix( table: pd.DataFrame, sample_names: Sequence[str], @@ -729,11 +1132,20 @@ def _prepare_variant_heatmap_matrix( ) ) + # Computed over the whole, unfiltered table (not just the rows that will + # survive filtering) so label-collision detection and canonical-row + # selection see the complete picture. Positional indices line up 1:1 + # with the `enumerate(table.itertuples(...))` loop that follows, since + # both walk the same unfiltered `table` in the same order. + variant_keys, colliding_base_labels, canonical_positions = compute_variant_labeling( + table + ) + records: List[Dict[str, Union[str, float, int]]] = [] record_index_by_label: dict[str, int] = {} qc_maps: dict[str, dict[str, str]] = {} - for row in table.itertuples(index=False): + for pos, row in enumerate(table.itertuples(index=False)): gene_value = getattr(row, "gene", "") if str(gene_value) in {"5' UTR", "3' UTR", "INTERGENIC"}: continue @@ -748,11 +1160,18 @@ def _prepare_variant_heatmap_matrix( continue gene_label, display_label, base_label = _resolve_variant_labels(row) + base_label, display_label = disambiguate_labels_if_colliding( + base_label, display_label, row, colliding_base_labels + ) change_type = ( str(getattr(row, "type_of_change", "")).strip().lower().lstrip("*") ) - if not include_joint and change_type.startswith("joint"): + if ( + not include_joint + and change_type.startswith("joint") + and pos != canonical_positions.get(variant_keys[pos]) + ): continue if included_patterns and not any( fnmatch.fnmatch(change_type, pattern) for pattern in included_patterns @@ -764,7 +1183,7 @@ def _prepare_variant_heatmap_matrix( if ( only_persistent and str(getattr(row, "persistence_status", "")).strip().lower() - != "new_persistent" + not in NEW_ENDS_PRESENT_STATUSES ): continue if ( @@ -1030,6 +1449,7 @@ def _write_interactive_heatmap_html( cli_command: Optional[str], sample_labels: Sequence[str] | None = None, plot_title: str | None = None, + filename_stem: str = "variant_allele_frequency_heatmap", ) -> None: if matrix.empty: return @@ -1390,7 +1810,7 @@ def _frequency_to_color(value: float) -> tuple[str, str]: """ - output_path = os.path.join(outdir, "variant_allele_frequency_heatmap.html") + output_path = os.path.join(outdir, f"{filename_stem}.html") with open(output_path, "w", encoding="utf-8") as handle: handle.write(html_doc) @@ -1425,6 +1845,7 @@ def generate_variant_heatmap( cli_command: Optional[str] = None, x_tick_labels: Sequence[str] | None = None, plot_title: str | None = None, + filename_stem: str = "variant_allele_frequency_heatmap", ): """Generate a heatmap of variant allele frequencies across passages.""" @@ -1456,7 +1877,7 @@ def generate_variant_heatmap( print("No variant data available for heatmap; skipping plot.") return - heatmap_path = os.path.join(outdir, "variant_allele_frequency_heatmap.pdf") + heatmap_path = os.path.join(outdir, f"{filename_stem}.pdf") fig_width, fig_height = _heatmap_figure_size( len(heatmap_data.index), len(heatmap_data.columns) @@ -1533,6 +1954,7 @@ def generate_variant_heatmap( cli_command, sample_labels=tick_labels, plot_title=heatmap_title, + filename_stem=filename_stem, ) except Exception as html_exc: # pragma: no cover - best-effort UX print(f"Warning: failed to generate interactive heatmap: {html_exc}") diff --git a/vartracker/annotation_processing.py b/vartracker/annotation_processing.py index 7a40543..cef7f5b 100644 --- a/vartracker/annotation_processing.py +++ b/vartracker/annotation_processing.py @@ -55,15 +55,27 @@ def _parse_attrs(attr: str) -> dict: return out -def gene_lengths_from_gff3(gff_path: str | Path) -> Dict[str, int]: +def _parse_gene_cds_by_contig( + gff_path: str | Path, +) -> Tuple[ + Dict[Tuple[str, str], int], + Dict[Tuple[str, str], int], + Tuple[int, int] | None, + int | None, + int | None, +]: """ - Parse a GFF3 and return {gene_name: CDS_length}, plus 5'/3' UTR and INTERGENIC=1. + Parse a GFF3 and return per-(contig, gene) CDS lengths and first-start + positions, plus whole-genome bounds for UTR computation. Robust to feature order (e.g., CDS before mRNA). Prefers: 1) CDS attribute 'gene' 2) Parent=gene: → gene 'Name' 3) Parent=transcript: → mRNA 'Name' or its parent gene 'Name' 4) As last resort, protein_id/tx id (rare). + + Shared by `gene_lengths_from_gff3` and `ambiguous_gene_names` so both draw + on the same single parse of the annotation. """ gff_path = Path(gff_path) @@ -71,10 +83,11 @@ def gene_lengths_from_gff3(gff_path: str | Path) -> Dict[str, int]: gene_id_to_name: Dict[str, str] = {} # gene-id -> gene symbol/name tx_id_to_gene_name: Dict[str, str] = {} # transcript-id -> gene name - # Accumulators - gene_cds_bp: Dict[str, int] = defaultdict(int) - gene_first_start: Dict[str, int] = {} - pending_by_tx: Dict[str, List[tuple[int, int]]] = defaultdict(list) + # Accumulators, keyed by (seqid, gene_name) so identically-named genes on + # different contigs are never conflated. + gene_cds_bp: Dict[Tuple[str, str], int] = defaultdict(int) + gene_first_start: Dict[Tuple[str, str], int] = {} + pending_by_tx: Dict[str, List[Tuple[int, int, str]]] = defaultdict(list) # For UTR computation (whole-genome) seq_range: Tuple[int, int] | None = None @@ -149,34 +162,94 @@ def gene_lengths_from_gff3(gff_path: str | Path) -> Dict[str, int]: gname = tx_id_to_gene_name.get(tid) if not gname: # Defer until after we see the mRNA - pending_by_tx[tid].append((cds_len, start_i)) + pending_by_tx[tid].append((cds_len, start_i, seqid)) continue # will add later once tx->gene is known if not gname: # Skip unresolvable CDS entries; they will be handled later continue - gene_cds_bp[str(gname)] += cds_len - if gname is not None: - current = gene_first_start.get(str(gname)) - gene_first_start[str(gname)] = ( - start_i if current is None else min(current, start_i) - ) + key = (seqid, str(gname)) + gene_cds_bp[key] += cds_len + current = gene_first_start.get(key) + gene_first_start[key] = ( + start_i if current is None else min(current, start_i) + ) # --- Resolve any pending CDS whose transcripts appeared after - for tid, lengths in pending_by_tx.items(): + for tid, entries in pending_by_tx.items(): gname = tx_id_to_gene_name.get(tid) if not gname: # Last resort: use the transcript id as a stand-in (or skip) gname = tid if not gname: continue - gene_cds_bp[str(gname)] += sum(length for length, _ in lengths) - if gname is not None: - current = gene_first_start.get(str(gname)) - first_start = min(start for _, start in lengths) - gene_first_start[str(gname)] = ( - first_start if current is None else min(current, first_start) + # All entries for a given transcript id share the same contig. + seqid = entries[0][2] + key = (seqid, str(gname)) + gene_cds_bp[key] += sum(length for length, _, _ in entries) + current = gene_first_start.get(key) + first_start = min(start for _, start, _ in entries) + gene_first_start[key] = ( + first_start if current is None else min(current, first_start) + ) + + return gene_cds_bp, gene_first_start, seq_range, first_cds_start, last_cds_end + + +def ambiguous_gene_names(gff_path: str | Path) -> set[str]: + """ + Return gene names that occur on more than one contig/replicon in a GFF3. + + Useful for other per-gene aggregations (e.g. results-table gene summaries) + that need to disambiguate the same gene names vartracker's own + `gene_lengths_from_gff3` disambiguates, even for genes with no CDS-length + contribution in a particular downstream table. + """ + gene_cds_bp, _, _, _, _ = _parse_gene_cds_by_contig(gff_path) + gname_to_seqids: Dict[str, set] = defaultdict(set) + for seqid, gname in gene_cds_bp: + gname_to_seqids[gname].add(seqid) + return {gname for gname, seqids in gname_to_seqids.items() if len(seqids) > 1} + + +def gene_lengths_from_gff3(gff_path: str | Path) -> Dict[str, int]: + """ + Parse a GFF3 and return {gene_name: CDS_length}, plus 5'/3' UTR and INTERGENIC=1. + + Gene names are keyed per-contig internally. If the same gene name occurs on + more than one contig/replicon (common for bacterial plasmids reusing names + like 'repA'), the returned keys are disambiguated as " ()" so + unrelated genes are never silently summed together. Gene names that occur + on only one contig are returned unqualified, matching prior behaviour. + """ + ( + gene_cds_bp, + gene_first_start, + seq_range, + first_cds_start, + last_cds_end, + ) = _parse_gene_cds_by_contig(gff_path) + + # --- Collapse (seqid, gene) keys into final gene labels. Gene names that + # only ever appear on a single contig keep their plain name (unchanged + # behaviour for single-contig/viral references and multi-segment genomes + # with unique per-segment gene names). Names shared across contigs are + # disambiguated so their lengths are never summed together. + gname_to_seqids: Dict[str, set] = defaultdict(set) + for seqid, gname in gene_cds_bp: + gname_to_seqids[gname].add(seqid) + + labelled_cds_bp: Dict[str, int] = defaultdict(int) + labelled_first_start: Dict[str, int] = {} + for (seqid, gname), length in gene_cds_bp.items(): + label = gname if len(gname_to_seqids[gname]) == 1 else f"{gname} ({seqid})" + labelled_cds_bp[label] += length + start = gene_first_start.get((seqid, gname)) + if start is not None: + current = labelled_first_start.get(label) + labelled_first_start[label] = ( + start if current is None else min(current, start) ) # --- Add UTRs if we have genome bounds @@ -184,16 +257,17 @@ def gene_lengths_from_gff3(gff_path: str | Path) -> Dict[str, int]: seq_start, seq_end = seq_range five_utr = max(0, first_cds_start - seq_start) # 1..first_cds_start-1 three_utr = max(0, seq_end - last_cds_end) # last_cds_end+1..end - gene_cds_bp["5' UTR"] = five_utr - gene_cds_bp["3' UTR"] = three_utr + labelled_cds_bp["5' UTR"] = five_utr + labelled_cds_bp["3' UTR"] = three_utr # Conventional placeholder - gene_cds_bp["INTERGENIC"] = 1 + labelled_cds_bp["INTERGENIC"] = 1 ordered = { - gene: gene_cds_bp[gene] + gene: labelled_cds_bp[gene] for gene in sorted( - gene_cds_bp.keys(), key=lambda g: gene_first_start.get(g, float("inf")) + labelled_cds_bp.keys(), + key=lambda g: labelled_first_start.get(g, float("inf")), ) } diff --git a/vartracker/main.py b/vartracker/main.py index ddbc2b1..6da2720 100644 --- a/vartracker/main.py +++ b/vartracker/main.py @@ -32,12 +32,19 @@ ProcessingError, FILE_COLUMNS, ) -from .vcf_processing import annotate_vcf, format_vcf, merge_consequences, process_vcf +from .vcf_processing import ( + NEW_ENDS_PRESENT_STATUSES, + annotate_vcf, + format_vcf, + merge_consequences, + process_vcf, +) from .analysis import ( process_joint_variants, generate_cumulative_lineplot, generate_gene_table, plot_gene_table, + parse_plot_genes_arg, generate_variant_heatmap, search_literature, ) @@ -75,6 +82,7 @@ ) from .data import parse_pokay as parse_pokay_module from .annotation_processing import ( + ambiguous_gene_names, gene_lengths_from_gff3, validate_reference_and_annotation, ) @@ -181,7 +189,11 @@ def add_legacy(name: str, dest: str, **kwargs) -> None: action="store_true", default=False, dest="heatmap_include_joint", - help="Include joint variants in heatmaps (default: exclude them)", + help=( + "Also show non-canonical joint/compound annotation-group rows in " + "heatmaps (a variant's canonical row is always shown by default, " + "even if joint)" + ), ) add_legacy("include-joint", "heatmap_include_joint", action="store_true") group.add_argument( @@ -189,7 +201,10 @@ def add_legacy(name: str, dest: str, **kwargs) -> None: action="store_true", default=False, dest="heatmap_only_persistent", - help="Only include variants with persistence_status == new_persistent", + help=( + "Only include new variants present at the final timepoint " + "(persistence_status new_persistent or new_intermittent)" + ), ) add_legacy("only-persistent", "heatmap_only_persistent", action="store_true") group.add_argument( @@ -385,7 +400,10 @@ def _add_shared_plot_filter_arguments(group: argparse._ArgumentGroup) -> None: "--persistent-only", action="store_true", default=False, - help="Only include variants with persistence_status == new_persistent", + help=( + "Only include new variants present at the final timepoint " + "(persistence_status new_persistent or new_intermittent)" + ), ) group.add_argument( "--new-only", @@ -674,6 +692,7 @@ def _configure_vcf_parser( include_input_csv: bool, input_csv_required: bool = False, include_consensus_options: bool = False, + consensus_group: argparse._ArgumentGroup | None = None, ) -> None: if include_input_csv: if input_csv_required: @@ -736,6 +755,24 @@ def _configure_vcf_parser( "one sample after filtering (default: error)" ), ) + analysis_group.add_argument( + "--local-csq", + action="store_true", + default=False, + help=( + "Pass bcftools csq the -l/--local-csq flag, so each variant is " + "annotated on its own rather than jointly with nearby co-occurring " + "variants ('joint_*' consequence types). Off by default, since " + "joint calling is the correct behaviour for genuinely linked " + "variants (e.g. adjacent-codon changes on the same haplotype). " + "Consider enabling this for datasets with many unphased, " + "sub-consensus variants clustered in the same gene, where " + "bcftools has no genotype evidence that co-occurring variants are " + "actually on the same haplotype and joint calling can otherwise " + "fragment a single variant's presence/absence trajectory across " + "samples." + ), + ) analysis_group.add_argument( "-d", "--min-depth", @@ -743,10 +780,17 @@ def _configure_vcf_parser( required=False, type=int, default=10, - help="Minimum depth threshold for variant QC (default: 10)", + help=( + "Minimum depth threshold for variant QC (default: 10). Below " + "this depth at a sample, genuine absence of a variant cannot be " + "distinguished from dropout/non-detection; that sample is " + "flagged 'F' in the per_sample_variant_qc output column rather " + "than treated as confident absence." + ), ) if include_consensus_options: - analysis_group.add_argument( + group = consensus_group or analysis_group + group.add_argument( "--consensus-snp-min-af", action="store", required=False, @@ -757,7 +801,7 @@ def _configure_vcf_parser( "considered for consensus (default: 0.25)" ), ) - analysis_group.add_argument( + group.add_argument( "--consensus-snp-thresh", action="store", required=False, @@ -768,7 +812,7 @@ def _configure_vcf_parser( "consensus base (default: 0.75)" ), ) - analysis_group.add_argument( + group.add_argument( "--consensus-indel-thresh", action="store", required=False, @@ -783,6 +827,31 @@ def _configure_vcf_parser( help="Only analyse samples with sample_number less than or equal to this value", default=None, ) + analysis_group.add_argument( + "--max-plot-genes", + action="store", + type=int, + default=30, + help=( + "Cap the gene-wise summary FIGURE to the top N genes, ranked by " + "number of newly emerged variants (total variants used as a " + "tiebreak). Does not affect the tabular/TSV output, which always " + "includes every annotated gene. Overridden by --plot-genes. " + "(default: 30)" + ), + ) + analysis_group.add_argument( + "--plot-genes", + action="store", + default=None, + help=( + "Comma-separated gene names, or a path to a file with one gene " + "name per line, to show on the gene-wise summary FIGURE. " + "Overrides --max-plot-genes. Genes absent from the reference " + "annotation or with no variants in this dataset are skipped " + "with a warning." + ), + ) analysis_group.add_argument( "--literature-csv", action="store", @@ -1034,6 +1103,7 @@ def _add_bam_subparser(subparsers): include_input_csv=True, input_csv_required=False, include_consensus_options=True, + consensus_group=snk_group, ) _move_action_group_after( bam_parser, "Snakemake options", "Vartracker Analysis Options" @@ -1463,6 +1533,21 @@ def _add_plot_heatmap_subparser(subparsers): _add_heatmap_option_arguments( heatmap_group, prefix="", add_legacy_prefixed_aliases=True ) + output_group = parser.add_argument_group("Output") + output_group.add_argument( + "--out", + default=None, + help=( + "Write the heatmap using this path as the base name (extension, if " + "any, is ignored) - e.g. --out plots/myheatmap writes " + "plots/myheatmap.pdf and plots/myheatmap.html" + ), + ) + output_group.add_argument( + "--outdir", + default=None, + help="Output directory for heatmap files (default: beside results.csv)", + ) parser.set_defaults(handler=_run_plot_heatmap_command) @@ -1645,7 +1730,10 @@ def _add_plot_genome_subparser(subparsers): "--persistent-only", action="store_true", default=False, - help="Only include variants with persistence_status == new_persistent", + help=( + "Only include new variants present at the final timepoint " + "(persistence_status new_persistent or new_intermittent)" + ), ) filter_group.add_argument( "--new-only", @@ -2148,6 +2236,7 @@ def _add_e2e_subparser(subparsers): e2e_parser, include_input_csv=False, include_consensus_options=True, + consensus_group=snk_group, ) _move_action_group_after( e2e_parser, "Snakemake options", "Vartracker Analysis Options" @@ -2222,6 +2311,15 @@ def _run_e2e_command(args): optional_empty={"bam", "vcf", "coverage", "reads2"}, ) + if not args.primer_bed: + print( + "\033[93mWarning:\033[0m no --primer-bed supplied. If this data was " + "generated with amplicon sequencing, primer-binding sites will not " + "be clipped and may produce spurious variant calls near " + "primer-overlap regions. If your library prep is not amplicon-based, " + "this does not apply and can be ignored." + ) + print(get_logo()) rulegraph_path = _normalise_rulegraph_path(args.rulegraph) @@ -2473,7 +2571,15 @@ def _run_plot_heatmap_command(args): if not results_csv.exists(): raise InputValidationError(f"Results CSV not found: {results_csv}") - outdir = results_csv.parent + filename_stem = "variant_allele_frequency_heatmap" + if args.out: + out_path = Path(args.out).expanduser().resolve() + outdir = out_path.parent + filename_stem = out_path.stem or filename_stem + elif args.outdir: + outdir = Path(args.outdir).expanduser().resolve() + else: + outdir = results_csv.parent outdir.mkdir(parents=True, exist_ok=True) table = pd.read_csv(results_csv, keep_default_na=False) @@ -2541,9 +2647,13 @@ def _run_plot_heatmap_command(args): plot_title=args.title, literature_hits=literature_df, literature_table_path=literature_path, + filename_stem=filename_stem, **_collect_heatmap_kwargs(args), ) - print(f"\nFinished: find results in {outdir}\n") + print( + f"\nFinished: wrote {outdir / f'{filename_stem}.pdf'} and " + f"{outdir / f'{filename_stem}.html'}\n" + ) return 0 except (InputValidationError, ProcessingError) as exc: print(f"\nERROR: {exc}\n") @@ -2723,6 +2833,7 @@ def _process_files( args.gff3, args.debug, args.multiallelic_overflow, + local_csq=args.local_csq, ) # Process VCF and extract variants @@ -2782,8 +2893,21 @@ def _process_files( os.path.join(args.outdir, "cumulative_mutations.pdf"), ) - gene_table = generate_gene_table(table, gene_lengths) - plot_gene_table(gene_table, pname, args.outdir) + ambiguous_genes = None + if gene_lengths is not None and getattr(args, "gff3", None): + ambiguous_genes = ambiguous_gene_names(args.gff3) + plot_genes = None + if getattr(args, "plot_genes", None): + plot_genes = parse_plot_genes_arg(args.plot_genes) + + gene_table = generate_gene_table(table, gene_lengths, ambiguous_genes) + plot_gene_table( + gene_table, + pname, + args.outdir, + max_plot_genes=getattr(args, "max_plot_genes", 30), + plot_genes=plot_genes, + ) plot_variant_turnover( *apply_shared_plot_filters( *prepare_plot_inputs(table)[:2], @@ -2822,9 +2946,11 @@ def _process_files( ].reset_index(drop=True) new_mutations.to_csv(os.path.join(args.outdir, "new_mutations.csv"), index=None) - persistent_mutations = table[table.persistence_status == "new_persistent"][ - ["gene", "variant", "amino_acid_consequence", "nsp_aa_change"] - ].reset_index(drop=True) + persistent_mutations = table[ + table.persistence_status.isin(NEW_ENDS_PRESENT_STATUSES) + ][["gene", "variant", "amino_acid_consequence", "nsp_aa_change"]].reset_index( + drop=True + ) persistent_mutations.to_csv( os.path.join(args.outdir, "persistent_new_mutations.csv"), index=None ) diff --git a/vartracker/plotting.py b/vartracker/plotting.py index 0ae53b0..07a9991 100644 --- a/vartracker/plotting.py +++ b/vartracker/plotting.py @@ -20,8 +20,11 @@ _coerce_frequency, _extract_numeric_position, _resolve_variant_labels, + compute_variant_labeling, + disambiguate_labels_if_colliding, ) from .core import InputValidationError, ProcessingError +from .vcf_processing import NEW_ENDS_PRESENT_STATUSES DEFAULT_TRAJECTORY_TOP_N = 12 DEFAULT_LIFESPAN_TOP_N = 20 @@ -113,22 +116,6 @@ def _project_name_from_results(table: pd.DataFrame, explicit_name: str | None) - return names[0] if len(names) == 1 else "" -def _build_variant_key(row: object, fallback_label: str) -> str: - chrom = str(getattr(row, "chrom", "")).strip() - start = str(getattr(row, "start", "")).strip() - ref = str(getattr(row, "ref", "")).strip() - alt = str(getattr(row, "alt", "")).strip() - variant_name = str(getattr(row, "variant", "")).strip() - - if chrom and start and ref and alt: - return f"{chrom}:{start}:{ref}>{alt}" - if chrom and start and variant_name: - return f"{chrom}:{start}:{variant_name}" - if variant_name: - return variant_name - return fallback_label - - def _display_label_prefix(label: object) -> str: return str(label).split(" (", 1)[0].strip() @@ -146,11 +133,21 @@ def prepare_plot_inputs( table: pd.DataFrame, ) -> tuple[pd.DataFrame, pd.DataFrame, list[str], list[int]]: sample_names, sample_numbers = _get_sample_axis(table) + dedup_table = table.drop_duplicates().reset_index(drop=True) + + all_variant_keys, colliding_base_labels, canonical_positions = ( + compute_variant_labeling(dedup_table) + ) + long_records: list[dict[str, object]] = [] - for row in table.drop_duplicates().itertuples(index=False): - _, display_label, base_label = _resolve_variant_labels(row) - variant_key = _build_variant_key(row, base_label) + for position, row in enumerate(dedup_table.itertuples(index=False)): + variant_key = all_variant_keys[position] + canonical_row = dedup_table.iloc[canonical_positions[variant_key]] + _, display_label, base_label = _resolve_variant_labels(canonical_row) + base_label, display_label = disambiguate_labels_if_colliding( + base_label, display_label, canonical_row, colliding_base_labels + ) names = _parse_slash_tokens(getattr(row, "samples", "")) numbers = _parse_slash_tokens(getattr(row, "sample_number", "")) freqs = [ @@ -184,13 +181,19 @@ def prepare_plot_inputs( { "variant_id": variant_key, "variant_label": display_label.replace("\n", " "), - "variant_name": str(getattr(row, "variant", "")).strip(), - "gene": str(getattr(row, "gene", "")).strip(), - "type_of_change": str(getattr(row, "type_of_change", "")).strip(), - "type_of_variant": str(getattr(row, "type_of_variant", "")).strip(), - "variant_status": str(getattr(row, "variant_status", "")).strip(), + "variant_name": str(getattr(canonical_row, "variant", "")).strip(), + "gene": str(getattr(canonical_row, "gene", "")).strip(), + "type_of_change": str( + getattr(canonical_row, "type_of_change", "") + ).strip(), + "type_of_variant": str( + getattr(canonical_row, "type_of_variant", "") + ).strip(), + "variant_status": str( + getattr(canonical_row, "variant_status", "") + ).strip(), "persistence_status": str( - getattr(row, "persistence_status", "") + getattr(canonical_row, "persistence_status", "") ).strip(), "sample_name": sample_name, "sample_number": int(float(sample_number)), @@ -221,6 +224,36 @@ def prepare_plot_inputs( ["sample_number", "variant_id", "sample_name"] ).reset_index(drop=True) + # Collapse duplicate (variant_id, sample_number) rows created when a + # variant has multiple bcftools-csq annotation-group rows (see + # _build_variant_key / process_vcf's row-canonicalisation): a row's + # presence/AF is masked to whichever samples that specific annotation + # group matched, so the same sample can appear as absent on one row and + # present on another for the same variant. Take the union of presence + # and the max AF across a variant's rows for each sample, so plots see + # one point per (variant, sample) instead of drawing/counting duplicates. + # Metadata columns are already uniform per variant_id at this point + # (sourced from each variant's canonical row above), so "first" here is + # well-defined rather than arbitrary. + long_df = long_df.groupby(["variant_id", "sample_number"], as_index=False).agg( + variant_label=("variant_label", "first"), + variant_name=("variant_name", "first"), + gene=("gene", "first"), + type_of_change=("type_of_change", "first"), + type_of_variant=("type_of_variant", "first"), + variant_status=("variant_status", "first"), + persistence_status=("persistence_status", "first"), + sample_name=("sample_name", "first"), + allele_frequency=("allele_frequency", "max"), + present=("present", "max"), + sample_qc=("sample_qc", "first"), + sample_pass_qc=("sample_pass_qc", "first"), + has_literature=("has_literature", "max"), + ) + long_df = long_df.sort_values( + ["sample_number", "variant_id", "sample_name"] + ).reset_index(drop=True) + summary = long_df.groupby("variant_id", as_index=False).agg( variant_label=("variant_label", "first"), variant_name=("variant_name", "first"), @@ -249,7 +282,9 @@ def prepare_plot_inputs( summary["duration"] = summary["last_seen"].fillna(0).astype(float) - summary[ "first_seen" ].fillna(0).astype(float) - summary["is_persistent_new"] = summary["persistence_status"].eq("new_persistent") + summary["is_persistent_new"] = summary["persistence_status"].isin( + NEW_ENDS_PRESENT_STATUSES + ) summary["is_nonsynonymous"] = ( ~summary["type_of_change"].str.lower().str.contains("synonymous", na=False) ) @@ -308,7 +343,7 @@ def apply_shared_plot_filters( if persistent_only: filtered_summary = filtered_summary[ - filtered_summary["persistence_status"].eq("new_persistent") + filtered_summary["persistence_status"].isin(NEW_ENDS_PRESENT_STATUSES) ] if new_only: @@ -1143,10 +1178,16 @@ def load_reference_feature_metadata(results_csv: str | Path) -> dict[str, object def _collapse_variants_for_genome_plot( table: pd.DataFrame, ) -> pd.DataFrame: + dedup_table = table.drop_duplicates().reset_index(drop=True) + + all_variant_keys, colliding_base_labels, canonical_positions = ( + compute_variant_labeling(dedup_table) + ) + records: list[dict[str, object]] = [] - for row in table.drop_duplicates().itertuples(index=False): + for position, row in enumerate(dedup_table.itertuples(index=False)): gene_label, display_label, base_label = _resolve_variant_labels(row) - variant_key = _build_variant_key(row, base_label) + variant_key = all_variant_keys[position] af_values = [ _coerce_frequency(token) for token in _parse_slash_tokens(getattr(row, "alt_freq", "")) @@ -1193,6 +1234,39 @@ def _collapse_variants_for_genome_plot( ) continue first = group.iloc[0].to_dict() + canonical_row = dedup_table.iloc[canonical_positions[variant_id]] + canonical_gene_label, canonical_display_label, canonical_base_label = ( + _resolve_variant_labels(canonical_row) + ) + canonical_base_label, canonical_display_label = ( + disambiguate_labels_if_colliding( + canonical_base_label, + canonical_display_label, + canonical_row, + colliding_base_labels, + ) + ) + first["variant_label"] = canonical_display_label.replace("\n", " ") + first["variant_name"] = str(getattr(canonical_row, "variant", "")).strip() + first["plot_gene"] = str(canonical_gene_label).strip() + first["raw_gene"] = str(getattr(canonical_row, "gene", "")).strip() + first["type_of_change"] = str( + getattr(canonical_row, "type_of_change", "") + ).strip() + first["type_of_variant"] = str( + getattr(canonical_row, "type_of_variant", "") + ).strip() + first["variant_status"] = str( + getattr(canonical_row, "variant_status", "") + ).strip() + first["persistence_status"] = str( + getattr(canonical_row, "persistence_status", "") + ).strip() + first["aa_position"] = _extract_numeric_position( + canonical_base_label.split(":", 1)[1] + if ":" in canonical_base_label + else canonical_base_label + ) merged_af_values: list[float] = [] for value in group["af_values"].tolist(): merged_af_values.extend([float(item) for item in cast(list[float], value)]) @@ -1242,7 +1316,7 @@ def _subset_collapsed_variants( if max_af is not None: subset = subset[subset["summary_af"] <= max_af] if persistent_only: - subset = subset[subset["persistence_status"].eq("new_persistent")] + subset = subset[subset["persistence_status"].isin(NEW_ENDS_PRESENT_STATUSES)] if new_only: subset = subset[subset["variant_status"].eq("new")] if subset.empty: diff --git a/vartracker/schemas.py b/vartracker/schemas.py index 8ea37b2..0c513fe 100644 --- a/vartracker/schemas.py +++ b/vartracker/schemas.py @@ -113,9 +113,18 @@ { "name": "persistence_status", "type": "string", - "description": "Persistence class based on first and last sample presence.", + "description": ( + "Persistence class based on presence at the first and last " + "sample, and whether presence was continuous in between. The " + "'_intermittent' classes mean the variant was present at both " + "the first/last (original) or reappeared by the last sample " + "(new), but was absent from at least one sample in between." + ), "units": "", - "values": "original_retained, original_lost, new_persistent, new_transient", + "values": ( + "original_retained, original_intermittent, original_lost, " + "new_persistent, new_intermittent, new_transient" + ), }, { "name": "presence_absence", @@ -141,7 +150,13 @@ { "name": "all_samples_pass_qc", "type": "boolean", - "description": "True if every sample passes per-sample variant QC.", + "description": ( + "True only if every sample is 'P' in per_sample_variant_qc. False " + "means at least one sample's presence/absence call could not be " + "confidently distinguished from dropout/non-detection - inspect " + "per_sample_variant_qc to see which sample(s) are affected before " + "treating this variant's presence/absence pattern as reliable." + ), "units": "", "values": "true, false", }, @@ -155,7 +170,18 @@ { "name": "per_sample_variant_qc", "type": "string (slash-separated)", - "description": "Per-sample QC flags (P/F) ordered by input.", + "description": ( + "Per-sample QC flags (P/F) ordered by input. A sample fails ('F') " + "when there is no variant-supporting read and site coverage is " + "below --min-depth, i.e. when genuine absence of the variant " + "cannot be distinguished from dropout/non-detection; 'P' means " + "absence (or presence) was confidently called at that sample. " + "E.g. 'P / P / F / P / P / P' identifies the third sample as the " + "one where QC failed. The default heatmap marks 'F' cells " + "visually: an unfilled black-bordered rectangle in the static " + "PDF, and a dark inset ring plus a 'QC=FAIL' hover tooltip in the " + "interactive HTML version." + ), "units": "", "values": "P, F", }, diff --git a/vartracker/vcf_processing.py b/vartracker/vcf_processing.py index 0c5f257..64d2370 100644 --- a/vartracker/vcf_processing.py +++ b/vartracker/vcf_processing.py @@ -11,6 +11,7 @@ import shutil from pathlib import Path from typing import TypedDict +from urllib.parse import unquote import pandas as pd import numpy as np @@ -563,8 +564,18 @@ def annotate_vcf( annotation, debug, multiallelic_overflow="error", + local_csq=False, ): - """Annotate a merged multi-sample VCF with bcftools csq.""" + """Annotate a merged multi-sample VCF with bcftools csq. + + local_csq: pass bcftools csq -l/--local-csq, so each variant is + considered on its own rather than jointly with nearby co-occurring + variants. Off by default - joint calling is correct for genuinely + linked variants, but can fragment a single variant's presence across + samples when co-occurrence is just an artefact of unphased, sub-consensus + calls (see analysis.py::process_joint_variants and the row-splitting + this can otherwise cause in process_vcf below). + """ if multiallelic_overflow not in {"error", "drop-lowest-af", "skip-site"}: raise RuntimeError( "multiallelic_overflow must be one of: error, drop-lowest-af, skip-site" @@ -690,8 +701,9 @@ def annotate_vcf( shutil.copyfile(skipped_vcf, output_file) return + local_flag = " -l" if local_csq else "" cmd = ( - f"bcftools csq -p R -f {reference} -g {annotation} --force " + f"bcftools csq -p R{local_flag} -f {reference} -g {annotation} --force " f"-Ov -o {annotated_only_vcf} {prepared_vcf}" ) @@ -709,12 +721,32 @@ def annotate_vcf( shutil.copyfile(annotated_only_vcf, output_file) -def calculate_variant_site_depths(cov_df, v, samples, min_depth: int): +def _mean_depth_in_range(arr, lo: int, hi: int): + """Mean depth over the 1-based inclusive position range [lo, hi]. + + Positions outside the array (beyond the coverage file's span) or marked + missing (-1, the fill value for positions absent from the coverage file) + are excluded rather than treated as zero coverage. Returns None if no + position in range has recorded depth. + """ + lo = max(lo, 1) + hi = min(hi, arr.shape[0] - 1) + if hi < lo: + return None + window = arr[lo : hi + 1] + valid = window[window >= 0] + if valid.size == 0: + return None + return valid.mean() + + +def calculate_variant_site_depths(depth_arrays, v, samples, min_depth: int): """ Calculate depth metrics for variant sites. Args: - cov_df (pd.DataFrame): Coverage data + depth_arrays (dict[str, np.ndarray]): Per-sample coverage depth, + indexed directly by 1-based genome position (index 0 unused). v: Variant object from cyvcf2 samples (list): List of sample names @@ -729,33 +761,23 @@ def calculate_variant_site_depths(cov_df, v, samples, min_depth: int): end = v.end if v.var_type == "snp": - site_depths = list(cov_df.loc[cov_df["pos"] == start]["depth"]) + site_depths = [] + for sample in samples: + arr = depth_arrays[sample] + if 0 <= start < arr.shape[0] and arr[start] >= 0: + site_depths.append(int(arr[start])) else: - site_depths = ( - cov_df[cov_df["pos"].between(start, end, inclusive="both")] - .groupby("sample")["depth"] - .mean() - .reset_index() - ) - site_depths["order"] = pd.Categorical( - site_depths["sample"], categories=samples, ordered=True - ) - site_depths = [ - int(round(x, 0)) for x in list(site_depths.sort_values("order")["depth"]) - ] - - window_depth = ( - cov_df[cov_df["pos"].between(start - 10, end + 10, inclusive="neither")] - .groupby("sample")["depth"] - .mean() - .reset_index() - ) - window_depth["order"] = pd.Categorical( - window_depth["sample"], categories=samples, ordered=True - ) - window_depth = [ - int(round(x, 0)) for x in list(window_depth.sort_values("order")["depth"]) - ] + site_depths = [] + for sample in samples: + mean_depth = _mean_depth_in_range(depth_arrays[sample], start, end) + if mean_depth is not None: + site_depths.append(int(round(mean_depth, 0))) + + window_depth = [] + for sample in samples: + mean_depth = _mean_depth_in_range(depth_arrays[sample], start - 9, end + 9) + if mean_depth is not None: + window_depth.append(int(round(mean_depth, 0))) try: dp_values = v.format("DP").tolist() @@ -799,27 +821,49 @@ def calculate_variant_site_depths(cov_df, v, samples, min_depth: int): return result +# A "new" variant present at the final timepoint is new_persistent when its +# presence was continuous from first appearance onward, or new_intermittent +# when it disappeared and reappeared along the way. Both are new variants +# that reached the final timepoint, just with a different path there - +# consumers that filter on "did this new variant persist to the end" should +# match this whole set, not new_persistent alone. +NEW_ENDS_PRESENT_STATUSES = frozenset({"new_persistent", "new_intermittent"}) + + def _summarise_sample_trajectory(allele_freqs, samples): """Derive trajectory metadata from per-sample allele frequencies.""" presence_absence = ["N" if x == "." else "Y" for x in allele_freqs] variant_status = "new" if allele_freqs[0] == "." else "original" + if "Y" in presence_absence: + first_present_idx = presence_absence.index("Y") + last_present_idx = rindex(presence_absence, "Y") + # A gap is an absence sandwiched between the first and last presence, + # i.e. the variant disappeared and later reappeared rather than + # persisting continuously - see original_intermittent/new_intermittent + # below. + has_gap = "N" in presence_absence[first_present_idx : last_present_idx + 1] + else: + first_present_idx = None + last_present_idx = None + has_gap = False + if allele_freqs[0] != "." and allele_freqs[-1] == ".": persistent_status = "original_lost" elif allele_freqs[0] != "." and allele_freqs[-1] != ".": - persistent_status = "original_retained" + persistent_status = "original_intermittent" if has_gap else "original_retained" elif allele_freqs[0] == "." and allele_freqs[-1] != ".": - persistent_status = "new_persistent" + persistent_status = "new_intermittent" if has_gap else "new_persistent" elif allele_freqs[0] == "." and allele_freqs[-1] == ".": persistent_status = "new_transient" else: persistent_status = "unknown" first_appearance = ( - samples[presence_absence.index("Y")] if "Y" in presence_absence else "None" + samples[first_present_idx] if first_present_idx is not None else "None" ) last_appearance = ( - samples[rindex(presence_absence, "Y")] if "Y" in presence_absence else "None" + samples[last_present_idx] if last_present_idx is not None else "None" ) return { @@ -1069,22 +1113,36 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): "Number of coverage files does not match number of VCF samples." ) - # Load coverage data - cov_list = [] + # Load coverage data. Coverage files follow the `samtools depth -aa` + # convention (one line per genome position, contiguous from position 1), + # so per-sample depth is stored as a plain numpy array indexed directly + # by 1-based position rather than as rows in a concatenated DataFrame - + # for a bacterial-scale genome the DataFrame form holds tens of millions + # of rows (with per-row "ref"/"sample" strings) in memory simultaneously, + # while the array form needs only a few bytes per position per sample. + depth_arrays: dict[str, "np.ndarray"] = {} total_cov_list = [] - for i, cov_file in enumerate(covs): + for sample, cov_file in zip(samples, covs): try: - cov = pd.read_csv(cov_file, sep="\t", names=["ref", "pos", "depth"]) + cov = pd.read_csv( + cov_file, + sep="\t", + names=["ref", "pos", "depth"], + usecols=["pos", "depth"], + dtype={"pos": "int64", "depth": "int32"}, + ) total_cov = "{:.2f}".format(((sum(cov.depth > 10) / cov.shape[0]) * 100)) total_cov_list.append(total_cov) - cov["sample"] = samples[i] - cov_list.append(cov) + + positions = cov["pos"].to_numpy() + depths = cov["depth"].to_numpy() + arr = np.full(int(positions.max()) + 1, -1, dtype=np.int32) + arr[positions] = depths + depth_arrays[sample] = arr except Exception as e: raise RuntimeError(f"Error reading coverage file {cov_file}: {str(e)}") - cov_df = pd.concat(cov_list).set_index("sample") - # Process variants for v in vcf: info = dict(list(v.INFO)) @@ -1094,7 +1152,7 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): allele_freqs = ["."] * max(1, len(samples)) trajectory = _summarise_sample_trajectory(allele_freqs, samples) - depths_qc = calculate_variant_site_depths(cov_df, v, samples, min_depth) + depths_qc = calculate_variant_site_depths(depth_arrays, v, samples, min_depth) all_samples_pass_qc = "F" not in depths_qc["variant_qc"] proportion_samples_passing_qc = ( sum(flag == "P" for flag in depths_qc["variant_qc"]) @@ -1111,6 +1169,8 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): produced_annotation_specific_row = False if sample_bcsq_map: + trajectories_by_alt: dict = {} + groups_by_alt: dict = {} for annotation_group in annotation_groups: annot = _representative_bcsq_annotation(annotation_group) anno = annot.split("|") @@ -1134,9 +1194,29 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): masked_trajectory = _summarise_sample_trajectory( masked_allele_freqs, samples ) + + if annotation_alt is not None: + groups_by_alt.setdefault(annotation_alt, []).append( + annotation_group + ) + trajectories_by_alt.setdefault(annotation_alt, []).append( + tuple(masked_trajectory["presence_absence"]) + ) + if "Y" not in masked_trajectory["presence_absence"]: continue + # Every joint/compound annotation-group row is kept, not + # just standalone ones: a compound consequence that only + # ever shows up in one sample can still be real biology + # (e.g. a newly emerging linked pair), and recurrence + # across samples can only ever confirm something after + # the fact - it cannot distinguish a genuine first + # occurrence from an artefact, so filtering on it would + # silently discard exactly the most novel findings along + # with the noise. `joint_variant` (set downstream by + # analysis.py::process_joint_variants) flags which rows + # these are, so nothing here is deleted, only labelled. result = _process_annotation( v, anno, @@ -1156,6 +1236,61 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): results.append(result) produced_annotation_specific_row = True + # A gene with many closely-spaced variants can lead bcftools + # to describe every sample's occurrence of this same allele + # as a *different* joint/compound consequence (whichever + # neighbouring variants happened to co-occur in that + # particular sample), with no sample ever getting a + # standalone description. Each annotation-group row above is + # then masked to only the sample(s) matching its own exact + # compound description - e.g. six samples can each land in a + # *different* one-sample-wide row, so the union of all rows + # covers every sample yet no single row ever shows the + # variant's true, continuous trajectory (and persistence + # classification, computed per-row, is wrong for every one + # of them). Whenever no already-emitted row's trajectory for + # this allele exactly matches the true, unmasked per-sample + # presence, synthesise one extra row carrying that true + # trajectory (the same genotype/AF-derived data used for the + # single-record-allele case above) alongside the existing + # joint-detail rows, so the variant's real presence and + # persistence_status is available from at least one row. + for alt, true_allele_freqs in alt_frequency_map.items(): + true_trajectory = _summarise_sample_trajectory( + true_allele_freqs, samples + ) + true_tuple = tuple(true_trajectory["presence_absence"]) + if "Y" not in true_tuple: + continue + if true_tuple in trajectories_by_alt.get(alt, []): + continue + + candidate_groups = groups_by_alt.get(alt) or annotation_groups + representative_group = min( + candidate_groups, + key=lambda g: len(_representative_bcsq_annotation(g)), + ) + annot = _representative_bcsq_annotation(representative_group) + anno = annot.split("|") + result = _process_annotation( + v, + anno, + true_trajectory["variant_status"], + true_trajectory["persistent_status"], + true_trajectory["presence_absence"], + true_trajectory["first_appearance"], + true_trajectory["last_appearance"], + all_samples_pass_qc, + proportion_samples_passing_qc, + depths_qc, + true_allele_freqs, + samples, + total_cov_list, + alt, + ) + results.append(result) + produced_annotation_specific_row = True + if not produced_annotation_specific_row: for annotation_group in annotation_groups: annot = _representative_bcsq_annotation(annotation_group) @@ -1290,10 +1425,11 @@ def _process_annotation( } else: # Regular annotation + gene = unquote(anno[1]) reformatted_aa = ( - reformat_csq_notation(anno[1], anno[5]) + reformat_csq_notation(gene, anno[5]) if len(anno) > 5 - else [anno[1] + ":" + anno[0], ""] + else [gene + ":" + anno[0], ""] ) aa_exploration = AminoAcidChange(reformatted_aa[0].replace("*", "")) @@ -1302,7 +1438,7 @@ def _process_annotation( "chrom": v.CHROM, "start": v.start + 1, "end": v.end, - "gene": anno[1], + "gene": gene, "ref": v.REF, "alt": selected_alt, "variant": v.REF + str(v.POS) + selected_alt,