cgDist 🧬

Ultra-fast SNP/indel-level distance calculator for core and whole genome MLST (cg/wgMLST) analysis

cgDist is a high-performance Rust implementation for calculating genetic distances in bacterial genomics, specifically designed for epidemiological outbreak investigations and phylogenetic analysis.

🚀 Features

⚡ Ultra-fast: Parallel processing with optimized algorithms
🎯 Precision: SNP/indel-level distance calculation
🔧 Flexible: Multiple hashing algorithms (CRC32, MD5, SHA256)
📊 Comprehensive: Built-in comparison tools and statistical analysis
📋 Per-pair data quality & dashboard: Optional long-format table of shared loci / missingness per pair, missingness confidence intervals, and a self-contained offline HTML dashboard for outbreak-cluster exploration
🧬 Recombination-candidate flagging: Per-locus mutation-density screen to flag loci as recombination candidates for downstream phylogenetic confirmation
💾 Efficient: LZ4 compression for fast caching
📈 Scalable: Memory-efficient processing of large datasets

📚 Table of Contents

Features
Installation
Quick Start
Usage
Per-Pair Quality Table & Analyst Dashboard
Recombination-Candidate Flagging
Cache Inspector
Custom Hashers Plugin System
API Documentation
Citation
Support
License

🔧 Installation

Prerequisites

Rust 1.88 or later (the minimum supported Rust version, MSRV, is also declared in Cargo.toml). The pinned dependency set requires a recent toolchain, so simply using the latest stable Rust is recommended. The easiest way to install or update Rust is via rustup.rs:
```
# Install rustup (skip if already installed)
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

# If rustup is already installed but Rust is older than 1.88, update with:
rustup update stable
```
Python 3.8+ (only for the validation scripts in validation_test/)
System build dependencies for parasail-rs: the alignment backend is built from C source via CMake, which requires a C compiler and zlib development headers. Install once per machine:
```
# Debian / Ubuntu / WSL
sudo apt install build-essential cmake zlib1g-dev

# RHEL / AlmaLinux / Rocky / CentOS / Fedora
sudo dnf install gcc gcc-c++ cmake zlib-devel

# macOS (Homebrew; Xcode Command Line Tools provide compiler + zlib)
xcode-select --install
brew install cmake
```
Windows users are encouraged to use the Docker image or WSL2, which provide a ready-to-build Linux environment. Native Windows builds additionally require zlib via vcpkg (vcpkg install zlib) or MSYS2.

From Source

# Clone the repository
git clone https://github.com/genpat-it/cgDist.git
cd cgDist

# Build release version
cargo build --release

# The binary will be available at ./target/release/cgdist

Install from crates.io (recommended)

cargo install cgdist

This fetches the latest published release from crates.io, builds it locally with your stable Rust toolchain, and installs the cgdist and recombination_candidate_analyzer binaries to ~/.cargo/bin/ (which should already be on your PATH after a default rustup install). Cache inspection is available via cgdist --inspector. The deprecated recombination_analyzer binary is also installed and forwards every argument to recombination_candidate_analyzer with a deprecation notice — existing scripts continue to work.

To pin a specific published version:

cargo install cgdist --version 0.2.1

For a fully reproducible build that uses exactly the dependency versions we tested against, add --locked (this reuses the published Cargo.lock instead of re-resolving to the newest compatible dependencies, which may otherwise require an even more recent Rust toolchain):

cargo install cgdist --locked

Install from GitHub (specific tag or unreleased commits)

To install directly from the GitHub repository — useful for installing an unreleased commit or for fully self-contained reproducibility when citing the manuscript:

# Specific release tag
cargo install --git https://github.com/genpat-it/cgDist --tag v0.2.1 cgdist

# Latest state on the default branch
cargo install --git https://github.com/genpat-it/cgDist cgdist

cgdist is a binary crate, so its Cargo.lock is committed to the repository to guarantee reproducible builds — this is the convention recommended in the official Cargo FAQ for binary crates.

Docker

A multi-arch (linux/amd64 + linux/arm64) image is published to GitHub Container Registry on every release:

# Pull the public image (no authentication required)
docker pull ghcr.io/genpat-it/cgdist:0.2.1
# or pin to the minor / major series:
# docker pull ghcr.io/genpat-it/cgdist:0.2
# docker pull ghcr.io/genpat-it/cgdist:latest   # tracks master HEAD

# Run with the image (mount your working directory at /data).
# Options are passed to cgdist:
docker run --rm -v $(pwd):/data ghcr.io/genpat-it/cgdist:0.2.1 \
    --schema /data/schema_dir --profiles /data/profiles.tsv \
    --output /data/distances.tsv --mode snps-indel-bases

# The other programs of the image run by name:
docker run --rm -v $(pwd):/data ghcr.io/genpat-it/cgdist:0.2.1 \
    cgdist-cache pull --from <store URL> --out /data/store

The image also works in Nextflow and other workflow managers, which run a shell in the container (/bin/bash ...), with no extra options. It runs as a non-root user: with Docker, set docker.runOptions = '-u $(id -u):$(id -g)' in nextflow.config so that it can write to the work directory.

To build the image locally instead of pulling (useful for development):

docker build -t cgdist:dev .
docker run --rm -v $(pwd):/data cgdist:dev --help

🚀 Quick Start

Basic Distance Calculation

# Calculate SNP distances from cgMLST profiles
cgdist --schema schema_dir/ --profiles profiles.tsv --output distances.tsv

