Metadata-Version: 2.4
Name: vartracker
Version: 2.3.0
Summary: Track the persistence (or loss) of mutations during long-term passaging
Project-URL: Homepage, https://github.com/charlesfoster/vartracker
Project-URL: Bug Reports, https://github.com/charlesfoster/vartracker/issues
Project-URL: Source, https://github.com/charlesfoster/vartracker
Project-URL: Documentation, https://github.com/charlesfoster/vartracker/blob/main/README.md
Author: Dr Charles Foster
Maintainer-email: Dr Charles Foster <charles.foster@unsw.edu.au>
License-Expression: MIT
License-File: LICENSE
Keywords: SARS-CoV-2,VCF,bioinformatics,genomics,longitudinal-analysis,mutation-tracking,passaging,variant-analysis
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Medical Science Apps.
Requires-Python: >=3.11
Requires-Dist: biopython<2,>=1.83
Requires-Dist: cyvcf2>=0.30.0
Requires-Dist: matplotlib>=3.5.0
Requires-Dist: numpy>=1.20.0
Requires-Dist: pandas>=1.3.0
Requires-Dist: seaborn>=0.11.0
Requires-Dist: snakemake>=9
Provides-Extra: dev
Requires-Dist: black; extra == 'dev'
Requires-Dist: flake8; extra == 'dev'
Requires-Dist: mypy<1.20,>=1.19; extra == 'dev'
Requires-Dist: pre-commit>=3.6; extra == 'dev'
Requires-Dist: pytest-cov; extra == 'dev'
Requires-Dist: pytest>=6.0; extra == 'dev'
Provides-Extra: test
Requires-Dist: pytest; extra == 'test'
Requires-Dist: pytest-cov; extra == 'test'
Description-Content-Type: text/markdown

