vdjtools — immune-repertoire analysis

TCR/BCR immune-repertoire analysis — a clean-room Python + C++ rewrite of the legacy Groovy/Java vdjtools, standardised on the AIRR schema and polars DataFrames with minimal object-orientation.

Built on the antigenomics ecosystem: seqtree (fuzzy search / e-value engine), vdjmatch (overlap + TCRnet), arda (AIRR annotation + markup repair).

Status: v3.0.0 — the native V(D)J model engine plus the full analytics suite (diversity, overlap/TCRnet, preprocessing, biomarkers, single-cell), longitudinal clonotype dynamics (paired expansion testing + the VDJtrack recapture model), CDR features, and legacy-format ingestion (MiXcr, MiGec, immunoSEQ, IMGT/HighV-QUEST, Vidjil, RTCR, TRUST4, arda). Clonotype columns follow the AIRR junction convention (junction_nt / junction_aa). The legacy v1.x tool lives on the legacy-1.x branch and its releases remain available under the repository tags (v0.0.1 … 1.2.1).

Install

pip install vdjtools

Prebuilt wheels ship for CPython 3.10–3.13 on Linux, macOS (Apple Silicon), and Windows; the native _core C++ extension is bundled (the source distribution compiles it on install).

That one command gives you everything vdjtools advertises — no extras to opt into. The three antigenomics engines it delegates to are base dependencies:

So preprocess.correct(), overlap.tcrnet(), biomarker.metaclonotypes() and model.reference.load_germline() all just work from a plain install — and downstream libraries can depend on plain vdjtools and rely on them being there. All three are imported lazily, so import vdjtools stays light; between them they add only requests on top of what vdjtools already needs. arda fetches its IMGT reference once, on first use.

Two things are still optional, because each has a working alternative or is a side integration:

pip install "vdjtools[overlap]"   # scikit-learn — only for cluster_samples(method="mds");
                                  # method="hclust" works out of the box (scipy)
pip install "vdjtools[sc]"        # single-cell extras: anndata (scverse bridge), pyyaml

MMseqs2 is needed only for arda's alignment/annotation path (model.stitch.annotate) — never for germline lookup, Pgen, generation, or the analytics. Install it via conda/brew if you need it.

Development

Uses uv — one repo-local .venv, no conda:

uv venv && source .venv/bin/activate
uv pip install -e ".[dev,test]"       # builds the _core C++ extension (scikit-build-core)

Or run the bootstrap script (portable across bash/zsh, uv-first with a python -m venv fallback):

bash setup.sh --dev-parents --tests   # or: zsh setup.sh

You need a C++ toolchain (Xcode CLT on macOS, build-essential on Linux) for the native _core extension. MMseqs2 is arda's aligner, needed only for the annotation path and the slow arda round-trip tests — brew install mmseqs2, or use the optional environment.yml conda env which bundles it.

Quickstart — recombination model engine

Precomputed models for all 7 human loci ship in the wheel — no OLGA or download needed:

from vdjtools.model import load_bundled, native
from vdjtools.model.generate import generate

model = load_bundled("TRB", source="olga")     # or source="learned" (fit to real repertoires)

native.pgen_nt(model, "TGTGCCAGCAGC...")        # nucleotide generation probability (native C++)
native.pgen_aa(model, "CASSLAPGATNEKLFF")       # amino-acid Pgen (codon-marginalised)
native.pgen_aa(model, "CASSLAPGATNEKLFF", mismatches=1)   # + the whole Hamming-1 ball
native.pgen_aa_batch(model, seqs, mismatches=1, threads=0)  # Pgen over many CDR3s, thread-parallel (~11×)
generate(model, 1000, seed=1)                    # sample a repertoire -> polars DataFrame
                                                # seed= is process-stable from 3.3.0

Where Pgen sums over recombination scenarios, model.viterbi takes the argmax — the single most likely one, which is the V/D/J boundary markup:

from vdjtools.model import best_scenario

sc = best_scenario(model, "TGTGCCAGCAGCTTAGGGACAGGGGGCTACGAGCAGTACTTC",
                   v="TRBV19*01", j="TRBJ2-7*01")     # ALLELE names, as the model's tables are

sc.v_end, sc.d_call, sc.d_start, sc.d_end, sc.j_start   # 0-based, half-open, in CDR3-nt space
# -> 8, 'TRBD1*01', 15, 24, 26

It reuses the same tables and the same loops as pgen_nt, so the chosen D obeys P(D|J) — a TRBD2–TRBJ1 pair is genomically impossible and cannot be called.

infer_nt goes the other way, reconstructing a nucleotide CDR3 from an amino-acid one — the VDJdb case, where a record carries (V, J, CDR3aa) and no nucleotides:

from vdjtools.model import infer_nt

sc = infer_nt(model, "CASSLGQAYEQYF", v="TRBV5-1*01", j="TRBJ2-3*01")
sc.cdr3_nt, sc.pgen, sc.margin        # sequence, its exact Pgen, and how far ahead of the runner-up

Germline positions are pinned to their segment; each free N-region position takes the nucleotide the insertion model prefers. It reproduces the exponential brute-force oracle exactly on every record the oracle can resolve (25/25 TRG, 19/19 TRA) — fixing the germline trim first and then picking the best codon per residue only manages 9/25 and 4/19, because a trim chosen before the codons pins a codon the true optimum would have trimmed away.

The search is native (the same Pi_L·Pi_R transfer matrix as pgen_aa, with max for the sums): 2.5 ms per human TRB CDR3, 0.5 ms per TRA — all 80k VDJdb records in about 3 minutes. v=/j= take one allele, several (a list or the comma-separated string an ambiguous v_call carries), or nothing at all, in which case the DP marginalizes over every gene at essentially no extra cost.

Matches OLGA's Pgen to machine precision across all 7 loci, and adds tandem-D (D-D) support that OLGA/IGoR lack. Learn a model from your own non-functional reads (out-of-frame or stop-codon — both escaped selection, which is all a generative model needs) with model.infer.infer_native.

Explore any model's recombination Bayes net interactively (entropy, mutual information, marginals):

pip install "vdjtools[examples]"
marimo edit examples/model_explorer.py

Interactive marimo notebooks (data auto-loads from HuggingFace, or a local ~/hf/ copy):

• examples/vaccination_tracking.py — clonotype tracking + the recapture model across yellow-fever / influenza / TBE vaccination time courses (vdjtools.dynamics).
• examples/aging.py — cohort-streaming diversity, clone-size and spectratype vs age.
• examples/ankspond_motif.py — the ankylosing-spondylitis TRBV9 "AS27" motif: disease vs HLA-B27 carriage.
• examples/biomarker_explorer.py — Emerson public-TCR association + co-occurrence.

Command line

pip install vdjtools installs the vdjtools command — the model engine (OLGA/IGoR-style) and the repertoire analytics (over sample files or a metadata table, like the legacy tool):

# recombination model engine — built-in models for all 7 loci (no download)
vdjtools models                                # list the bundled models
vdjtools generate -m TRB -n 1000 -o gen.tsv    # sample sequences   (cf. olga-generate_sequences)
vdjtools pgen seqs.tsv -m TRB -o pgen.tsv      # Pgen per CDR3       (cf. olga-compute_pgen)
vdjtools pgen seqs.tsv -m TRB --mismatches 1   # + the Hamming-1 ball; --v-col/--j-col to condition

# model workshop — a model is a directory, or LOCUS[:source[:organism]]
vdjtools model check TRB:learned                     # audit vs its germline; exits 1 on an error
vdjtools model template --locus TRB -o tmpl/         # scaffold from arda, or your own --germline-v/-j
vdjtools model learn clones.tsv -t tmpl/ -o fitted/  # EM on your sequences (--init template = fine-tune)
vdjtools model log fitted/                           # log-likelihood per iteration
vdjtools model diversity TRB:olga                    # entropy + total diversity estimate
vdjtools model compare TRB:olga TRB:learned --by gene --dot diff.pdf
vdjtools model loglik seqs.tsv TRB:learned           # log-likelihood, free parameters, AIC, BIC
vdjtools model extend fitted/ --locus TRB -o bigger/ # add a larger allele library
vdjtools model export TRB:olga --long -o marginals.tsv

# data — convert any format to the canonical table (TSV, or Parquet by extension), preprocess
vdjtools convert mixcr.txt.gz -o clones.parquet   # MiXcr/immunoSEQ/AIRR/… → canonical Parquet
vdjtools downsample clones.parquet 100000 -o ds.tsv
vdjtools filter clones.parquet --coding --min-freq 1e-4 -o coding.tsv
vdjtools pool s1.tsv s2.tsv s3.tsv --join --min-samples 2 -o joint.tsv

# repertoire analytics — sample files, or a cohort via -m/--metadata + --base-dir
vdjtools diversity      sampleA.tsv sampleB.tsv -o diversity.tsv
vdjtools overlap        *.tsv -o overlap.tsv
vdjtools segment-usage  *.tsv --segment v -o usage.tsv
vdjtools spectratype    *.tsv -o spectra.tsv
vdjtools diversity      -m metadata.txt --base-dir samples/ --threads 8 -o div.tsv   # parallel cohort
vdjtools spectratype    --cohort cohort_parquet/ -o spectra.tsv                       # one streamed pass