# Use different distance mode
cgdist --schema schema_dir/ --profiles profiles.tsv --output distances.tsv --mode snps-indel-bases

# Use different hashing algorithm
cgdist --schema schema_dir/ --profiles profiles.tsv --output distances.tsv --hasher-type sha256

# Enable cache for faster recomputation
cgdist --schema schema_dir/ --profiles profiles.tsv --output distances.tsv --cache-file cache.lz4

# Specify number of threads
cgdist --schema schema_dir/ --profiles profiles.tsv --output distances.tsv --threads 16

Validation / Smoke Test

A self-contained validation suite with a small embedded test dataset (3 loci, 10 samples, ~3 KB) is provided in validation_test/. It verifies algorithmic correctness across all four distance modes (Hamming, SNPs, SNPs+InDel-events, SNPs+InDel-bases), checks the mathematical invariant cgDist ≥ Hamming, and confirms Parasail alignment integration.

# After installing or building cgDist (see Installation)
cd validation_test
pip install -r requirements.txt        # one-time: installs pandas

# run_validation.py is self-contained: it locates the cgdist binary
# ($CGDIST_BIN, ../target/release, or PATH), regenerates the four distance
# matrices, and validates them.
python3 run_validation.py

Expected output: 🎉 ALL VALIDATION TESTS PASSED! See validation_test/README.md for details on the test design, expected distances, and how to regenerate the fixture from scratch.

The validation suite also runs automatically in CI on every push and pull request (see .github/workflows/ci-and-docker.yml).

Configuration File

A configuration file is optional: every parameter accepted by cgdist also has a CLI flag. The configuration file simply lets you persist commonly-used settings without retyping them. CLI flags always override TOML values when both are provided.

A canonical example is shipped at examples/cgdist-config.toml; a Hamming-mode variant is at examples/hamming-config.toml. Both files use the same flat key structure (no [sections]), and the same key names as the corresponding CLI flags (the only normalization is that CLI flag dashes become underscores in TOML — e.g. --hasher-type becomes hasher_type).

You can also generate a fresh annotated sample with:

cgdist --generate-config > cgdist-config.toml

A minimal example (alignment-based mode):

profiles = "profiles.tsv"
schema = "schema/"
output = "distances.tsv"
hasher_type = "crc32"
mode = "snps"            # legacy alias snps-indel-events == snps-indel-contiguous (deprecated)
format = "tsv"
missing_char = "-"
threads = 1              # default; set to 0 for auto-detect
hamming_fallback = false # opt-in (see Hamming Fallback section below)

# Use a configuration file
cgdist --config cgdist-config.toml

# CLI overrides example: config says threads=1, but the CLI wins → 16 threads used
cgdist --config cgdist-config.toml --threads 16

CLI vs TOML precedence

When both the configuration file and the command line specify the same parameter, the command-line value wins. Internally this is implemented by loading the TOML first, then overlaying any CLI flag that the user explicitly set. The same rule applies to switches: e.g. if the TOML says hamming_fallback = false but you pass --hamming-fallback on the command line, the fallback will be enabled for that run.

📊 Usage

Command Line Options

cgdist [OPTIONS]

MAIN OPTIONS:
    --schema <PATH>            Path to FASTA schema directory or schema file
    --profiles <PATH>          Path to allelic profile matrix (.tsv or .csv)
    --output <FILE>            Output distance matrix file
    --mode <MODE>              Distance mode [default: snps]
                               Options: snps, snps-indel-contiguous, snps-indel-bases, hamming
                               (legacy alias snps-indel-events == snps-indel-contiguous, deprecated)
    --format <FORMAT>          Output format [default: tsv]
                               Options: tsv, csv, phylip, nexus

FILTERING OPTIONS:
    --min-loci <N>             Minimum shared loci for distance calculation [default: 0]
    --sample-threshold <VAL>   Sample quality filter (0.0-1.0) [default: 0.0]
    --locus-threshold <VAL>    Locus quality filter (0.0-1.0) [default: 0.0]
    --include-samples <REGEX>  Include only samples matching regex pattern
    --exclude-samples <REGEX>  Exclude samples matching regex pattern
    --include-loci <REGEX>     Include only loci matching regex pattern
    --exclude-loci <REGEX>     Exclude loci matching regex pattern
    --include-samples-list <FILE>  Include samples from file (one per line)
    --exclude-samples-list <FILE>  Exclude samples from file (one per line)
    --include-loci-list <FILE>     Include loci from file (one per line)
    --exclude-loci-list <FILE>     Exclude loci from file (one per line)

ALIGNMENT OPTIONS:
    --alignment-mode <MODE>    Alignment mode [default: dna]
                               Options: dna, dna-strict, dna-permissive, custom
    --match-score <N>          Custom match score (enables custom mode)
    --mismatch-penalty <N>     Custom mismatch penalty (enables custom mode)
    --gap-open <N>             Custom gap open penalty (enables custom mode)
    --gap-extend <N>           Custom gap extend penalty (enables custom mode)
    --save-alignments <FILE>   Save detailed alignments to TSV file
    --save-cigar <FILE>        Save one compact CIGAR row per aligned pair (TSV)

