vdjmatch — TCR peptide–MHC annotation

vdjmatch annotates selected T-cell receptors and repertoires against VDJdb, preserving query identities and reporting peptide–MHC candidates, competing evidence and optional background calibration. Native sequence search uses seqtree; repertoire ingestion reuses published vdjtools readers. The minimum vdjtools version is 5.0.1, which validates explicit native thread requests against the process CPU allocation.

The Python implementation follows the legacy Java/Groovy package, preserved on the legacy-java branch. See the tutorial, workflow guide and score explanation.

Pick an entry point

Task	Entry point	Guide
Annotate a repertoire or cohort	vdjmatch annotate	Tutorial
Compare retained geometries from raw inputs	sh examples/compare.sh NEW_DIRECTORY	CLI examples
Calibrate the final selected-call event	reference-calls, then calibrate-calls	CLI reference
Choose a pinned or offline reference	--vdjdb LOCAL_REFERENCE or vdjmatch update	Workflow guide

Runnable offline examples

Browse the examples gallery or the notebooks gallery for rendered previews, source downloads and raw inputs.

The examples directory contains tiny raw synthetic inputs, a pinned offline reference, annotation and comparison shell scripts, and output-free annotation and comparison notebooks. Activate your installed environment and run from the repository root:

sh examples/annotate.sh annotation-demo
sh examples/compare.sh comparison-demo

Both use one native thread, reject existing output directories and fail on CLI errors. The comparison explicitly runs exact-excluded detailed stem4 and advanced vcdr12; notebooks analyze only fresh CLI tables. DEMO_A/DEMO_B labels are synthetic software examples, not binding truth. See the example guide.

Features

• Latest/pinned GitHub releases and offline legacy, rich or AIRR reference inputs.
• Source-aware junction ingestion, stable row identities and explicit ingestion counts.
• Native batched searches by locus, with retained reference observations.
• Targeted stem4 kernel annotation; original and advanced scorers remain available.
• All peptide–MHC candidates, NED-v1 ranking and competing-label support.
• Optional finite-background calibration with a stated counting unit and search predicate.
• Paired evidence from linked reference complexes plus unlinked same-target chain marginals.
• CLI hit/candidate/call/summary/ingestion tables and a reproducibility manifest.
• Separate precursor-frequency APIs and their established worked examples.

Targeted annotation

$ vdjmatch annotate sample.tsv --threads 4 -o annotation
$ vdjmatch annotate airr_directory/ --threads 4 -o cohort
$ vdjmatch annotate --metadata samples.tsv --vdjdb reference.zip --threads 4 -o cohort

The reference is fetched automatically; --vdjdb selects a pinned offline table or archive. A metadata TSV has path and optional sample_id columns; relative paths resolve beside the metadata file. Directory discovery accepts TSV, CSV, text, gzip and Parquet tables. Sample names must be unique.

AIRR input uses anchor-inclusive junction_aa and V calls. The default scores all peptide–MHC targets with the stem4 BLOSUM62 kernel ladder. --locus paired also retains independent TRA/TRB evidence from unlinked references. Scores are relative evidence; whole-selection null calibration is required for false-call recommendations. The original scorer remains available through historical-density; advanced alternatives remain available through that command and tcrdist-neighbours.

Query format	Required sequence/gene fields	Identity/linkage
Single AIRR	junction_aa, v_call, locus (or locus inferred from V)	Optional sequence_id, query_id, clonotype_id, clone_id
Simple repertoire	junction_aa, v, locus; optional j	Optional clonotype_id, clone_id or cloneid
Linked AIRR	One TRA/TRB row per linked receptor, junction_aa, v_call, locus	Shared cell_id, pair_id, clone_id or cloneid; cell/pair linkage takes priority
Wide paired repertoire	junction_aa_alpha, junction_aa_beta, v_alpha, v_beta; optional j_alpha, j_beta	Optional clone_id or cloneid; otherwise row identity links the pair

--locus auto infers single or paired populations; --locus TRA, TRB or paired selects explicitly. Missing chains and invalid rows retain dispositions; multiple chains in one clone/locus require resolution before annotation. clonotype_id and clone_id are retained as separate identities. Bare AIRR cdr3_aa excludes conserved anchors and cannot replace junction_aa. For amino-acid clonotype inputs, collapse by the supplied V, junction and J calls, summing counts before annotation. Distinct J calls remain distinct clonotypes even when their search geometry is shared.

