vdjmatch — control-calibrated TCR antigen-specificity annotation

Fast, control-calibrated annotation of T-cell receptor antigen specificity.

vdjmatch annotates clonotypes in large AIRR repertoires against VDJdb by fuzzy CDR3 search, reporting a control-calibrated E-value (BLAST-style significance against a background repertoire) and enriched antigen-specificity labels. It is a Python rewrite of the legacy Java/Groovy vdjmatch, built on the seqtree search core.

Status: early alpha (PyPI 0.2.0), under active development on dev. The "2.0" line is the Python rewrite of the legacy Java/Groovy vdjmatch (1.x), which is preserved on the legacy-java branch (tags 1.1.4–1.3.1).

Features

• Fetch the latest VDJdb release and annotate AIRR Rearrangement / Cell (paired α/β) samples.
• Extremely fast, multithreaded search of million-scale repertoires (via seqtree).
• Control-calibrated E-values — single-chain and paired α/β (vdjmatch.evalue).
• Custom substitution matrices, including segment-specific (V / NDN / J) scoring; the TCR-specific VDJAM matrix is bundled.
• Rich per-hit output: ranked hits, CIGAR + alignment match/gap, alignment scores, E-values.
• Epitope-level enrichment summaries; pairwise sample overlap.
• T-cell precursor frequency for an epitope (vdjmatch precursor, optional extra) — how much repertoire mass can see it, how much of the cognate set no database has catalogued, and how many precursor cells that implies.
• A small command-line interface (vdjmatch update / match / precursor) and a polars-native Python API.

Precursor frequency

How much repertoire mass can see an epitope, and how many cognate clonotypes exist — seen and unseen. Needs the optional extra (pip install 'vdjmatch[precursor]', which pulls vdjtools for the recombination model):

$ vdjmatch precursor --vdjdb --mhc-class MHCI --min-junctions 10 \
      --n-eff 1e8 --selection auto -o precursor.txt

from vdjmatch import precursor as P

model = P.load_model("TRB")
junctions = ["CASSIRSSYEQYF", "CASSLGQAYEQYF", ...]     # junctions, not IMGT CDR3s

P.union_mass(model, junctions)      # exact union of the 1-mm balls, and the overlap a sum invents
P.occupancy(model, junctions, n_eff=1e8, selection=P.SELECTION_BY_CHAIN["TRB"])
# -> {"S": ..., "F": ..., "n_seen": ..., "n_unseen": ..., "seen_fraction": ...}

S is the effective cognate-set size, F the fraction of the naive repertoire, and n_seen saturates in depth — a cognate set concentrated on a few high-Pgen junctions is exhausted at shallow sequencing, a broad one keeps accumulating clonotypes. SELECTION_BY_CHAIN carries the measured per-chain depth factors (TRB 4.62, TRA 1.07); they differ and should not be pooled.

Two things worth knowing before trusting a number:

• --source picks the recombination model, and it matters for a mass over a set. The default olga is a bit-faithful import of OLGA's published models — faithful to OLGA's deletion-bin grid too, under which 5.7% of human TRA junctions in VDJdb score Pgen exactly zero (up to 15% for some epitopes) and vanish from the total without an error. learned and arda are refits that do not inherit it; arda is the only set carrying mouse. The n_zero_pgen output column reports the loss per group and the CLI warns above 1%.
• --species selects records, --organism selects the model, and they must agree. A mismatch used to run and return a plausible number that was wrong by an epitope-dependent factor; since 0.3.0 it is refused.

The method, its benchmarks and the paper live outside this repository: repseq/2026-precursor-freq (benchmarks, result tables, the 16-study literature compendium) and repseq/2026-precursor-freq-ms (the manuscript). This repository holds only the software, its CLI, tests, docs and the worked example (docs/notebooks/precursor.ipynb).

Install (development)

python -m venv .venv
source .venv/bin/activate.fish
pip install -e .[test,bench]

seqtree (the search engine) is installed from PyPI as a dependency.

Python API

One ergonomic entry point — list[CDR3] → hits up to polars df → df + annotation columns (ids/labels preserved), single- or paired-chain, against any VDJdb version or a custom reference:

import vdjmatch
ann = vdjmatch.Annotator.latest()                      # or .version("2026-06-11-ZENODO") / .from_path(...) / .from_frame(df)
ann.hits(["CASSIRSSYEQYF", "CASSLAPGATNEKLFF"])        # → long per-hit polars frame
ann.annotate(df, cdr3="junction_aa", locus="locus")    # → df + vdjmatch_{epitope,mhc_class,score,n_hits}
ann.annotate_paired(cell_df, cdr3a="cdr3_alpha_aa", cdr3b="cdr3_beta_aa")
vdjmatch.annotate(["CASS..."])                         # module-level shortcut (cached default reference)