CODING / CUSTOM DISTANCE OPTIONS:
    --mode nonsyn-snps         Count only nonsynonymous SNPs
    --mode custom --weights S  Per-locus weighted sum, e.g. "nonsyn=1,indel_events=1"
                               keys: allele, snps, indel_events, indel_bases,
                               syn, nonsyn, frame_disrupted (integer weights)
    --translation-table <N>    NCBI genetic code for syn/nonsyn [default: 11]
    --no-first-codon-as-met    Do not read GTG/TTG/... at codon 1 as Met
    --coding-stats             Store syn/nonsyn counts in the cache for later use

PROTEIN-LEVEL OPTIONS:
    --mode aa-hamming          Loci whose alleles encode different proteins
    --mode aa-substitutions    Amino-acid substitutions
    --mode aa-substitutions-indel-events    ... + InDel events
    --mode aa-substitutions-indel-residues  ... + inserted/deleted residues
    --aa-matrix <NAME|FILE>    blosum30..blosum100, pam10..pam500, or a file [default: blosum62]
    --aa-gap-open <N>          Gap of L residues costs open+(L-1)*extend [default: 11]
    --aa-gap-extend <N>        [default: 1]
    --protein-cache-file <F>   Cache of protein-pair results (.lz4)
    --verify-alignments <F>    Re-check fraction F (0-1) of new alignments against
                               parasail's original kernel; stop on any difference

PERFORMANCE OPTIONS:
    --threads <N>              Number of threads [default: 1; pass 0 for auto-detect]
    --cache-file <FILE>        Cache file path (.lz4 extension)
    --cache-note <TEXT>        Note to save with cache
    --cache-only               Build cache only without computing distance matrix
    --cache-dir <DIR>          Cache store (one compact file per locus; read + written)
    --cache-layer <SRC>        Read-only cache store: directory, .cgpack or URL (repeatable)
    --force-recompute          Force recomputation ignoring cache
    --hasher-type <TYPE>       Allele hasher type [default: crc32]
                               Options: crc32, sha256, md5, sequence, hamming

REPORTING OPTIONS (opt-in; the distance matrix is unaffected):
    --emit-pairs <FILE>        Also write a long-format per-pair quality table
    --report-ci                Add missingness confidence-interval columns to --emit-pairs
    --ci-level <VAL>           Confidence level for --report-ci, in (0, 1) [default: 0.95]
    --report <FILE.html>       Also write a self-contained HTML analyst dashboard
    --recomb-threshold <PCT>   Per-locus mutation density (%) for the dashboard's
                               recombination view (needs enriched cache) [default: 3.0]

CACHE ENRICHMENT OPTIONS:
    --enrich-lengths           Enrich cache with nucleotide sequence lengths from schema
    --enrich-output <FILE>     Output file for enriched cache [default: overwrites input cache]

OTHER OPTIONS:
    --missing-char <CHAR>      Missing data character [default: -]
    --no-hamming-fallback      Disable Hamming fallback for SNPs mode
    --stats-only               Show matrix statistics only
    --benchmark                Measure alignment processing speed
    --benchmark-duration <N>   Benchmark duration in seconds [default: 15]
    --dry-run                  Validate inputs without computation
    --inspector <FILE>         Inspect cache file
    --config <FILE>            Path to TOML configuration file
    --generate-config          Generate sample configuration file
    --help                     Display usage information

Supported Input Formats

Schema (FASTA directory):

Individual FASTA files per locus
Each file contains allele sequences
File names correspond to locus names

Profiles (allelic profiles):

TSV: Tab-separated values
CSV: Comma-separated values
Format: Sample name | Locus1 | Locus2 | ... | LocusN
Missing data represented by configurable character (default: -)

Cache files:

LZ4: Compressed cache files (.lz4 or .bin extension)
Automatic compression/decompression

Output Formats

TSV: Tab-separated distance matrix (default)
CSV: Comma-separated distance matrix
PHYLIP: Phylogenetic analysis format
NEXUS: Nexus format for phylogenetic tools

Optional extra outputs (see Per-Pair Quality Table & Analyst Dashboard): --emit-pairs (long-format per-pair TSV) and --report (self-contained HTML).

📋 Per-Pair Quality Table & Analyst Dashboard

A distance between two samples is only as reliable as the number of loci it rests on. cgDist can report this alongside the matrix. Both outputs are opt-in and computed from the same (cached) per-locus distances, so they add little run time and never change the distance matrix.

# Distance matrix + per-pair quality table with 95% missingness CI
# + offline HTML dashboard
cgdist --schema schema_dir --profiles profiles.tsv --output distances.tsv \
    --mode snps --cache-file cache.lz4 \
    --emit-pairs pairs.tsv --report-ci --report report.html

Per-pair table (`--emit-pairs`)

One row per unordered pair (upper triangle), tab-separated, after three # header lines (command, timestamp, version):

Column	Meaning
`sample_i`, `sample_j`	Sample IDs
`distance`	Same value as in the matrix; `NA` if the pair fails `--min-loci`
`shared_loci`	Loci called in both samples (the distance is computed over these)
`total_loci`	Loci in the analysis after filtering
`missing_frac`	`1 - shared_loci / total_loci`

With --report-ci four columns are added:

Column	Meaning
`dist_norm`	Estimated full-schema distance: observed distance + expected contribution of the unshared loci
`ci_low`, `ci_high`	Missingness confidence interval at `--ci-level` (default 0.95)
`ci_reliable`	`false` in the low-information regime (< 10 differing shared loci and `missing_frac` > 0.15), where no method can guarantee coverage

