██    ██  █████  ██████  ████████ ██████   █████   ██████ ██   ██ ███████ ██████
██    ██ ██   ██ ██   ██    ██    ██   ██ ██   ██ ██      ██  ██  ██      ██   ██
██    ██ ███████ ██████     ██    ██████  ███████ ██      █████   █████   ██████
 ██  ██  ██   ██ ██   ██    ██    ██   ██ ██   ██ ██      ██  ██  ██      ██   ██
  ████   ██   ██ ██   ██    ██    ██   ██ ██   ██  ██████ ██   ██ ███████ ██   ██

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
• Installation
• Quick Start
• Output
• What does vartracker do?
• Limitations
• Citation
• License
• Contributing
• 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:

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

If you prefer conda:

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:

# 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:

# 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):

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:

# 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):

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):

# 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:

conda install -c bioconda bcftools samtools tabix fastp bwa lofreq

Development Installation

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

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:

  # 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:

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:

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:

vartracker --help

Typical commands

# 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:

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); 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). 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).
• --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) 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):

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:

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:

# 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:

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)
• 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). --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): 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), 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:

vartracker schema results

To write the schema to a file instead, use:

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:

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

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 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
2. Search existing issues
3. Create a new issue with detailed information about your problem