# the portable signature — one fixed, named, positional feature vector per sample
vdjtools signature      --preset classify -m metadata.txt --base-dir samples/ -o sig.tsv
vdjtools signature      --preset compact *.tsv -t 0 -o vsig.parquet      # -t 0 = every core
vdjtools presets                                          # the named feature sets, ranked
vdjtools presets classify                                 # what one preset is, and when to use it
vdjtools signature      --describe --preset classify      # the column dictionary; reads no input

# longitudinal — paired within-donor expansion test between two timepoints
vdjtools dynamics day0.tsv day15.tsv -o tracked.tsv

Native vdjtools, AIRR Rearrangement, Parquet, and third-party inputs are auto-detected; every command writes to -o — TSV, or Parquet when the path ends in .parquet / .pq — or to stdout (so it pipes). Cohort commands parallelise over samples with -t/--threads or stream a pre-ingested Parquet cohort with --cohort. Run vdjtools <command> --help for options.

Analytics (Python API)

Every reader returns one canonical polars clonotype frame (AIRR junction columns), and every analysis function takes and returns such frames — so results chain together and drop straight into plotting. A tour of the analysis modules (full runnable walkthrough in the User guide):

from vdjtools import io as vio, stats, features, overlap, preprocess

# load (auto-detects MiXcr / immunoSEQ / AIRR / native / … and converts), or a whole cohort:
sample = vio.read("clones.tsv")
cohort = vio.read_samples(vio.read_metadata("metadata.txt"), base_dir="samples/")

# diversity, rarefaction, segment usage, spectratype
stats.diversity_stats(sample)                 # observed, Chao1, Shannon, inverse-Simpson, d50, …
stats.inext(sample, q=(0, 1, 2))              # Hill-number rarefaction/extrapolation + bootstrap CIs
stats.segment_usage(sample, "v")              # V (or "j") usage;  stats.spectratype(sample)

# CDR3 physicochemistry & k-mers
features.physchem_profile(sample, region="all")

# repertoire overlap & TCRnet (fuzzy/similarity/TCRnet via the [overlap] engine)
overlap.overlap_metrics(sampleA, sampleB)     # F / D / Jaccard / Morisita-Horn …
overlap.tcrnet(sample)                         # per-clonotype neighbourhood enrichment

# preprocessing: downsample to a common depth, error-correct, filter, pool
preprocess.downsample(sample, 100_000)
preprocess.correct(preprocess.filter_functional(sample))

# cross-batch V/J-usage bias: batch-correct usage, then resample the clonotype table
usage = preprocess.correct_vj_usage(cohort, batch_col="batch", transform="sigmoid")  # Vlasova 2026
fixed = preprocess.apply_vj_correction(sampleA, usage, sample_id="A0")

The portable signature — one repertoire in, a fixed named positional feature vector out, on a scale a downstream model can consume without fitting a scaler of its own. This is the statistics half (vsig); the geometry half (rsig, features of the prototype-sum embedding) is mirpy's mir.signature, and the two concatenate on sample_id into one contract:

from vdjtools.signature import vsig, vsig_cohort, columns, describe

v = vsig({"TRB": sample}, tier="standard")    # {column: value}, in frozen layout order
describe("standard")                          # the column dictionary

Every feature carries a variance-stabilising transform chosen from its support — Haldane–Anscombe logit for a proportion, Anscombe arcsine for a share, CLR (k−1 parts) for a composition, log for a count — so that a read count, an isotype fraction and a principal component can sit in one matrix. core ⊂ standard ⊂ full are exact index subsets of one frozen column order. A locus that was not sequenced is nan plus a mask: column, never a zero.

Longitudinal tracking — which clonotypes changed between two timepoints, and the VDJtrack recapture model (Pavlova, Zvyagin & Shugay 2024):

from vdjtools import dynamics

# paired within-donor test: emergent / expanded / persistent / contracted / vanishing
tracked = dynamics.test_pair(day0, day15)                  # depth handled per-pair (effective N)
grouped = dynamics.test_metaclonotypes(day0, day15, scope="1,0,0,1")  # 1-Hamming CDR3 ball first
called  = dynamics.expansion_test(day0, day15)             # edgeR NB-exact caller (log2FC + p)

# VDJtrack size-bucket recapture model — recapture fraction per clone-size class (Beta bands);
# split by a group column + capture_test() for the group effect (see examples/vaccination_tracking.py)
rates = dynamics.capture_rates(pre, post)

Incidence-based clonotype association (Emerson 2017 / Howie 2015 / De Witt 2018 / Vlasova 2026) — a choice of test, condition, and co-occurrence — and single-cell paired-chain Pgen:

from vdjtools import biomarker, sc
from vdjtools.biomarker import association, condition