The interval accounts only for the loci a pair does not share: it answers "how large could this distance be if the missing loci had been called?" It combines an exact Beta-Binomial interval on the number of additional differing loci with a normal compound interval, and was validated against a Monte-Carlo posterior predictive. When shared_loci == total_loci the interval collapses to the observed distance.

HTML dashboard (`--report`)

A single self-contained HTML file (no external assets; works offline and can be e-mailed or archived with the analysis) with:

dataset summary and per-sample data quality (mean shared loci, nearest neighbour);
the pairwise-distance distribution and a clusters-vs-threshold curve;
an interactive single-linkage clustering explorer at a chosen threshold, showing the missingness CI of each linking edge;
a recombination view (per-pair and per-sample recombinant-loci load, top pairs) when the cache is enriched with sequence lengths (--enrich-lengths); a locus counts as recombinant when its mutation density exceeds --recomb-threshold percent (default 3).

Memory note: both outputs hold one record per pair in memory (n·(n−1)/2 pairs), so for very large datasets (tens of thousands of samples) prefer the distance matrix alone or run them on a subset.

📍 SNP and InDel Positions (`--save-cigar`)

The cache stores only counts (SNPs, InDel events, InDel bases) per allele pair. To get where the differences are, --save-cigar <file> writes one row per pair aligned in the run:

locus  hash1  hash2  cigar  snps  indel_events  indel_bases  alignment_score
wgMLST-00047350  250050033  2240303468  152=1X215=1X26=...  10  0  0  2742.00

The CIGAR uses parasail's extended alphabet. hash1's allele is the query and hash2's the reference:

op	meaning	advances
`=`	same base	query and reference
`X`	SNP	query and reference
`I`	base in the query, gap in the reference	query
`D`	base in the reference, gap in the query	reference

To locate each difference, walk the CIGAR while counting query and reference positions. In the example above, the first SNP is at query and reference position 153 (1-based). As in parasail, bases are compared case-insensitively. The file is about 35x smaller than --save-alignments, which stores both gapped sequences. Both files contain only the pairs aligned in that run: pairs already in the cache are not re-aligned, so use --force-recompute (or no cache) to get every pair.

Inspecting one allele pair (`cgdist-diff`)

cgdist-diff takes two allele hashes, as found in hashed profiles and in the cache keys, and lists every difference with its position and bases. It aligns the two alleles exactly as cgdist does:

cgdist-diff --schema schema_dir/ --locus cgMLST-00081096 \
    --hash1 4106171463 --hash2 1998779011 --cache-file cache.lz4

allele 1   hash 1998779011  cgMLST-00081096_9  735 bp
allele 2   hash 4106171463  cgMLST-00081096_29  738 bp
counts     snps=36  indel_events=1  indel_bases=3  score=1353
cigar      23=1X101=1X...2=3D3=1X5=1X
SNP        allele1     24 G  ->  allele2     24 A
...
INS #1     after allele1    725  allele2 726-728  +AAT (3 bp)
SNP        allele1    729 A  ->  allele2    732 C
cache      snps=36  indel_events=1  indel_bases=3  -> MATCH

Allele 1 is the one with the smaller hash, as in the cache. INS/DEL are relative to allele 1.

Each difference is also annotated at protein level. The genetic code is configurable: --translation-table <N> accepts any NCBI table (default 11, Bacterial), and --no-first-codon-as-met stops an alternative start codon (GTG, TTG, ...) from being read as Met at position 1. Effects use Sequence Ontology terms, the vocabulary of snpEff, Ensembl VEP and bcftools csq:

term	meaning
`synonymous_variant`	same amino acid
`stop_retained_variant`	stop codon changed to another stop
`start_retained_variant`	initiator codon changed to another start codon
`missense_variant`	different amino acid
`stop_gained` / `stop_lost`	premature stop / stop codon lost
`start_lost`	initiator codon no longer a start codon
`coding_sequence_variant`	SNP in a codon shifted or split by an InDel (not comparable codon-to-codon)
`frameshift_variant`, `inframe_insertion`, `inframe_deletion`	InDels

SNPs are classified by the change of the whole codon in both alleles, so several SNPs in one codon are evaluated together (haplotype-aware, as in bcftools csq). Protein changes are numbered on allele 1 (p.Glu243Asp, p.Arg8=). The summary counts synonymous and nonsynonymous SNPs and says whether the two proteins are identical.

The two proteins are also aligned, with --aa-matrix, --aa-gap-open and --aa-gap-extend as in the aa-* distance modes. The summary gives their substitutions and InDels, and --protein-diffs lists each amino-acid difference (AA_SUB, AA_INS, AA_DEL) with protein positions, for example AA_SUB protein1 243 E -> protein2 244 D p.Glu243Asp.

🧪 Synonymous / Nonsynonymous and Custom Distances

The SNPs of every aligned pair can be classified at codon level, with the same rules as cgdist-diff, into synonymous (syn), nonsynonymous (nonsyn: missense, stop gained or lost, start lost) and frame-disrupted (frame_disrupted: in codons shifted or split by an InDel). The three always add up to snps. Distance modes that use them:

# only nonsynonymous SNPs
cgdist ... --mode nonsyn-snps --cache-file cache.lz4

# any weighted combination, per locus
cgdist ... --mode custom --weights "nonsyn=1,frame_disrupted=1,indel_events=1"

With --mode custom, the contribution of a locus whose alleles differ is the sum of weight × count over allele (1 per differing locus), snps, indel_events, indel_bases, syn, nonsyn and frame_disrupted. The built-in modes are particular cases:

mode	weights
`hamming`	`allele=1`
`snps`	`snps=1`
`snps-indel-contiguous`	`snps=1,indel_events=1`
`snps-indel-bases`	`snps=1,indel_bases=1`
`nonsyn-snps`	`nonsyn=1`

The counts are stored in the cache as optional fields (syn, nonsyn, frame_disrupted), together with the genetic code they were computed with (--translation-table, --no-first-codon-as-met). Pairs cached without them, or with another genetic code, are realigned when needed. The realignment must reproduce the cached SNP/InDel counts, otherwise cgdist stops with an error. Older cgdist versions read such caches unchanged. --coding-stats stores the counts even when the mode does not use them.

🧫 Protein-Level Distances

Alleles are translated on the fly from the DNA schema: the genetic code is set with --translation-table, and the terminal stop is removed. Each allele gets a protein hash, so alleles that differ only by synonymous mutations share a protein and need no alignment. Distinct proteins are aligned globally with parasail (--aa-matrix, --aa-gap-open, --aa-gap-extend). Amino-acid substitutions and InDels are counted as for DNA:

cgdist ... --mode aa-substitutions-indel-events --protein-cache-file prot.lz4

On the L. monocytogenes test set (300 samples), 26,926 DNA alleles encode 13,453 distinct proteins, and 234,439 DNA allele pairs reduce to 67,763 protein pairs.

Protein results are kept in a separate cache (--protein-cache-file). Its metadata records the genetic code, the matrix and the gap penalties, and cgdist refuses the cache if they differ from the run's. In --mode custom the keys aa_allele, aa_subs, aa_indel_events and aa_indel_residues can be combined with the DNA keys, e.g. --weights "aa_subs=1,aa_indel_events=1,syn=0".

📦 Cache Stores (`cgdist-cache`)

A cache store is a directory with one compact binary file per locus plus a manifest.json. It holds alignment parameters, hasher, genetic code, checksums and counts, and can be precomputed, shared and updated incrementally. The per-locus format stores each allele hash once and encodes pair statistics column-wise before LZ4. It takes about 1-2 bytes per allele pair, against about 38 in a .lz4 cache.

# precompute every allele pair of a schema (incremental and resumable)
cgdist-cache build --schema schema_dir/ --out lm_store --threads 64

# use it: pairs are read from the store, only missing ones are aligned
cgdist ... --cache-layer lm_store            # read-only layer (dir, .cgpack or URL)
cgdist ... --cache-dir my_store              # local store, receives new pairs

# publish as a single file (e.g. Zenodo) and download only what is needed
cgdist-cache pack --store lm_store --out lm_dna.cgpack
cgdist-cache pull --from https://host/lm_dna.cgpack --out ~/.cache/cgdist/lm --profiles profiles.tsv

cgdist-cache info   --store lm_dna.cgpack    # parameters, loci, pairs, size
cgdist-cache verify --store lm_dna.cgpack    # checksums and decoding of every locus
cgdist-cache verify --store lm_dna.cgpack --schema schema_dir/ --realign 0.01
                                             # + recompute 1% of the pairs from the schema
cgdist-cache stats  --store lm_dna.cgpack --out stats.json   # per-locus summaries, histograms
cgdist-cache import --cache cache.lz4 --out lm_store   # convert a cgdist cache

build aligns with cgdist's own engine, so a store holds exactly what cgdist would compute. It only aligns the pairs a locus is missing, which makes an update after new alleles cheap, and saves the manifest after every locus. --verify-alignments re-checks alignments against parasail's original kernel, and --coding-stats also stores the synonymous/nonsynonymous counts. On the L. monocytogenes schema (1,748 loci, 388,648 alleles), all 54,387,751 pairs took 169 s on 64 threads, for a 56 MB store.
pull accepts a directory, a .cgpack file, or a URL of either. For a .cgpack URL it downloads only the manifest and the needed loci (HTTP Range); a server without Range support is read once in full. Every locus is checked against its sha256, and loci already present are skipped, so repeating a pull updates the store.
A store is used only if its hasher and alignment parameters match the run; coding counts are used only if the genetic code matches.

The format is specified in docs/STORE_FORMAT.md; scripts/cgds_reader.py is an independent reference reader (Python standard library only) and scripts/test_store_format.py checks it against cgdist.

New alleles

Pairs involving alleles that are not in a store (e.g. alleles inferred by chewBBACA in your samples) are aligned from --schema and kept in --cache-dir. Give cgdist the schema folder used for allele calling: if a pair cannot be aligned (an allele of the profiles is missing from the schema) cgdist stops with an error and writes no distances, since such a pair would count as 0 and undercount the distance (the pairs aligned so far are still saved to the cache). --allow-unaligned writes the distances anyway, as cgdist 0.1.x did. Locus files record a digest of every allele sequence, so an allele with the same CRC32 as a different store allele is never taken from the store. cgdist-cache pull merges downloaded loci with local pairs, and cgdist-cache build --schema <your schema> --out <pulled store> adds only the pairs of the new alleles.

Protein stores

The same format holds protein-level results for the aa-* modes and aa_* weights. The allele table lists the distinct proteins of each locus (protein hash and amino-acid length) and the pair columns hold amino-acid substitutions, InDel events and InDel residues. The manifest records the genetic code, substitution matrix and gap penalties, which must match the run.