The default writes .queries.tsv, .selection.tsv and .manifest.json, retaining sample and source identities, competing scores, ties and no-support rows. --details additionally writes candidates and per-target evidence; its memory grows with queries × targets. --epitope PEPTIDE restricts the requested targets. --exclude-exact is for honest benchmarks; ordinary annotation retains exact hits. --kernel-scorer selects the advanced all-bounded, vonly-bounded, v-rare, all-rare or vonly-rare alternatives; v-bounded is the targeted default. The retained vcdr12 geometry is --v-regions cdr12 --stem-depth 0 --distance-matrix tcrdist --junction-ends trim3; it uses CDR1/CDR2 V similarity. Matrix extras use --distance-matrix pam250, mj or vdjam, with --stem-depth 0 --junction-ends trim3; gaps and ladder cutoffs use that matrix's units. The retained sig geometry adds --position-weighting significance-trimmed to the trim3 TCRdist geometry; sigb62 also uses --distance-matrix blosum62. These two extras require compact output and retain the recorded rounding before V-loop addition. Calibration does not transfer between geometries. Additional recorded geometries use --geometry-extra NAME. This fixed recipe sets its matrix, trim3 ends and stem term; compact output is required.

NAME	Retained geometry
flank	BLOSUM62 costs for the three residues at each junction end
framev, framemid, framej	Clipped TCRdist with the recorded frame gap prior
freegap	Clipped TCRdist without the gap-position prior
jw, stem4j	J-family mismatch cost, with zero or four stem levels
jlen, vjlen	J mismatch and length difference; vjlen also adds V-family mismatch
norm, jnorm	V-weighted population normalization, without or with J mismatch
sw, nw	Native local or global affine alignment with BLOSUM62

For example: vdjmatch annotate repertoire.tsv --geometry-extra stem4j --threads 4 --output-prefix results/stem4j. Frame recipes preserve the recorded callback argument ordering. The historical norm/jnorm denominator retains its pre-exclusion same-V target count; with --exclude-exact these are reproduction diagnostics, not a fully punctured benchmark model. Every extra needs its own calibration before confidence calls. --threads N budgets native compute threads; searches batch queries by population and do not create a worker process per receptor or input file. Mixed single/paired cohorts share marginal chain searches. V-loop tables are prepared once per invocation, including the selected matrix and V regions. The library also enforces this ceiling for the Polars table pool. If Polars was already imported with a larger pool, start a fresh process with POLARS_MAX_THREADS=N set before importing Polars; the API fails explicitly rather than exceeding the requested budget. Input metadata whose names conflict with generated output fields is retained under input_ prefixes.

from vdjmatch import annotate_repertoires

result = annotate_repertoires("airr_directory", reference="reference.zip",
                             output_prefix="cohort", threads=4)
queries, selections = result["queries"], result["selection"]

The original Annotator and annotate APIs remain available. Agent usage is documented in the repository vdjmatch skill. The synthetic library/CLI batch resource gate runs with test/bench extras: python -m pytest tests/test_annotation_resources.py -q.

For single-chain release calibration, reference-calls --scorer stem4 runs the complete selected score with full junction exclusion. Add --leave-epitope-out to remove each query's truth peptide from every model count, including shared and sparse reference labels:

$ vdjmatch reference-calls --vdjdb reference.zip --locus TRB --scorer stem4 \
    --reference-population curated-benchmark --leave-epitope-out --threads 4 \
    --output-prefix withheld
$ vdjmatch calibrate-calls withheld.calibration-cases.tsv --target-fpr 0.01 \
    --input-manifest withheld.manifest.json \
    --output-prefix calibration

This produces an empirical operating point for that pinned pipeline and null population. Sample recommendations also require matched generated nulls and reference policy; paired evidence requires its own joint calibration. The reference-calls original scorer and curated reference policy remain its defaults.

For paired nulls, use reference-calls --locus paired --scorer stem4 --leave-epitope-out. Both linked and full unlinked chain references contribute. The selected kernel_any_axis_score is the maximum of the linked, alpha and beta scores; each axis retains its own selected labels and ties.

$ vdjmatch reference-calls --vdjdb reference.zip --locus paired --scorer stem4 \
    --leave-epitope-out --threads 4 --output-prefix paired-withheld
$ vdjmatch calibrate-calls paired-withheld.calibration-cases.tsv --target-fpr 0.01 \
    --score-column kernel_any_axis_score --input-manifest paired-withheld.manifest.json \
    --output-prefix paired-calibration

This fits the declared any-axis event for the supplied null population. Repertoire recommendations also require matched generated paired nulls. Keep the reference population and exact-hit policy identical when applying a recipe; single-chain or linked-only calibration cannot be applied to the paired event.

Precursor frequency

How much repertoire mass can see an epitope, and how many cognate clonotypes exist — seen and unseen. Use the documented precursor setup (pip install 'vdjmatch[precursor]') for these examples:

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

Checking it against what laboratories measured