First-hit (adaptive) E-value. Significance is evaluated at each query's nearest VDJdb hit (up to 5 edits, ≤2 ins, ≤2 del): the control's neighbour count grows with the radius, so a distance-1 hit is significant while a distance-5-only hit is not — noise is rejected without a fixed scope (vdjmatch.evalue.first_hit). Edit-distance and BLOSUM+possig-penalty "nearest" both supported.

Command line

vdjmatch update                          # fetch + cache the latest VDJdb release
vdjmatch match sample.tsv                # annotate an AIRR rearrangement sample
vdjmatch match -o run/out --match-v --threads 8 *.tsv

match writes three TSV files per sample (prefix set by -o):

Key options: --scope s,i,d,t (search budget, default 1,0,0,1), --matrix {vdjam,none}, --match-v / --match-j, --no-evalue. Run vdjmatch match -h for the full list.

Scoring: what works (and what doesn't)

An empirical study on VDJdb (see appendix/vdjmatch_scoring.tex; regenerated on the 2026-06-11-ZENODO release with composition controls and balanced metrics) settles the scoring question honestly:

• Hamming distance 1 is the signal:noise optimum — macro purity (per-epitope mean) falls 0.49 → 0.07 across edit distance 1–5, a 56× → 2.5× enrichment over chance (reproducing the original VDJdb observation, Shugay et al. NAR 2018). The search-ball radius, not the substitution matrix, is the dominant lever. (On the dense 2026 release this only shows up once a few 10× mega-studies are capped and the random tail is treated as an admixed control — naive pooled purity reads a flat ~0.9.)
• Central (NDN) substitutions carry the specificity signal — a mismatch in the CDR3 core most often changes specificity (P(same epitope) ≈ 0.31) while near-anchor mismatches are germline noise (≈ 0.75); the NDN core is also ~31–34% glycine (insertion/D-gene signature).
• No amino-acid matrix clearly beats BLOSUM62 — and a genetic-code null ties it. BLOSUM62 ≈ PAM250 ≈ structural > Hamming > data-derived VDJAM. Strikingly, VDJAMr — a matrix built from the genetic code alone (how mutationally accessible one AA is from another; loo_vdjam.codon_dissim) — matches BLOSUM62 (0.581 vs 0.564 @≤2), so TCR CDR3 substitution structure is generative, not chemical. A published TCR-specific matrix, tcrBLOSUM, does not beat BLOSUM62 either (0.555 vs 0.567; bench/tcrblosum_refute.py) — its same-epitope counts don't transfer.
• Position-weighting BLOSUM62 does beat it. Encoding the central-substitution finding as a seqtree positional matrix (PositionalMatrix.from_weights(BLOSUM62, …), centre ~2× the borders) raises leave-one-out retrieval (balanced PR-AUC) above flat BLOSUM62 in 7/8 held-out epitopes (0.564 → 0.598 @≤2; 6/8 @≤4). The weight is end-anchored (offset from each germline anchor, Beta-Binomial-smoothed) and BLOSUM severity discriminates exactly where the weight is high (the core). For CDR3, where a mismatch falls matters more than which residue it is; the first-order statistic is still the control-calibrated E-value.
• The V gene is a strong, near-binary prior — partly recovered at near-exact germline identity. Same-V neighbours share the epitope ~45–64% of the time vs ~6–17% cross-V (ratio up to ~7×). Loose CDR1/CDR2 similarity barely predicts it (point-biserial r ≈ 0), but cross-V co-specificity rises monotonically as germline CDR1+CDR2 approach identity — ~11% at ≥6 mismatches to ~32% at edit-0: a real ~3× lift, but only ~half the same-V level (32% vs 64%), recovered whole-loop rather than via a sparse pseudosequence (per-position lift flat except CDR2 pos 5; edit-0 bin n≈25, noisy). So the germline loops carry part of the V prior at near-exact tolerance; the rest is gene-identity-specific (bench/vregion_decompose.py, bench/vpseudo.py; vdjmatch.match.vgene).

Benchmark

The standing benchmark is the 2026-06-11-ZENODO VDJdb release (mirrored on the isalgo/airr_benchmark HF dataset, fetched via db.fetch_hf). For the scoring studies it is composition-controlled (_bench.long_list: keep epitopes ≥30 clonotypes, cap mega-epitopes to a random 3000, drop spectratype-anomalous spike studies), with imbalance-robust metrics (bench/metrics.py: ROC-AUC + balanced PR/F1). Its highest-confidence subset — the shortlist of clonotype–epitope pairs in ≥2 independent references (db.replicated, ~3000 TRB + ~2000 TRA) — is the gold standard (as in mhcmatch), kept separate from the long-list. Leave-one-out NN annotation (bench/shortlist_accuracy.py, subs 1, exact self excluded) reaches ~44–47% top-1, modestly above single-reference controls (the dramatic 3× gap on a sparse older export was an artefact — the dense release makes controls findable too).

License

GPL-3.0-or-later (it builds on seqtree, which is GPL-3.0-or-later).