cgdist-cache build --protein --schema schema_dir/ --out lm_prot --threads 64 \
    [--aa-matrix blosum62 --aa-gap-open 11 --aa-gap-extend 1 --translation-table 11]
cgdist ... --mode aa-substitutions --protein-cache-layer lm_prot     # read-only
cgdist ... --mode aa-substitutions --protein-cache-dir my_prot       # read/write

Proteins are translated from the run's schema, so novel alleles are translated too and only protein pairs missing from the stores are aligned. The protein length stored for each hash must equal the translated length, otherwise the store is refused (protein-hash collision or different sequences). For the whole L. monocytogenes schema, the 17,821,184 pairs of distinct proteins took 116 s on 24 threads and 14 MB. A DNA store is never accepted as a protein store, and a protein store never as a DNA store.

🧬 Recombination-Candidate Flagging

cgDist includes a companion screen that flags candidate recombinant loci based on per-locus mutation density. This is not a recombination detector: confirmation of recombination requires downstream phylogeny-aware tools (e.g. Gubbins, ClonalFrameML, fastGEAR). The flagging output identifies which loci warrant that follow-up.

Features

Mutation Density Analysis: Flags loci with high SNP/indel density per alignment as recombination candidates
Hamming Distance Filtering: Focuses analysis on genetically related sample pairs
Pairwise Flagging Summary: Per sample-pair count of flagged loci
EFSA Loci Support: Compatible with standardized loci sets for food safety applications
Distance Matrix Correction: Recomputes distances excluding flagged loci

How to flag candidates

Recombination-candidate flagging is performed by the recombination_candidate_analyzer binary, which reads an enriched cache (a cgDist cache that also stores per-allele sequence lengths) and flags loci whose per-locus nucleotide-difference density exceeds a threshold. This is the workflow described in the paper (Supplementary §S6).

# Step 1: build an enriched cache alongside the distance matrix.
# On a freshly created cache, --enrich-lengths records the sequence lengths
# in place, so a single cgdist run is enough (requires cgdist >= 0.1.4; older
# versions dropped the lengths in this case - see CHANGELOG 0.1.4).
cgdist --schema schema_dir/ --profiles profiles.tsv --output distances.tsv \
    --mode snps-indel-bases --cache-file cache.bin --enrich-lengths

# Step 2: flag candidate loci and write a corrected distance matrix.
recombination_candidate_analyzer \
    --enriched-cache cache.bin \
    --profiles profiles.tsv \
    --distance-matrix distances.tsv \
    --output-matrix corrected_distances.tsv \
    --candidate-recombination-log candidate_recombination_loci.tsv \
    --threshold 3.0          # mutation-density % (default: 3.0; e.g. 5.0 for a stricter screen)

Note — the binary recombination_analyzer and the flag --recombination-log are kept as deprecated aliases for backward compatibility (a deprecation notice is printed when invoked). Existing scripts continue to work.

Inputs

cgDist consumes standard cgMLST outputs. Profiles and schemas can be generated, for example, with ChewBBACA (Silva et al. 2018) or downloaded from Chewie-NS (Mamede et al. 2020).

Enriched cache (.bin): a cgDist cache built with --cache-file + --enrich-lengths
Allelic profiles (TSV/CSV): sample-locus-allele matrix
Distance matrix (TSV): the original distance matrix from cgdist
EFSA loci (optional): TSV file listing loci of interest

Outputs

Corrected distance matrix (--output-matrix): distances recomputed with the flagged candidate loci excluded
Candidate flagging log (--candidate-recombination-log, TSV): one row per flagged sample-pair/locus, with SNP/InDel counts, Hamming distance, average allele length, mutation-density %, and recombination-excess %

Parameters

--threshold: mutation-density percentage above which a locus is flagged (default: 3.0)
--candidate-recombination-log: output flagging log path. Legacy alias --recombination-log is also accepted.
--output-matrix: corrected distance-matrix output path

Complete Workflow Example

# Step 1: build an enriched cache + distance matrix
cgdist --schema schema/ --profiles samples.tsv --output distances.tsv \
    --mode snps-indel-bases --cache-file cache.bin --enrich-lengths

# Step 2: flag candidate loci and write a corrected distance matrix
recombination_candidate_analyzer \
    --enriched-cache cache.bin \
    --profiles samples.tsv \
    --distance-matrix distances.tsv \
    --output-matrix corrected_distances.tsv \
    --candidate-recombination-log candidate_recombination_loci.tsv \
    --threshold 3.0

Interpretation Guidelines

High SNP Density: > 3% flags a locus as a recombination candidate (confirm with phylogeny-aware tools)
High Indel Events: May indicate mobile genetic elements; warrants downstream inspection
Pairwise Patterns: Multiple flagged loci between the same sample pair suggests related strains
Hamming Filtering: Ensures focus on epidemiologically relevant comparisons

Performance Considerations

Memory Usage: ~4-8GB for typical bacterial datasets (1000+ samples)
Processing Time: 2-5 minutes for 21M cache entries on modern hardware
Scalability: Linear with cache size, efficient for large epidemiological studies

Scientific Applications

Outbreak Investigation: Flag candidate recombination loci in transmission chains for downstream confirmation
Evolutionary Analysis: Identify candidate horizontal gene transfer events
Food Safety: Screen for recombination signatures in foodborne pathogens
Antimicrobial Resistance: Flag candidate resistance gene transfer events
Population Genomics: Identify loci that may bias clonal-frame distance estimates