F(e) estimates a naive repertoire mass, and that quantity has been measured directly — by tetramer and dextramer enrichment of unexposed donors, cord blood and naive mice. Both sides of that comparison ship, and are deliberately kept apart:

from vdjmatch.precursor import load_compendium, load_estimates

comp = load_compendium()                                  # experimental: 147 published measurements,
                                                          # 14 studies, 31 originating laboratories
est = load_estimates(species="human", chain="TRB", model="arda")   # derived: what the model predicts
est.select("epitope_seq", (est["union"] * 1e6).alias("predicted"), "measured_per_1e6")

Both fetch from isalgo/airr_benchmark (vdjmatch/precursor_freq/) and cache next to the VDJdb releases. Merging them into one column would conflate a measurement with a prediction of it.

Beyond a point estimate, the same two inputs — the Pgen spectrum and the catalogued junction count — answer more than one question, and the answers are not interchangeable:

u = P.union_mass(model, junctions, r=1, count_members=True)
u["n_union"] - u["n_seqs"]                     # uncatalogued neighbours: a census, finite and exact
P.observed_mass(model, junctions) / u["union"] # share of the neighbourhood's mass already catalogued
u["overlap"]                                   # what naive per-sequence summing would double-count

Unseen mass converges; unseen count does not. The Horvitz–Thompson weight p/π → 1/N as p → 0, which is what makes coverage_corrected_mass work without guessing the shape of the tail. But 1/π diverges, so a richness extrapolation on real TCR data returns 10¹¹–10¹⁴ unseen junctions — not a number. unseen_junctions marks richness_reliable false when the rarest observed junction was captured too seldom, and the ball census above is the finite alternative.

An interactive marimo notebook walks the whole comparison — the scatter against measurement, the Pgen spectrum inside one cognate set, the ball census, and the step from a mass to a number of cells:

$ pip install 'vdjmatch[precursor,notebook]'
$ marimo edit docs/notebooks/precursor_compendium.py

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

Install

python -m pip install vdjmatch

For development, use an isolated environment:

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

Python API

import polars as pl
from vdjmatch import Annotator

ann = Annotator.latest()  # GitHub latest; .version(tag) pins a release
# ann = Annotator.from_path("reference.zip")  # existing local input, offline
queries = pl.DataFrame({
    "sequence_id": ["query-a", "query-b"],
    "junction_aa": ["CASSIRSSYEQYF", "CWWWWF"],
    "locus": ["TRB", "TRB"],
})
calls = ann.annotate(queries, threads=1)  # every input row remains
candidates = ann.candidates(queries, threads=1)  # all peptide–MHC alternatives
calibrated = ann.annotate(queries, calibrate=True, species="human", threads=1)

AIRR junction_aa contains the conserved anchors; AIRR cdr3_aa excludes them. True AIRR input without a junction is rejected. Legacy VDJdb cdr3 can already contain the junction. Custom ingestion requires an explicit junction convention; the software never fabricates anchors.

The API uses uncalibrated ranking unless calibrate=True or an explicit control index is supplied. NED ranks candidates; E is an expected background match count, and p_enrichment tests the stated reference set. None is a posterior probability of the winning specificity. Exact ranking ties abstain as ambiguous.

ann.annotate_paired(cell_df, cdr3a="cdr3a", cdr3b="cdr3b") requires both chains to match the same reference complex and peptide–MHC label. Paired background calibration states an independent-chain null explicitly. Agreement between separate alpha and beta top labels alone does not establish paired-reference evidence.

Command line

vdjmatch update --asset default --cache reference-inputs
vdjmatch match --input-format airr --threads 1 -o results/run sample.tsv
vdjmatch match --vdjdb reference.zip --no-evalue -o results/offline sample.tsv
vdjmatch match --vdjdb reference.zip --paired --link cell_id -o results/cells cells.tsv
vdjmatch match --vdjdb reference.zip --search-mode ball -o results/balls sample.tsv

CLI matching requests calibration by default; --no-evalue selects uncalibrated ranking. The command emits hits, candidates, row-preserving calls, descriptive summaries and ingestion diagnostics as TSV, plus a JSON manifest.

Use --input-format custom --sequence-convention junction for an ambiguous custom table. --scope s,i,d,t sets substitution/insertion/deletion/total budgets. Hard --match-v or --match-j requests require --no-evalue with sequence-only backgrounds. Native search defaults to one thread. See vdjmatch match -h and the CLI reference.

Ranked mode scans the complete eligible reference with a native restricted gap-block scorer, retaining at most K unique junctions or genuine paired keys before expanding observations. Its global fixed-K enrichment test is separate from descriptive epitope ranks. Ball mode reports radii1–5 separately. Neither selects the best K/radius after seeing the results. These new modes require seqtree1.0.2; see the CLI reference for scoring and null assumptions.

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