# feature vs condition: Fisher / chi2 / Bayesian / permutation; binary, per-HLA-allele, or CMH-stratified
association(cohort, condition.binary(meta, "cmv"), test=["fisher", "bayes_bf"])
association(cohort, condition.stratified(meta, "cmv", "hla"), stratum_col="_stratum")  # CMV | HLA (CMH)

# feature vs feature: in-silico α-β pairing / same-chain co-specificity (θ lift + Fisher + FDR)
biomarker.cooccurrence(cohort, chain_a="TRA", chain_b="TRB", evalue=True)

sc.paired_pgen(sc.pair_chains(sc.read_10x("filtered_contig_annotations.csv")))  # pgen_alpha·pgen_beta

Performance

The Pgen / generation / EM / diversity hot paths are a native C++ (pybind11) core; everything else is polars. Amino-acid Pgen matches OLGA to machine precision (1e-15) across all 7 loci while being several times faster, and the built-in models keep the resident set small. Single thread, Apple M3 (arm64), bundled human TRB model:

Nucleotide Pgen (via the same transfer-matrix DP as the aa path — an in-frame CDR3 is an aa query with one codon fixed per position) is exact vs OLGA across all loci. Batched Pgen / 1-mismatch over many CDR3s parallelises over sequences (native.pgen_aa_batch, ~11× on 16 cores, bitwise-identical to the serial result); the EM E-step parallelises over reads (~6.7× on 8 threads); diversity/rarefaction run on a native iNEXT kernel (bootstrap + parallel batch). Memory stays light — ~63 MB resident for import vdjtools plus one loaded model, ~123 MB with all seven bundled models resident. Reproduce with ~/vcs/projects/2026-vdjtools-benchmark/bench/bench_pgen.py and the test_*_benchmark.py suites (RUN_BENCHMARK=1).

Capabilities (see the User guide and the API reference)

• IO — canonical clonotype frame on AIRR junction columns (junction_nt / junction_aa); readers for native vdjtools, AIRR Rearrangement TSV, and Parquet, plus format-detecting converters for MiXcr (v1/2 + v3/4, incl. C-gene / BCR isotype), MiGec, Adaptive immunoSEQ (v1/v2), IMGT/HighV-QUEST, Vidjil, RTCR, TRUST4, and arda AIRR output (vdjtools.io.convert); metadata-driven batch + hive-partitioned cohorts.
• Model — native V(D)J recombination model: generation probability (Pgen — nt, aa, 1-mismatch, V/J-agnostic, thread-parallel batch), sequence generation, and EM inference, all in a native (pybind11) core. Supersedes OLGA and IGoR: arda-driven scenario enumeration, polars marginal tables, read-parallelised EM, and tandem-D (D-D) support. Concordant with OLGA across all 7 loci; precomputed OLGA + real-data-learned models bundled (load_bundled).
• Model workshop (user guide) — build a model on your own V(D)J germline library (from_germline, FASTA + anchors) and fit it to your own sequences; export and re-import every marginal as tables; check a model against its germline (check_model — functional genes stuck at P=0, unreachable deletion mass, incomplete conditionals); read the EM training log and per-iteration likelihood; compare two models (per-event Jensen-Shannon / total variation, gene usage, a bnlearn-style comparison graph) and their Pgen distributions; compute log-likelihood, AIC and BIC of a clonotype set under a model; fine-tune it, extend it with a larger allele library, and re-weight V/J usage for protocol bias. Plus information content per recombination event and a total diversity estimate (human TRB: ~52 bits per rearrangement, ~45 bits per sequence, ~3·10¹³ effective sequences).
• Stats — diversity (Chao1/Shannon/Simpson/…), spectratype, V/J/VJ usage.
• Features — CDR physicochemical profiles, k-mer / V+k-mer summaries.
• Overlap — sample overlap and TCRnet (via vdjmatch/seqtree), similarity-aware overlap, clustering.
• Preprocess — downsampling, error-correction, VJ-usage batch-effect correction, pooling/joining.
• Biomarker — incidence association (Fisher / χ² / Bayesian / permutation) vs binary / HLA-allele / CMH-stratified conditions; α-β & same-chain co-occurrence pairing; metaclonotypes.
• Dynamics — longitudinal clonotype tracking between timepoints: the paired within-donor expansion test (emergent / expanded / persistent / contracted / vanishing), the VDJtrack size-bucket recapture model, metaclonotype-grouped testing, and an edgeR NB-exact caller (vdjtools.dynamics).
• Single-cell — AIRR Cell / 10x interoperability, chain pairing + QC, paired α/β Pgen, and a to_anndata bridge into the scverse ecosystem (writes .h5ad / .zarr via AnnData).

License

GPL-3.0-or-later.