🔍 Cache Inspector

Inspect a cgDist cache — version, hasher type, distance mode, alignment parameters, per-locus entry counts, and any saved note — with the built-in inspector flag:

cgdist --inspector cache.lz4

This reads the cache format written by --cache-file (LZ4-compressed JSON), including length-enriched caches built with --enrich-lengths (the file extension, e.g. .lz4 or .bin, does not matter — the format is the same).

Use it to:

Validate a cache file's integrity before reuse
Audit which alignment parameters and distance mode a cache was built with
Inspect cache size and per-locus entry distribution
Troubleshoot cache compatibility issues

🔌 Custom Hashers Plugin System

cgDist provides a powerful plugin architecture for implementing custom hashing algorithms. This is particularly useful for specialized applications or compatibility with other tools.

Implementing a Custom Hasher

Create a new hasher by implementing the AlleleHasher trait:

use cgdist::hashers::{AlleleHasher, AlleleHash};

/// Example: Simple nucleotide composition hasher
#[derive(Debug)]
pub struct CompositionHasher;

impl AlleleHasher for CompositionHasher {
    fn hash_sequence(&self, sequence: &str) -> AlleleHash {
        // Count nucleotides: A, T, G, C
        let mut counts = [0u8; 4]; // A, T, G, C
        for nucleotide in sequence.chars() {
            match nucleotide.to_ascii_uppercase() {
                'A' => counts[0] += 1,
                'T' => counts[1] += 1,
                'G' => counts[2] += 1,
                'C' => counts[3] += 1,
                _ => {} // Ignore ambiguous bases
            }
        }
        
        // Create hash from composition: AAAAATTTTGGGGCCCC format
        let hash_string = format!("A{}T{}G{}C{}", 
            counts[0], counts[1], counts[2], counts[3]);
        AlleleHash::String(hash_string)
    }
    
    fn parse_allele(&self, allele_str: &str, missing_char: &str) -> Result<AlleleHash, String> {
        if allele_str == missing_char {
            Ok(AlleleHash::Missing)
        } else {
            // Parse composition string or return as-is
            Ok(AlleleHash::String(allele_str.to_string()))
        }
    }
    
    fn name(&self) -> &'static str {
        "composition"
    }
    
    fn description(&self) -> &'static str {
        "Nucleotide composition-based hasher (A/T/G/C counts)"
    }
    
    fn validate_sequence(&self, sequence: &str) -> Result<(), String> {
        // Only allow ATGC nucleotides
        for ch in sequence.chars() {
            match ch.to_ascii_uppercase() {
                'A' | 'T' | 'G' | 'C' | 'N' => {}
                _ => return Err(format!("Invalid nucleotide: {}", ch)),
            }
        }
        Ok(())
    }
}

Registering Your Custom Hasher

use cgdist::hashers::HasherRegistry;

fn main() {
    let mut registry = HasherRegistry::new();
    
    // Register your custom hasher
    registry.register_hasher("composition", Box::new(CompositionHasher));
    
    // Use it like any built-in hasher
    let hasher = registry.get_hasher("composition").unwrap();
    let hash = hasher.hash_sequence("ATCGATCG");
    
    println!("Hash: {}", hash); // Output: A2T2G2C2
}

Advanced Custom Hasher Examples

1. K-mer Based Hasher

#[derive(Debug)]
pub struct KmerHasher {
    k: usize,
}

impl KmerHasher {
    pub fn new(k: usize) -> Self {
        Self { k }
    }
}

impl AlleleHasher for KmerHasher {
    fn hash_sequence(&self, sequence: &str) -> AlleleHash {
        let mut kmers = Vec::new();
        let seq_bytes = sequence.as_bytes();
        
        if seq_bytes.len() >= self.k {
            for i in 0..=(seq_bytes.len() - self.k) {
                let kmer = std::str::from_utf8(&seq_bytes[i..i + self.k])
                    .unwrap_or("")
                    .to_string();
                kmers.push(kmer);
            }
        }
        
        kmers.sort();
        let hash_string = kmers.join("|");
        AlleleHash::String(hash_string)
    }
    
    // ... implement other required methods
}

2. Custom Numeric Hasher

#[derive(Debug)]
pub struct CustomNumericHasher;

impl AlleleHasher for CustomNumericHasher {
    fn hash_sequence(&self, sequence: &str) -> AlleleHash {
        // Convert sequence to custom numeric representation
        let mut hash_value = 0u32;
        for (i, nucleotide) in sequence.chars().enumerate() {
            let base_value = match nucleotide.to_ascii_uppercase() {
                'A' => 0,
                'T' => 1,
                'G' => 2,
                'C' => 3,
                _ => 0, // Default for ambiguous
            };
            // Simple polynomial rolling hash
            hash_value = hash_value.wrapping_mul(4).wrapping_add(base_value);
        }
        AlleleHash::Crc32(hash_value)
    }
    
    fn parse_allele(&self, allele_str: &str, missing_char: &str) -> Result<AlleleHash, String> {
        if allele_str == missing_char {
            Ok(AlleleHash::Missing)
        } else {
            match allele_str.parse::<u32>() {
                Ok(value) => Ok(AlleleHash::Crc32(value)),
                Err(_) => Err(format!("Invalid numeric allele: {}", allele_str)),
            }
        }
    }
    