[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
[![PyPI version](https://img.shields.io/pypi/v/vartracker.svg)](https://pypi.org/project/vartracker/)

```
██    ██  █████  ██████  ████████ ██████   █████   ██████ ██   ██ ███████ ██████
██    ██ ██   ██ ██   ██    ██    ██   ██ ██   ██ ██      ██  ██  ██      ██   ██
██    ██ ███████ ██████     ██    ██████  ███████ ██      █████   █████   ██████
 ██  ██  ██   ██ ██   ██    ██    ██   ██ ██   ██ ██      ██  ██  ██      ██   ██
  ████   ██   ██ ██   ██    ██    ██   ██ ██   ██  ██████ ██   ██ ███████ ██   ██


```
# vartracker

A bioinformatics pipeline to summarise variants called against a reference in a longitudinal study design. Written to investigate longitudinal sequencing data from long-term passaging of SARS-CoV-2. However, with appropriate reference data it can be expanded to other pathogens too.

**Author:** Dr Charles Foster

## Table of Contents

- [Features](#features)
- [Installation](#installation)
- [Quick Start](#quick-start)
- [Output](#output)
- [What does vartracker do?](#what-does-vartracker-do)
- [Limitations](#limitations)
- [Citation](#citation)
- [License](#license)
- [Contributing](#contributing)
- [Support](#support)

## Features

- Track mutation persistence across longitudinal samples
- Comprehensive variant analysis including amino acid consequences
- Built-in SARS-CoV-2 reference data and annotations
- Integration with functional mutation databases (literature)
- Automated plotting and statistical analysis
- Support for both SNPs and indels
- Quality control metrics for variants

## Installation

The simplest, and preferred, installation route is via a conda-compatible package manager (pixi, conda, or mamba).

## Conda/Mamba
`vartracker` and its external bioinformatics dependencies can be installed from the package channels directly:

```bash
mamba create -n vartracker -c conda-forge -c bioconda vartracker
mamba activate vartracker
```

If you prefer `conda`:

```bash
conda create -n vartracker -c conda-forge -c bioconda vartracker
conda activate vartracker
```

### Pixi

You can also install `vartracker` via `pixi` either globally or inside an existing workspace:

```bash
# global installation
pixi global install vartracker

# add into an existing workspace
pixi workspace channel add conda-forge
pixi workspace channel add bioconda
pixi add vartracker
```

### Biocontainers

Every Bioconda package is available as a container image for usage with your preferred container runtime. An example command to pull `vartracker` with `docker`:

```bash
# latest build
docker pull quay.io/biocontainers/vartracker:latest
# specific tag
docker pull quay.io/biocontainers/vartracker:<tag>
```

### Alternative: PyPI (Python-only)

If you want a Python-only install (requires Python 3.11 or newer), you can still install from PyPI. In that case you must provide the required external bioinformatics tools yourself (see below):

```bash
pip install vartracker
```

### External Dependencies

vartracker shells out to a handful of bioinformatics tools. Make sure they are discoverable on `PATH` before running the CLI.
Minimum tested versions are tracked in `docs/DEPENDENCIES.md`.

- **bcftools** and **tabix** – required for all modes
- **samtools**, **lofreq**, **fastp**, **bwa**, and **snakemake** – required for the `bam` and `end-to-end` Snakemake workflows

If you only plan to run `vartracker vcf` against pre-generated VCFs, the first pair is sufficient. The additional tools are needed whenever you ask vartracker to align reads or call variants for you.
Consensus genome generation in the `bam` and `end-to-end` workflows uses `bcftools` and `samtools`; it does not require `bedtools`.

Note: the pinned micromamba environment installs `tabix`/`bgzip` via `htslib`.

#### Installing bcftools and tabix

**On macOS:**
```bash
# Using Homebrew
brew install bcftools htslib samtools fastp bwa
# lofreq is available via bioconda (requires conda/mamba)
conda install -c bioconda lofreq

# Using MacPorts
sudo port install bcftools htslib samtools fastp bwa
```

**On Linux (Ubuntu/Debian):**
```bash
sudo apt-get update
sudo apt-get install bcftools tabix samtools fastp bwa
# lofreq is easiest to install via bioconda on Debian-based systems:
conda install -c bioconda lofreq
```

**On Linux (CentOS/RHEL/Fedora):**
```bash
# CentOS/RHEL with EPEL
sudo yum install epel-release
sudo yum install bcftools htslib samtools fastp bwa

# Fedora
sudo dnf install bcftools htslib samtools fastp bwa
# Install lofreq via bioconda on RPM-based systems:
conda install -c bioconda lofreq
```

**Using conda:**
```bash
conda install -c bioconda bcftools samtools tabix fastp bwa lofreq
```

### Development Installation

For development or to get the latest version (requires Python 3.11+):

```bash
git clone https://github.com/charlesfoster/vartracker.git
cd vartracker
pip install -e .[dev]
pre-commit install
```

### Docker: self-build

Build a container image that bundles Python, vartracker, and all external bioinformatics tools:

```bash
  # released version on Bioconda
  docker build -t vartracker:release .
  # development version
  docker build -t vartracker:dev -f Dockerfile.dev .
```

Docker is a self-contained reproducible option. If you publish the image, record the digest and
set it when running to include it in the run manifest:

```bash
export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.3.0
export VARTRACKER_CONTAINER_DIGEST=sha256:...
```

Run workflows by mounting your data directory into the container. The command below analyses an input CSV located in the current directory and writes results beside it:

```bash
docker run --rm -v "$(pwd)":/workspace vartracker \
  vcf /workspace/inputs/vcf_inputs.csv \
  --outdir /workspace/results
```

## Quick Start

After installation, `vartracker` will be available as a command-line tool:

```bash
vartracker --help
```

### Typical commands

```bash
# Analyse pre-called VCFs plus coverage files
vartracker vcf path/to/vcf_inputs.csv --outdir results/vcf_run

# Run BAMs through the Snakemake workflow, then summarise variants
vartracker bam path/to/bam_inputs.csv \
  --snakemake-outdir work/bam_pipeline \
  --outdir results/bam_summary

# Start from raw reads (FASTQ) and run the full pipeline
vartracker end-to-end path/to/read_inputs.csv \
  --cores 12 \
  --outdir results/e2e_summary

# Re-plot a heatmap from an existing vartracker results file
vartracker plot heatmap results/results.csv \
  --aa-exclude "*frameshift*" \
  --x-labels sample-number \
  --literature-csv results/sample.literature_database_hits.full.csv \
  --title "Variant allele frequencies"

# Plot whole-dataset turnover from an existing results file
vartracker plot turnover results/results.csv

# Plot collapsed variant frequencies along the genome
vartracker plot genome results/results.csv

# Zoom to a gene region, optionally using amino-acid coordinates
vartracker plot genome results/results.csv --gene F --aa-scale

# Plot selected variant trajectories from an existing results file
vartracker plot trajectory results/results.csv \
  --variants "S:D614G,S:E484K,S:N501Y"

# Plot takeover-style trajectories using AF thresholds
vartracker plot trajectory results/results.csv \
  --thresholds 0.5,0.9 \
  --crossing-only

# Generate a template spreadsheet for a directory of files
vartracker prepare spreadsheet --mode e2e --dir data/passaging --out inputs.csv

# Build a reference FASTA+GFF3 bundle from GenBank accessions
vartracker prepare reference --accessions CY114381,CY114382 --outdir refs/flu --prefix flu_ref

# Exercise the bundled smoke-test dataset
vartracker vcf --test
vartracker bam --test
vartracker end-to-end --test
```

All modes understand `--test`, which copies the example dataset from `vartracker/test_data`
into a temporary directory, resolves relative paths, and runs the appropriate workflow.

Temporary LoFreq note:
- In `bam` and `end-to-end` mode, `vartracker` currently caps `lofreq call-parallel` at 8 threads even if `--cores` is higher.
- This is a temporary workaround for an older Bioconda LoFreq build that can fail during `call-parallel` final filtering when many shards produce an excessively long merged VCF header.
- The cap will be revisited once an updated LoFreq build is available through Bioconda.

LoFreq primer-overlap rescue:
- Amplicon schemes can create a specific LoFreq false-negative mode: after primer clipping, reads from one strand may be soft clipped at primer-overlap sites, so a genuine near-fixed variant can fail LoFreq's default strand-bias filter.
- `bam` and `end-to-end` therefore run LoFreq with `--no-default-filter`, then apply the normal `lofreq filter` step so standard LoFreq PASS calls are unchanged.
- With the default `--lofreq-primer-rescue auto`, the rescue step runs only when `--primer-bed` is supplied. In other words, `auto` means "use primer rescue when an amplicon primer scheme has been explicitly provided."
- In `end-to-end` mode, the same `--primer-bed` is used for `samtools ampliconclip` and for rescue. In `bam` mode, vartracker does not clip the input BAMs; the primer BED is used only to identify primer-overlap sites for rescue.
- Rescue candidates must be single-ALT SNPs that overlap a primer interval, fail LoFreq's default strand-bias filtering, and pass conservative near-fixed thresholds (`AF>=0.95`, `DP>=100`, `DP4 alt count>=95`, `QUAL>=100`, `DP4 ref count<=20`, minor ALT strand fraction `<=0.05`). Indels, multi-ALT records, lower-frequency variants, non-primer-overlap variants, and variants filtered for non-strand-bias reasons are not rescued by this rule.
- The raw LoFreq calls are retained as `<sample>_variants.raw.vcf.gz` and listed in the updated spreadsheet as `raw_vcf`.
- Rescued variants are marked with `FILTER=RESCUED_PRIMER_OVERLAP`, `INFO/PRIMER_OVERLAP`, and `INFO/RESCUED_BY=overlap_primer_interval`; per-sample details are written to `<sample>_variants.rescued.tsv` and listed in the updated spreadsheet as `lofreq_rescued_tsv`.
- Variants called by raw LoFreq but filtered out of the final VCF are written to `<sample>_variants.filtered_out.tsv` with the LoFreq filter reason and core metrics. This is useful for auditing high-frequency calls that fail strand-bias or other LoFreq filters.
- Use `--lofreq-primer-rescue off` to disable rescue even when a primer BED is supplied, or `--lofreq-primer-rescue on` to require rescue and fail if `--primer-bed` is missing. The rescue thresholds can be adjusted with the `--lofreq-rescue-*` options.

Example amplicon run with primer rescue:

```bash
vartracker end-to-end inputs.csv \
  --primer-bed primers.bed \
  --ampliconclip-tolerance 1 \
  --outdir results/e2e_amplicon
```

### Input Spreadsheets

Every CLI mode reads the same canonical columns:

- `sample_name` (required) – display name for the sample
- `sample_number` (required) – passage/order index used in longitudinal plots
- `reads1`, `reads2` – FASTQ paths (required for `end-to-end`, optional elsewhere). The pipeline runs in single-end mode (leave the `reads2` column empty) but the results are less well tested.
- `bam` – BAM file aligned against the SARS-CoV-2 reference
- `vcf` – bgzipped VCF containing variant calls with depth (`DP`) and allele-frequency tags
- `coverage` – per-base coverage TSV with columns `reference<TAB>position<TAB>depth`

Mode-specific expectations:

- **VCF mode** requires `vcf` and `coverage`, while leaving `reads*`/`bam` empty.
- **BAM mode** requires `bam` and will fill `vcf` + `coverage` during the workflow.
- **End-to-end mode** requires `reads1` (and optionally `reads2`); remaining fields are generated.

The `bam` and `end-to-end` workflows also write two consensus FASTA columns to
the updated Snakemake spreadsheet, plus LoFreq audit columns:
`consensus` for a simple consensus, `iupac_consensus` for an IUPAC-aware
consensus, `raw_vcf` for raw LoFreq calls, `lofreq_rescued_tsv` for the
per-sample primer-overlap rescue table, and `lofreq_filtered_out_tsv` for raw
LoFreq records excluded from the final VCF. SNPs below
`--consensus-snp-min-af` are ignored, SNPs from `--consensus-snp-min-af` up to
`--consensus-snp-thresh` stay as reference bases in the simple consensus
and become REF+ALT ambiguity codes in the IUPAC consensus, and SNPs at or above
`--consensus-snp-thresh` become ALT bases. Indels are controlled
separately by `--consensus-indel-thresh` in both consensus modes. Low-depth bases
are masked as `N`, except for called deletion intervals so true deletions are not
converted to low-depth masks.

Relative paths are resolved with respect to the CSV location, so you can store the sheet alongside
your sequencing artefacts. The `prepare spreadsheet` subcommand can scaffold a CSV and highlight missing files.

Coverage files can be produced with `samtools depth -aa sample.bam > sample_depth.txt` or
`bedtools genomecov -ibam sample.bam -d`. The file name suffix does not matter; vartracker checks
for both `.depth.txt` and `_depth.txt` patterns when preparing its internal test dataset.

### Mode-specific options

- `vartracker vcf` – accepts core analysis options such as `--min-snv-freq`, `--min-indel-freq`,
  `--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`,
  `--consensus-snp-thresh`, and `--consensus-indel-thresh`.
- `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 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 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 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 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`.
- `--qc`: comma-separated `all_samples_pass_qc` patterns to include. Accepted values include `true`, `false`, `pass`, and `fail`.
- `--min-prop-passing-qc`: minimum fraction of samples that must pass per-sample QC.
- `--min-persistence`: minimum number of included samples in which the variant must be present.
- `--min-max-af`: minimum maximum allele frequency across included samples.
- `--min-sample-af`: minimum allele frequency that must be reached in at least one included sample.
- `--sample-subset`: comma-separated sample-name patterns to plot.
- `--hide-singletons`: hide variants present in only one included sample.
- `--min-depth`: minimum site depth a variant must reach in at least one included sample.
- `--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 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.

Standalone plot output:
- `--out`: write to an exact file path.
- `--outdir`: write beside the results CSV or into the chosen directory using deterministic names such as `variant_trajectory_plot.pdf` or `variant_genome_plot.pdf`.
- `--format`: choose `pdf`, `png`, or `svg`.
- `--dpi`: set raster output resolution.

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 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.
- `--include-indels`: opt in to plotting indels too. This may be ambiguous or hard to interpret.
- The standalone genome plot auto-discovers `reference_features.json` beside `results.csv`; workflow runs generate this sidecar automatically.

Trajectory threshold mode:
- `--thresholds`: draw horizontal AF threshold lines, e.g. `0.5,0.9`.
- `--crossing-only`: keep only variants crossing at least one supplied threshold.
- `--label-threshold-crossers`: label only threshold-crossing variants to reduce clutter.
- `--crossing-rule`: choose whether threshold equality counts (`at_or_above`) or requires a strict exceedance (`strictly_above`).

Standalone plot examples:
- `vartracker plot genome results.csv`
- `vartracker plot genome results.csv --gene F`
- `vartracker plot genome results.csv --gene F --aa-scale`
- `vartracker plot genome results.csv --gene F --cds-scale --focus-coords "184-210,586-630"`
- `vartracker plot genome results.csv --focus-coords "150-300,900-1800"`
- `vartracker plot genome results.csv --focus-coords "62-69,196-210;31-42,323-332,379-399;254-277"`
- `vartracker plot genome results.csv --gene F --aa-scale --focus-coords "Ø:62-69,196-210;I:31-42,323-332,379-399;II:254-277"`
- `vartracker plot genome results.csv --gene F --aa-scale --focus-coords "Ø:62-69,196-210;I:31-42,323-332,379-399" --show-intersections`
- `vartracker plot genome results.csv --focus-region-file fusion_regions.json`
- `vartracker plot genome results.csv --gene F --aa-scale --focus-coords "50-120,180-220"`
- `vartracker plot turnover results.csv`
- `vartracker plot trajectory results.csv --variants "S:D614G,S:E484K"`
- `vartracker plot trajectory results.csv --thresholds 0.5,0.9`
- `vartracker plot trajectory results.csv --thresholds 0.5,0.9 --crossing-only`
- `vartracker plot trajectory results.csv --thresholds 0.5,0.9 --crossing-only --label-threshold-crossers`
- `vartracker plot lifespan results.csv --top-n 20 --persistent-only`

Note:
- The standalone `plot` commands require `results.csv` files written by current vartracker versions, which now include a slash-separated `sample_number` column for stable passage ordering.

- `vartracker prepare spreadsheet` – specify `--mode` (`vcf`, `bam`, or `e2e`), `--dir` to scan, `--out` for the CSV,
  and `--dry-run` to preview without writing a file.
- `vartracker prepare reference` – build a merged FASTA/GFF3 bundle from GenBank nucleotide accessions.
  Use `--accessions` or `--accession-file`, plus `--outdir`. Optional flags: `--prefix`, `--force`,
  `--keep-intermediates`, `--skip-csq-validation`.

### Using Literature Database

To search mutations against functional databases:

1. **Set up a literature database (optional):**
```bash
parse_pokay pokay_database.csv
```
   This command automatically downloads the required literature files from the
   pokay repository into `pokay_literature/NC_045512` (override with
   `--download-dir`) and writes the processed CSV for downstream analysis.

2. **Run vartracker with literature search:**
```bash
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
<file>`). 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

```
usage: main.py [-h] [-V] {vcf,bam,end-to-end,e2e,prepare,schema} ...

positional arguments:
  {vcf,bam,end-to-end,e2e,prepare,schema}
    vcf                 Analyse VCF inputs
    bam                 Run the BAM preprocessing workflow
    end-to-end (e2e)    Run the end-to-end workflow (Snakemake + vartracker)
    prepare             Prepare inputs and references for vartracker
    schema              Print schemas for results tables or literature CSV input

options:
  -h, --help            show this help message and exit
  -V, --version         show program's version number and exit
```

Use `vartracker <subcommand> --help` to inspect the full list of mode-specific arguments.

### Prepare reference from accessions

Use this workflow to build a `bcftools csq`-ready reference bundle from nucleotide accessions:

```bash
# Comma-separated accessions
# Example with influenza A segments
vartracker prepare reference \
  --accessions CY114381,CY114382,CY114383,CY114384,CY114385,CY114386,CY114387,CY114388 \
  --outdir refs/influenza_a \
  --prefix influenza_a_ref

# One accession per line in a file
vartracker prepare reference \
  --accession-file accessions.txt \
  --outdir refs/
```

Required external tools:

- `bcftools` for csq smoke validation

Outputs:

- `<outdir>/<prefix>.fa`
- `<outdir>/<prefix>.gff3`
- `<outdir>/<prefix>.fa.fai`
- `<outdir>/prepare_metadata.json`

Validation notes:

- Unless `--skip-csq-validation` is supplied, vartracker writes a dummy coding-region VCF variant
  and runs `bcftools csq` against the generated FASTA/GFF3.
- Validation fails fast if `bcftools csq` exits non-zero or if the output VCF does not contain `BCSQ`.

Troubleshooting:

- Accession fetch failures: verify accession spelling and network access to NCBI efetch.
- SeqID mismatch errors: confirm FASTA headers and GFF3 seqids match exactly.
- csq validation failure: inspect the stderr snippet in the error output and confirm `bcftools`
  version and annotation structure.

### Installation Test

After installation you can verify the workflows using the bundled
demonstration dataset:

```bash
vartracker vcf --test --outdir vartracker_vcf_test_results
vartracker bam --test --outdir vartracker_bam_test_results
vartracker end-to-end --test --outdir vartracker_e2e_test_results
```

Each command copies the example dataset, resolves relative paths, checks for
the required external tools, and writes a self-contained set of results.

## Output

vartracker produces several output files:

- **results.csv**: Comprehensive variant analysis with all metrics
- **results_metadata.json**: Output schema version and results metadata
- **`<sample>_variants.raw.vcf.gz`** (`bam`/`end-to-end`): Raw LoFreq calls before default filtering and primer-overlap rescue
- **`<sample>_variants.rescued.tsv`** (`bam`/`end-to-end`): LoFreq primer-overlap rescue audit table, empty when rescue is disabled or no variants are rescued
- **`<sample>_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 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
- **variant_allele_frequency_heatmap.pdf**: Heatmap of variant allele frequencies across passages
- **literature_database_hits.*.csv**: Functional annotation results (if literature search used)
- **run_metadata.json**: Provenance manifest capturing inputs, tool versions, and run status

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:

```bash
vartracker schema results
```

To write the schema to a file instead, use:

```bash
vartracker schema results --out docs/output_schema.csv
vartracker schema results --out docs/output_schema.json --format json
```

To print the expected literature CSV structure for `--literature-csv`, use:

```bash
vartracker schema literature
```

## What does vartracker do?

The pipeline performs the following analysis:

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:
   - Gene location and amino acid consequences
   - Variant type (SNP/indel) and change type (synonymous/missense/etc.)
   - Persistence across samples (new/original, persistent/transient)
   - Quality control metrics
   - Amino acid property changes
   - Allele frequency dynamics

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.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.

Also cite relevant methods or data sources, for example:

- Foster CSP, et al. Long-term serial passaging of SARS-CoV-2 reveals signatures of convergent evolution. Journal of Virology. 2025;99: e00363-25. doi:10.1128/jvi.00363-25
- Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10. doi:10.1093/gigascience/giab008
- Danecek P, McCarthy SA. BCFtools/csq: haplotype-aware variant consequences. Bioinformatics. 2017;33: 2037–2039. doi:10.1093/bioinformatics/btx100
- Wilm A, Aw PPK, Bertrand D, Yeo GHT, Ong SH, Wong CH, et al. LoFreq: a sequence-quality aware, ultra-sensitive variant caller for uncovering cell-population heterogeneity from high-throughput sequencing datasets. Nucleic Acids Res. 2012;40: 11189–11201. doi:10.1093/nar/gks918
- Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34: i884–i890. doi:10.1093/bioinformatics/bty560
- Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv. 2013 [cited 13 Apr 2021]. Available: https://arxiv.org/abs/1303.3997v2
- Mölder F, Jablonski KP, Letcher B, Hall MB, Tomkins-Tinch CH, Sochat V, et al. Sustainable data analysis with Snakemake. F1000Res. 2021;10: 33. doi:10.12688/f1000research.29032.2

## License

This project is licensed under the MIT License - see the [LICENSE](LICENSE) file for details.

## Contributing

Contributions are welcome! Please feel free to submit a Pull Request.

## Support

If you encounter any issues or have questions:

1. Check the [documentation](https://github.com/charlesfoster/vartracker)
2. Search existing [issues](https://github.com/charlesfoster/vartracker/issues)
3. Create a new issue with detailed information about your problem