    fn name(&self) -> &'static str {
        "custom-numeric"
    }
    
    fn description(&self) -> &'static str {
        "Custom polynomial rolling hash for sequences"
    }
}

Integration with cgdist CLI

To use custom hashers with the cgdist command-line tool, you can:

Fork and modify: Add your hasher to the registry in src/main.rs
Configuration file: Load hashers from a configuration file
Dynamic loading: Use Rust's plugin system (advanced)

Example integration in main.rs:

fn create_registry() -> HasherRegistry {
    let mut registry = HasherRegistry::new();
    
    // Add your custom hashers here
    registry.register_hasher("composition", Box::new(CompositionHasher));
    registry.register_hasher("kmer3", Box::new(KmerHasher::new(3)));
    registry.register_hasher("custom-numeric", Box::new(CustomNumericHasher));
    
    registry
}

Use Cases for Custom Hashers

Legacy Compatibility: Match existing tool formats
Domain-Specific: Specialized algorithms for specific organisms
Research: Experimental hashing strategies
Performance: Optimized for specific hardware or datasets
Compliance: Meet specific regulatory or institutional requirements

Best Practices

Deterministic: Ensure same sequence always produces same hash
Collision-Resistant: Minimize hash collisions for your use case
Performance: Consider computational overhead
Validation: Implement robust input validation
Documentation: Provide clear usage examples and limitations

The plugin architecture makes cgDist highly extensible while maintaining backward compatibility with existing workflows.

Running the Custom Hasher Example

See the complete working example:

# Run the custom hasher demonstration
cargo run --example custom_hasher

# Output shows different hashers applied to test sequences:
# 🔌 cgDist Custom Hasher Examples
# ===================================
# 
# 📊 Available Hashers:
#   • crc32: Fast CRC32 checksum (chewBBACA compatible)
#   • composition: Nucleotide composition-based hasher (A/T/G/C counts)
#   • kmer3: K-mer composition hasher (sorted k-mers)
#   • polynomial: Polynomial rolling hash for sequences
# 
# 🧬 Testing hasher: composition
#    Description: Nucleotide composition-based hasher (A/T/G/C counts)
#    ATCGATCGATCG → A3T3G3C3
#    AAATTTGGGCCC → A3T3G3C3
#    ATGCATGCATGC → A3T3G3C3

This example demonstrates practical implementation patterns for:

Composition-based hashing: Count nucleotide frequencies
K-mer analysis: Extract and sort sequence k-mers
Polynomial hashing: Mathematical sequence encoding
Error handling: Validation and missing data management

📖 API Documentation

Rust API

use cgdist::{DistanceCalculator, Config};

// Create calculator with custom config
let config = Config::new()
    .hasher("crc32")
    .threads(8)
    .cache_enabled(true);

let calculator = DistanceCalculator::new(config);

// Calculate distances
let distances = calculator.calculate_from_file("sequences.fasta")?;

Python Integration

import subprocess
import pandas as pd

# Run cgdist from Python
result = subprocess.run([
    'cgdist',
    '--schema', 'schema_dir/',
    '--profiles', 'profiles.tsv',
    '--output', 'distances.tsv',
    '--mode', 'snps-indel-bases'
], capture_output=True, text=True)

# Check for errors
if result.returncode != 0:
    print(f"Error: {result.stderr}")
else:
    # Load results
    distances = pd.read_csv('distances.tsv', sep='\t', index_col=0)
    print(f"Distance matrix shape: {distances.shape}")
    print(distances.head())

📜 Citation

If you use cgDist in your research, please cite:

de Ruvo, A.; Castelli, P.; Bucciacchio, A.; Mangone, I.; Mixão, V.; Borges, V.; Flammini, M.; Radomski, N.; Di Pasquale, A. (2026). cgDist: Nucleotide-level distance calculation from cgMLST allelic profiles. NAR Genomics and Bioinformatics, 8(3), lqag090. DOI: 10.1093/nargab/lqag090

@article{deruvo2026cgdist,
  title = {cgDist: Nucleotide-level distance calculation from cgMLST allelic profiles},
  author = {de Ruvo, Andrea and Castelli, Pierluigi and Bucciacchio, Andrea and Mangone, Iolanda and Mixão, Verónica and Borges, Vítor and Flammini, Michele and Radomski, Nicolas and Di Pasquale, Adriano},
  year = {2026},
  journal = {NAR Genomics and Bioinformatics},
  volume = {8},
  number = {3},
  pages = {lqag090},
  doi = {10.1093/nargab/lqag090},
  note = {Software: https://github.com/genpat-it/cgDist}
}

The earlier preprint remains available on bioRxiv (10.1101/2025.10.16.682749).

📞 Support

Issues: GitHub Issues
Discussions: GitHub Discussions
Email: a.deruvo@izs.it

📄 License

This project is licensed under the MIT License - see the LICENSE file for details.

Made with ❤️ for the bioinformatics community

Name		Name	Last commit message	Last commit date
.githooks		.githooks
.github/workflows		.github/workflows
docs		docs
examples		examples
scripts		scripts
src		src
validation_test		validation_test
.gitignore		.gitignore
CHANGELOG.md		CHANGELOG.md
CITATION.cff		CITATION.cff
Cargo.lock		Cargo.lock
Cargo.toml		Cargo.toml
Dockerfile		Dockerfile
LICENSE		LICENSE
README.md		README.md
docker-entrypoint.sh		docker-entrypoint.sh

Latest commit

History

Folders and files

Repository files navigation