API reference#

The public API is re-exported from the top-level vdjmatch package — vdjmatch.api.Annotator and the module-level vdjmatch.annotate() shortcut. The submodules below expose the building blocks (VDJdb access, the search index, E-values, I/O).

vdjmatch.api#

High-level annotation API.

One ergonomic entry point over the VDJdb reference + seqtree search engine, from the simplest (list[CDR3] -> hits) to the fully general (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() # pinned HF release (or .version(tag)) ann.hits([“CASSIRSSYEQYF”, “CASSLAPGATNEKLFF”]) # list -> long per-hit polars frame ann.annotate(df, cdr3=”junction_aa”, locus=”locus”) # df -> df + vdjmatch_* columns ann.annotate_paired(cell_df, cdr3a=”cdr3_alpha_aa”, cdr3b=”cdr3_beta_aa”) # paired alpha+beta

vdjmatch.annotate([“CASSIRSSYEQYF”, …]) # module-level shortcut (cached default ref)

class vdjmatch.api.Annotator(index)[source]#

Bases: object

A reusable, indexed VDJdb reference. Build once, annotate many query sets.

Parameters:

index (VdjdbIndex)

classmethod latest(*, species='HomoSapiens', source='hf')[source]#

Pinned HF benchmark release (source="hf") or the latest GitHub release.

Parameters:

  • species (str | None)

  • source (str)

Return type:

Annotator

classmethod version(tag, *, species='HomoSapiens')[source]#
Parameters:

  • tag (str)

  • species (str | None)

Return type:

Annotator

classmethod from_path(path, *, species=None)[source]#
Parameters:

species (str | None)

Return type:

Annotator

classmethod from_frame(vdj, *, species=None)[source]#

Custom VDJdb-like frame (normalized columns: gene,cdr3,v,j,epitope,mhc_*,...).

Parameters:

  • vdj (polars.DataFrame)

  • species (str | None)

Return type:

Annotator

property loci: list[str]#
hits(cdr3s, *, locus='TRB', scope='1,0,0,1', match_v=False, match_j=False, **kw)[source]#

list[CDR3] (or one string) -> long per-hit frame (query/db CDR3, epitope, score, edits).

Parameters:
Return type:

polars.DataFrame

annotate(data, *, cdr3='cdr3', v=None, j=None, locus=None, scope='1,0,0,1', prefix='vdjmatch_', match_v=False, match_j=False)[source]#

Annotate a list[CDR3] or a polars frame, returning the input with appended {prefix}epitope/mhc_class/score/n_hits columns. All input rows/columns/ids/labels are kept (we annotate unique CDR3s per locus and join back). locus may name a column (per-row locus) or be a single locus string (default TRB).

Parameters:
Return type:

polars.DataFrame

annotate_paired(data, *, cdr3a='cdr3_alpha_aa', cdr3b='cdr3_beta_aa', scope='1,0,0,1', prefix='vdjmatch_')[source]#

Paired alpha+beta: annotate each chain, then call the epitope supported by both chains (intersection; score = sum of the two chain scores). Rows/ids/labels preserved. (The control-calibrated paired E-value lives in match.PairedVdjdbIndex.annotate_pairs.)

Parameters:

  • data (polars.DataFrame)

  • cdr3a (str)

  • cdr3b (str)

  • scope (str)

  • prefix (str)

Return type:

polars.DataFrame

vdjmatch.api.annotate(data, **kw)[source]#

Module-level shortcut: annotate against the default (latest pinned) reference. See Annotator.annotate().

Return type:

polars.DataFrame

vdjmatch.db#

VDJdb fetch, cache, and parse.

vdjmatch.db.fetch_latest(asset='slim', cache=None, pin=None, force=False)[source]#

Download (and cache) a VDJdb table from the latest release (or pin tag).

Parameters:

  • asset (str) – which table to extract — "slim" (default), "full", or "default".

  • cache (str | PathLike | None) – cache directory (default ~/.cache/vdjmatch or $VDJMATCH_CACHE).

  • pin (str | None) – pin a specific release tag for reproducibility; None = latest.

  • force (bool) – re-download even if cached.

Returns:

Path to the extracted table file, named vdjdb-<tag>.<asset>.txt in the cache.

Return type:

Path

vdjmatch.db.fetch_hf(tag='2026-06-11-ZENODO', asset='default', cache=None, force=False, repo='isalgo/airr_benchmark')[source]#

Fetch a pinned VDJdb table from the HuggingFace mirror (isalgo/airr_benchmark).

The benchmark releases are mirrored gzipped at vdjdb/vdjdb-<tag>/<table>.txt.gz. This is the reproducible benchmark source (GitHub release assets can change); asset selects the table ("default" = vdjdb.txt, the canonical long per-record format). Decompresses to the same cache name as fetch_latest() so the two are interchangeable.

Parameters:
Return type:

Path

vdjmatch.db.load(source=None, *, asset='slim', species=None, gene=None, mhc_class=None, min_score=0, paired_only=False, pin=None)[source]#

Parse a VDJdb table into a normalized polars frame (see schema.CANONICAL).

source is a path to an existing table (e.g. a local snapshot); if None the latest release is fetched (asset/pin). Optional filters: species (e.g. "HomoSapiens"), gene ("TRA"/"TRB"), mhc_class ("MHCI"/"MHCII"), min_score (VDJdb confidence), paired_only (keep only rows with a non-zero complex_id).

Parameters:
Return type:

polars.DataFrame

vdjmatch.db.replicated(df, min_refs=2)[source]#

High-confidence shortlist: clonotype-epitope associations independently reported in at least min_refs distinct references. Returns unique (gene, cdr3, v, j, epitope, mhc_class) rows whose (gene, cdr3, v, j, epitope) key is backed by >= min_refs distinct reference_id, with an n_refs column. Reference-replicated entries are the most trustworthy labels for benchmarking against the latest VDJdb release (cf. the mhcmatch shortlist).

Parameters:

  • df (polars.DataFrame)

  • min_refs (int)

Return type:

polars.DataFrame

vdjmatch.match#

Matching engine, scope parsing, scoring, CIGAR, paired-chain index.

class vdjmatch.match.VdjdbIndex(by_gene)[source]#

Bases: object

Searchable VDJdb partitioned by gene (TRA/TRB). Build once, reuse across samples.

Parameters:

by_gene (dict[str, tuple[Index, pl.DataFrame, pl.DataFrame]])

classmethod build(df, species=None)[source]#
Parameters:

  • df (polars.DataFrame)

  • species (str | None)

Return type:

VdjdbIndex

property genes: list[str]#
index_for(gene)[source]#
Parameters:

gene (str)

Return type:

Index | None

annotate(queries, params, *, gene, threads=0, match_v=False, match_j=False, align=False, region_aware=False, progress=False, chunk=2000)[source]#

Annotate single-gene query clonotypes; returns a long per-hit frame.

queries schema: cdr3, v, j, locus, count (one gene). Output: query_*, db_* (epitope/mhc/…), n_subs/n_ins/n_dels, score, and (if align) cigar/match.

progress shows a tqdm bar over the query batch; search_batch has no callback, so the queries are chunked (chunk per step). When progress is False the whole batch is searched in a single call (no behaviour change).

Parameters:

  • queries (polars.DataFrame)

  • params (seqtree.SearchParams)

  • gene (str)

  • threads (int)

  • match_v (bool)

  • match_j (bool)

  • align (bool)

  • region_aware (bool)

  • progress (bool)

  • chunk (int)

Return type:

polars.DataFrame

class vdjmatch.match.PairedVdjdbIndex(a_idx, b_idx, a_to_cplx, b_to_cplx, cplx_epitope, n_pairs)[source]#

Bases: object

Per-chain seqtree indices over the paired VDJdb complexes, with α/β CDR3 → complex maps.

classmethod build(vdjdb, species=None)[source]#

Build from a full VDJdb frame (needs complex_id pairing; TRA/TRB rows).

Parameters:

  • vdjdb (polars.DataFrame)

  • species (str | None)

Return type:

PairedVdjdbIndex

annotate_pairs(pairs, control_a, control_b, params, threads=0)[source]#

Joint E-value per query pair. pairs needs cdr3a, cdr3b (+ optional epitope ground truth). Returns per-pair n_joint, E, p_joint (Poisson), p_fisher, and the predicted epitope (modal among joint matches).

Parameters:

  • pairs (polars.DataFrame)

  • control_a (seqtree.Index)

  • control_b (seqtree.Index)

  • params (seqtree.SearchParams)

  • threads (int)

Return type:

polars.DataFrame

vdjmatch.match.parse_scope(spec)[source]#

Parse "s,i,d,t" / "s,id,t" / "s" → (subs, ins, dels, total).

Parameters:

spec (str)

Return type:

tuple[int, int, int, int]

vdjmatch.match.search_params(scope='1,0,0,1', *, engine='seqtm', matrix='', pos_matrix=None, max_penalty=0, mode='all', gap_open=1, gap_extend=1)[source]#

Build a seqtree SearchParams from a scope spec and scoring options.

pos_matrix is a settable attribute (not a constructor arg) in seqtree; result-count limiting (top-k) is done in Python, as seqtree has no max_hits. When a scaled matrix is used with indel scopes, set gap_open/gap_extend to the matrix scale (~a typical substitution penalty) so gaps aren’t absurdly cheap relative to substitutions.

Parameters:
Return type:

seqtree.SearchParams

vdjmatch.match.load_vdjam(path=None, scale=100)[source]#

Build a seqtree SubstitutionMatrix from a VDJAM-format similarity table.

Parameters:

  • path (str | PathLike | None) – VDJAM aa.1 aa.2 score TSV (default: bundled vdjam.txt).

  • scale (int) – integer scale factor applied to the float similarities.

Return type:

seqtree.SubstitutionMatrix

vdjmatch.evalue#

Control-calibrated E-values (single-chain first-hit + paired α/β).

vdjmatch.evalue.query_evalues(target, control, queries, params, threads=0, exclude_exact=False)[source]#

Per-query E-values → polars frame with columns query_cdr3, n_target, n_control, E, p_any, p_enrichment, rule_of_three (one row per input query, input order preserved).

exclude_exact=True punctures distance-0 (self/duplicate) hits on both sides — use when queries may be members of the target/control (e.g. VDJdb-vs-VDJdb verification).

Parameters:

  • target (seqtree.Index)

  • control (seqtree.Index)

  • queries (list[str])

  • params (seqtree.SearchParams)

  • threads (int)

  • exclude_exact (bool)

Return type:

polars.DataFrame

vdjmatch.evalue.background(locus='TRB', species='human', size=None, cache_dir=None)[source]#

Load a deduplicated background repertoire Index for the given locus/species.

Bundled: human TRB. Others (human TRA, mouse TRA/TRB) download via seqtree[control].

Parameters:

  • locus (str)

  • species (str)

  • size (int | None)

  • cache_dir (str | None)

Return type:

seqtree.Index

vdjmatch.aggregate#

Epitope-level enrichment aggregation.

vdjmatch.aggregate.epitope_summary(hits)[source]#

Epitope-level summary: unique matched query clonotypes, reads, and best alignment score, grouped by (epitope, mhc_class, antigen_species). Sorted by unique clonotypes desc.

Parameters:

hits (polars.DataFrame)

Return type:

polars.DataFrame

vdjmatch.aggregate.best_call(hits, evals=None)[source]#

Per query clonotype, the predicted epitope = most-supported (most DB records, then highest VDJdb confidence, then best alignment). If evals (query_evalues output) is supplied, attach the query’s E / p_enrichment as the annotation confidence.

Parameters:

  • hits (pl.DataFrame)

  • evals (pl.DataFrame | None)

Return type:

pl.DataFrame

vdjmatch.io#

Query repertoire I/O (AIRR rearrangement + paired cell / TCRvdb).

vdjmatch.io.read_rearrangement(path, dedup=True, valid_aa=True)[source]#

Single-chain query → schema cdr3, v, j, locus, count, pair_id. If dedup (default), collapse to unique (cdr3, v, j, locus) summing count (required for E-values).

Parameters:
Return type:

polars.DataFrame

vdjmatch.io.read_cell(path, link=None, valid_aa=True)[source]#

Paired long-form: one row per chain, α/β linked by link (default: auto = cell_id / clone_id, resolved to pair_id). Pivots to one row per cell with cdr3a/va/ja and cdr3b/vb/jb (locus prefix TRA→a, TRB→b).

Parameters:
Return type:

polars.DataFrame

vdjmatch.io.read_tcrvdb(path, valid_aa=True)[source]#

TCRvdb wide-form CSV: α/β CDR3s in one row. → cdr3a, va, ja, cdr3b, vb, jb, epitope, mhc (epitope/mhc kept as ground-truth labels for benchmarking).

Parameters:
Return type:

polars.DataFrame

vdjmatch.cluster#

Pairwise sample overlap.

vdjmatch.cluster.overlap(a, b=None, scope='1,0,0,1', matrix=None, threads=0)[source]#

Fuzzy-matching pairs between two CDR3 sets (or within one if b is None).

Returns a long frame a_idx, a_cdr3, b_idx, b_cdr3, score, n_subs of all within-scope pairs. Within-set mode (b is None) drops the trivial self-pairs (i == j).

Parameters:
Return type:

polars.DataFrame

vdjmatch.cluster.overlap_metrics(a, b, scope='1,0,0,1', threads=0)[source]#

Summary overlap metrics between two repertoires: number of matched pairs, and the fraction of each set with at least one fuzzy match in the other.

Parameters:
Return type:

dict[str, float]

vdjmatch.precursor#

Note

Needs the optional extra: pip install 'vdjmatch[precursor]'. It pulls vdjtools for the recombination model; importing vdjmatch itself never does.

T-cell precursor frequency for an epitope: how much repertoire mass can see it.

The estimand is

F(e) = sum over the cognate set C_e of pi(tau)

— the probability that a random naive-repertoire junction recognises epitope e. This is the continuous quantity behind “immunogenic”: recognition is a spectrum, and F(e) is where on it a given epitope sits.

Seven estimators of the same F(e). One is empirical and needs no model at all; the rest ride on the recombination model, in increasing order of commitment.

event_ratio

F(e) counted directly off repertoire data, as a ratio of independent recombination events to independent recombination events — no Pgen, no model. The unit is (donor, V, J, junction_nt): the same nucleotide junction in two donors is two events, two rearrangements that converged. Because numerator and denominator are counted with the same key over the same donors, repertoire depth divides out and there is no coverage correction to make. This is the estimand itself rather than a proxy for it, which makes it the direct empirical check on everything below. Its own weak point is the other sampling — the size of the cognate set VDJdb holds — which it shares with every set total here; the per-junction median form removes it, see the function’s docstring.

observed_mass

Sum of Pgen over the junctions actually recorded for the epitope. A strict lower bound, and a biased one: VDJdb samples cognate TCRs size-biased by Pgen (a TCR enters the record roughly in proportion to its repertoire frequency), so the observed members are systematically the high-Pgen ones and the deficit does not shrink with more studies at the same depth.

coverage_corrected_mass

observed_mass with the size-biased deficit put back. A capture curve increasing in Pgen is fitted to the per-junction donor/study multiplicities, then the mass is Horvitz–Thompson reweighted by each member’s inclusion probability. Unlike the bound it is not monotone in depth by construction, and it degenerates loudly rather than silently — see the function’s own docstring for what it assumes and where it breaks.

unseen_junctions

The richness counterpart of coverage_corrected_mass, from the same fitted capture curve: how many cognate junctions were never catalogued, and how rare each one is on average. Richness and mass are not interchangeable — because the unseen members are systematically the low-Pgen ones, the missing count can be enormous while the missing mass stays small.

union_mass

Mass of the union of Hamming-r balls around the observed junctions. Cognate TCRs are near-duplicates by construction, so their balls overlap and the sum of per-sequence ball masses double-counts; the union is the correct object. Exact, and without enumerating the union. The returned overlap quantifies that double-counting directly, which is worth reporting rather than hiding — it is a measurement of how tight the specificity group is. (ball_mass is the same quantity by enumeration, kept as its oracle.)

shell_profile

ball_mass resolved by exact edit distance, so the empirically measured cognacy-retention profile alpha_r can be applied per shell instead of assuming every ball member is still cognate. F ≈ sum_r alpha_r * mass(shell r).

motif_mass

Pgen of a degenerate motif (a VDJdb cluster PWM is V/J/length-pinned, hence exactly a per-position residue set — see load_cluster_motifs()). This computes the mass of a set directly, so unlike the others it does not suffer the observed-sample coverage bias at all.

Which one to use#

question

estimator

“what is F(e), measured?”

event_ratio (needs donor-resolved nucleotide repertoires)

“what can I defend without any assumption?”

observed_mass (report it as a bound)

“how much did the sampling miss?”

coverage_corrected_mass (needs >= 2 capture units and some recaptures)

“how many TCRs were never catalogued?”

unseen_junctions

“how tight is this specificity group?”

union_mass -> the overlap field

“best point estimate from observed TCRs”

shell_profile -> retained

“an estimate with no sampling bias at all”

motif_mass (needs a cluster PWM)

“how much mass is still missing?”

cross_check -> missing_fraction

“how many cells, in this donor?”

precursor_frequency -> cells

“how many clonotypes, seen and unseen?”

occupancy -> S, n_seen, n_unseen

cross_check() is the scientifically load-bearing one. motif_mass and the observed-sample estimators measure the same F(e) by independent routes with different biases, so their disagreement is an estimate of the missing mass. event_ratio() then adjudicates from outside the model entirely.

Two objects, easily confused#

coverage_corrected_mass and event_ratio both take per-junction donor counts, and they mean opposite things:

  • in coverage_corrected_mass the count is a multiplicity — how many independent units re-reported the same object, i.e. recaptures, from which what was never captured is inferred;

  • in event_ratio distinct donors are distinct objects, each one its own recombination event.

Feeding the same donor counts to both is silently wrong in one of them. “Seen again” is not “happened again”.

One-substitution matching is part of the estimator#

Everything here defaults to a closed one-substitution ball — pgen(..., mismatches=1), ball_mass(r=1), event_ratio(r=1) — and that is a definition, not a tuned parameter. At radius 0 the repertoire is too sparse for either route to be estimable, and the measured consequence is large: the Pogorelyy replication scores 0.51–0.61 at radius 0 against 0.76–0.86 at radius 1. Use radius 0 only to demonstrate that sparsity.

Requires the optional vdjtools dependency (the recombination model):

pip install 'vdjmatch[precursor]'

Nothing here reimplements Pgen — the DP, the closed Hamming-1 ball and the degenerate/masked DP all live in vdjtools, and the neighbourhood enumeration lives in seqtree. Importing vdjmatch itself never touches vdjtools; the import happens on first use inside this subpackage.

vdjmatch.precursor.load_model(locus='TRB', source='olga', organism='human')[source]#

Bundled recombination model.

vdjtools ships three sets, and they differ in ways that matter for a mass over a set of junctions, where a junction scoring exactly zero is silently dropped from the total:

"olga"

A bit-faithful import of OLGA’s published models — native Pgen matches OLGA’s own to machine precision. Human, seven loci. Faithful includes faithful to OLGA’s defects: its per-locus deletion grid puts mass on trims an allele is too short to reach, and the DP never visits those, so Pgen through the affected alleles is an underestimate or exactly zero. On VDJdb this costs 5.7% of human TRA junctions (3,000 of 52,363) and 41 of 111,331 on TRB.

"learned"

Refit from real 5’RACE reads on arda germline, so it does not inherit that grid. Human, seven loci. Loses 270 TRA junctions (0.5%) and 2 on TRB.

"arda"

The same refit on the arda IMGT allele namespace, and the only set carrying a non-human organism: human for seven loci plus mouse TRA/TRB. Loses 19 TRA junctions and 2 on TRB; neither mouse locus loses any.

"olga" remains the default because exact agreement with the reference implementation is the right property for a published null. For a mass over a set — which is what everything in this package computes — prefer ``”learned”`` or ``”arda”``, and for anything comparing human with mouse use "arda" on both so the model family is held fixed.

Sparse gene usage is not a proxy for this: mouse arda TRA has P(V) = 0 for 60% of its V genes and loses no junctions, because marginalising over V and J routes around genes the fit never saw, while olga TRA has one zero-usage V gene and loses 5.7%.

Parameters:
vdjmatch.precursor.check_junctions(seqs)[source]#

Split seqs into (junctions, suspect) by the conserved anchors.

CDR3 is not junction. VDJdb’s column is named cdr3 but holds junctions — Cys104 and Phe118/Trp included. An anchor-stripped IMGT CDR3 scores exactly 0.0 with no error, so a silently mis-typed input reports a precursor frequency of zero rather than failing. Callers should check before scoring and report the dropped count rather than letting it vanish.

vdjmatch.precursor.pgen(model, junctions, v=None, j=None, mismatches=0, threads=0)[source]#

Per-junction Pgen — the vector behind every mass in this package.

mismatches=1 returns the closed Hamming-1 ball mass of each junction in closed form (inclusion–exclusion inside vdjtools, no enumeration), which is the frequency proxy Pogorelyy et al. used. v/j are per-junction allele-resolution call lists or None to marginalise; marginal and conditioned are different quantities, so never mix them in one comparison.

Parameters:

  • mismatches (int)

  • threads (int)

Return type:

list[float]

vdjmatch.precursor.observed_mass(model, junctions, v=None, j=None, threads=0)[source]#

Sum of Pgen over the given junctions – a strict lower bound on F(e).

v/j are per-junction allele-resolution call lists (TRBV27*01, not TRBV27), or None to marginalise over V/J. The two are different quantities – the marginal is larger – so never mix them within one comparison.

Parameters:

threads (int)

Return type:

float

vdjmatch.precursor.coverage_corrected_mass(model, junctions, multiplicity, n_units=None, threads=0)[source]#

F(e) with the size-biased sampling deficit put back — the estimator, not the bound.

multiplicity[i] is the number of independent capture units (donors, or failing that studies) that re-reported junctions[i] for this epitope. Take it from donor/study counts, never from VDJdb record counts: rows duplicate the same donor’s TCR across curation passes, so a record count is not a recapture and inflating it collapses the correction. n_units is the total number of capture units for the epitope (defaults to max(multiplicity)).

Warning

multiplicity counts recaptures of one object, not objects. event_ratio() uses the opposite convention on the same-looking numbers: there each donor carrying a junction is a separate recombination event, a distinct object. Both are coherent, they are not the same quantity, and passing one function’s counts to the other is silently wrong rather than an error. If you are counting how many independent rearrangements exist, you want event_ratio(); if you are inferring what the sampling never saw, you want this.

Model. Each cognate junction is captured by one unit with probability p_i = 1 - exp(-theta * pi_i) — Poisson sampling at a rate proportional to Pgen, which is exactly the size-biasing that makes observed_mass biased. theta is fitted by maximising the zero-truncated Binomial(n_units, p_i) likelihood over the observed junctions (unobserved ones contribute nothing, hence the truncation). The corrected mass is then Horvitz–Thompson:

F_hat = sum_i pi_i / (1 - (1 - p_i)^n)

which is unbiased for the total over the whole cognate set, unobserved members included, because each observed member is upweighted by its own inclusion probability. The same inclusion probabilities give the richness, S_hat = sum_i 1 / (1 - (1 - p_i)^n), from which unseen_junctions() reads off how many members were never catalogued and how rare they are.

Returns {"observed", "corrected", "coverage", "theta", "gt_coverage", "f1", "n_units", "n_seqs", "n_zero", "n_total", "n_unseen", "unseen_mass", "mean_unseen_pgen", "degenerate", "reason"}; coverage = observed / corrected.

Why not Good–Turing. The textbook flat coverage 1 - f1/N is known-bad on TCR data — Laydon et al., PLoS Comput Biol 2014;10:e1003646 (PMID 24945836) measure 61.7% median error on real TCR abundance data, against 43.8% for Chao1bc and 42.8% for ACE, because the capture-probability distribution is far too heterogeneous for the uniform-multinomial assumption behind it. Here the inclusion probability is known rather than inferred — Pgen is the sampling probability — and the consequential behaviour is at the rare end: p/(1 - exp(-N*p)) -> 1/N as p -> 0, a constant. Every rare junction contributes the same weight however rare it is, so the estimator never has to guess the shape of the tail it cannot see, which is precisely the guess that sinks Good–Turing. The flat number is still returned as gt_coverage so the two can be compared, but it is a diagnostic, not the estimate.

Richness and mass are not interchangeable. Because the unseen junctions are systematically the low-Pgen ones, the missing count can be enormous while the missing mass stays small. F(e) is made of mass; n_unseen answers a different question and should be quoted as one.

What it assumes. (1) Capture is independent across units — false where two studies share donors or one study’s TCRs were re-curated into another, which biases theta upward and the correction downward. (2) The capture curve is the one-parameter saturating form above; real capture also depends on assay sensitivity and HLA typing of the cohort, which this folds into a single theta. (3) Pgen is the right size variable — it is the one the size-bias argument names, but clonal expansion in the source samples adds a second, unmodelled one. (4) The cognate set is fixed; a junction that is cognate in one donor and not another breaks the frame entirely.

Where it breaks, loudly. When every junction is a singleton there is no recapture information at all, the likelihood is maximised as theta -> 0 and the Horvitz–Thompson sum diverges: the function then sets degenerate=True, names the reason and returns the observed bound unchanged rather than an infinity or a ZeroDivisionError. The same happens with fewer than two capture units, and if the fit hits the search boundary.

Do not compose this with the neighbourhood route. ball_mass() puts back mass that lies near what was seen, and most of what the sampling missed is near what it saw, because cognate junctions are near-duplicates by construction. Applying both counts the same missing mass twice.

Parameters:

  • n_units (int | None)

  • threads (int)

Return type:

dict

vdjmatch.precursor.unseen_junctions(model, junctions, multiplicity, n_units=None, threads=0)[source]#

How many cognate junctions were never catalogued, and how rare each one is.

The richness counterpart of coverage_corrected_mass(), read off the same fitted capture curve: S_hat = sum_i 1/pi_i is the Horvitz–Thompson estimate of the cognate set’s true size, so n_unseen = S_hat - n_observed and the mass they carry is corrected - observed.

Returns {"n_observed", "n_total", "n_unseen", "observed_mass", "unseen_mass", "mean_unseen_pgen", "mean_observed_pgen", "rarity_ratio", "min_inclusion", "richness_reliable", "degenerate", "reason"}. rarity_ratio is how many times rarer an average unseen junction is than an average observed one — the number that makes “the missing count is huge, the missing mass is not” concrete.

Read the mass; treat the count as an order of magnitude at best. unseen_mass converges (the weight p/pi tends to a constant as p -> 0) but n_unseen does not (1/pi diverges), so a group whose rarest observed junction was barely captured extrapolates to an arbitrarily large count. richness_reliable is False when that is happening — see MIN_INCLUSION_FOR_RICHNESS. For a finite census of uncatalogued candidates, use the ball instead: union_mass(model, junctions, r)["n_union"] - n_observed.

Degenerates exactly where coverage_corrected_mass() does and for the same reasons; it reports that it cannot answer rather than returning a number.

Parameters:

  • n_units (int | None)

  • threads (int)

Return type:

dict

vdjmatch.precursor.closed_ball_mass(model, junctions, r=1, threads=0)[source]#

Mass of the closed Hamming-r ball around each junction, in closed form.

No enumeration at any radius. r=1 is vdjtools’ own pgen(..., mismatches=1). For r >= 2 the same masked transfer-matrix DP gives the ball by an alternating sum over wildcarded motifs:

m(B_r(a)) = sum_{k=0..r} (-1)^(r-k) * C(L-k-1, r-k) * sum_{|S|=k} m(W_S)

where W_S is a with the positions in S freed. The coefficient counts how many W_S a sequence at distance d from a falls into, and is constructed so the alternating sum leaves 1 for every d <= r and 0 beyond. Cost is sum_k C(L,k) DP passes rather than |B_r| Pgen calls — 106 against 33,117 at L=14, r=2.

V/J are deliberately not accepted: a substituted neighbour need not keep the centre’s V/J assignment, so conditioning the ball on the centre’s call would be wrong. This marginalises.

Parameters:
Return type:

list[float]

vdjmatch.precursor.ball_mass(model, junctions, r=1, threads=0)[source]#

Mass of the union of Hamming-r balls, by enumeration — the oracle for union_mass().

Same return shape as union_mass() minus the component fields. Materialises the whole deduplicated union as Python strings and scores every member, so it is exact but scales as 19L per centre at r=1 and ~180 L^2/2 at r=2. Prefer union_mass() for real work and keep this to regression-test it.

V/J are deliberately not accepted: a substituted neighbour need not keep the centre’s V/J assignment, so conditioning the ball on the centre’s call would be wrong. This marginalises.

Parameters:
Return type:

dict

vdjmatch.precursor.union_mass(model, junctions, r=1, threads=0, max_members=20000000, count_members=True)[source]#

Mass of the union of Hamming-r balls — exact, and without enumerating the union.

Returns {"union", "naive_sum", "overlap", "n_seqs", "n_union", "n_multiply_covered", "n_components", "n_clustered"} where overlap = 1 - union/naive_sum is the share of the naive per-sequence sum that double-counting would have invented, and n_union counts the distinct sequences the union holds — n_union - n_seqs of which are candidate cognate junctions no database has catalogued. n_union is None when counting them would exceed max_members; the masses are unaffected, since they never enumerate.

How. The naive sum counts every x in the union cov(x) = #{a : d(x,a) <= r} times, so

m(union) = sum_a m(B_r(a)) - sum_{x : cov(x) >= 2} (cov(x) - 1) * Pgen(x)

exactly, with no inclusion–exclusion and hence no truncation error (truncating I–E at pairs and triples is not safe: a component of four mutually-close junctions has a non-empty four-way term). The first sum is closed-form (closed_ball_mass()). The multiply-covered set is union_{a != b} B_r(a) ∩ B_r(b), which is empty unless d(a,b) <= 2r, so only centres inside one connected component of the 2r graph can contribute — and within a component only the members covered twice need a Pgen call. Singleton components cost nothing at all, which matters: 41.8% of VDJdb human TRB junctions are singletons.

V/J are not accepted, for the same reason as closed_ball_mass(): a substituted neighbour need not keep the centre’s V/J.

Memory. Finding the multiply-covered members counts neighbours one connected component at a time, so peak memory is the largest component’s ball union rather than the whole set’s — for a set whose junctions are all mutually distant it is zero. A single component above max_members raises rather than thrashing; split that epitope, or drop r.

Parameters:
Return type:

dict

vdjmatch.precursor.shell_profile(model, junctions, r=1, alpha=0.1, threads=0)[source]#

union_mass() resolved by exact edit distance, with cognacy retention applied per shell.

A ball at radius r treats a junction r substitutions from an observed cognate TCR as fully cognate, which it is not. Shell k is the set of sequences whose distance to the nearest observed junction is exactly k, and the retained estimate is

F ~= sum_k alpha**k * mass(shell k)

with alpha the per-edit cognacy retention, default ALPHA_PER_EDIT (0.1, Mayer & Callan 2023). alpha=1 reproduces the raw union; alpha=0 collapses to observed_mass().

The shells are obtained by differencing unions, mass(shell k) = union(r=k) - union(r=k-1), which is exact because the min-distance shells partition the ball. So nothing is enumerated and there is no memory ceiling — the r=2 profile that used to cost ~9.9M materialised strings for 300 junctions costs r+1 calls to union_mass().

Returns {"shells": [{"r", "n", "mass", "alpha"}...], "retained", "union", "n_union", "n_seqs", "alpha", "overlap"}. A shell’s n is None when the union was too large to census; its mass never is, because the masses do not enumerate.

Parameters:
Return type:

dict

class vdjmatch.precursor.RecombinationEvent(donor, v, j, junction_nt, junction_aa)[source]#

Bases: object

One independent recombination event: the key is (donor, v, j, junction_nt).

Not a sequence, not a clonotype, not a person. The same junction_nt observed in two donors is two events — two rearrangements that happened to converge on the same nucleotide string — and must never be deduplicated across donors. The same nucleotide junction rearranged onto a different V is likewise a separate event.

junction_aa is carried only for matching against a cognate set; it is deliberately not part of the key, because convergent recombination reaches one amino-acid junction by many nucleotide paths and collapsing them would undercount exactly the events this estimator exists to count.

Construction validates the pair so a mis-typed call fails loudly: junction_nt must be ACGTN-only and exactly 3x the length of junction_aa, and donor must be non-empty. Passing an amino-acid junction as the nucleotide key, or pooling donors by leaving donor blank, are the two mistakes that silently produce a wrong answer rather than an error.

Parameters:
donor: str#
v: str#
j: str#
junction_nt: str#
junction_aa: str#
property key: tuple[str, str, str, str]#
vdjmatch.precursor.event_ratio(cognate, events, r=1, denominator=None, max_members=2000000)[source]#

F(e) counted directly off repertoire data — events over events, no Pgen at all.

F_hat(e) = #{distinct (donor, V, J, junction_nt) matching C_e within r mismatches}
           ----------------------------------------------------------------------
           #{distinct (donor, V, J, junction_nt) in the whole dataset}

This is the estimand — the probability that a random naive-repertoire rearrangement recognises e — not a proxy for it, so it is a direct empirical check on the Pgen route rather than a correlate of it.

Two samplings, and only one of them divides out. Repertoire depth does: numerator and denominator are counted with the same key over the same donors, so sequencing more deeply changes both and nothing needs coverage-correcting. That is the advantage over coverage_corrected_mass(), which needs recaptures and is undefined when every junction is a singleton. Database depth does not: cognate is whatever VDJdb happens to hold for the epitope, and a bigger set mechanically matches more events. Measured on 138 epitopes, the set total f_hat is predicted at out-of-fold R^2 ~0.5 by log cognate-set size alone — observed_mass() and union_mass() carry exactly the same defect, for exactly the same reason, since all three are set totals.

The size-invariant form is a median over junctions, and it needs no new function: pass a mapping of singletons, event_ratio({tau: [tau] for tau in cognate}, events), and take the median of the per-junction f_hat. That estimates the per-junction event rate, which is a property of the epitope rather than of how many of its TCRs were catalogued, keeps the events-over-events property exactly, and throws away no data. It is also the empirical counterpart of the median rearrangement probability that Pogorelyy et al. correlated — the quantity the replication gate is built on — which is why the two agree. Prefer it to dividing the set total by len(cognate): cognate junctions are near-duplicates, their one-substitution balls overlap, and dividing a union by a count under-corrects tight groups and over-corrects diverse ones.

cognate is either one epitope’s junction list or a {epitope: junctions} mapping; the mapping form streams events once and shares the denominator across epitopes. events is an iterable of RecombinationEvent. r is the match radius in amino-acid substitutions. denominator overrides the counted total, for callers that pre-filtered the event stream; it must have been counted with the same key.

Returns {"denominator", "matched", "f_hat", "epitopes": {e: {"matched", "f_hat"}}, "n_cognate", "r"}.

`r = 1` is part of the estimator, not a tuning knob. At r = 0 the repertoire is far too sparse for this ratio to be estimable: a cognate set of a few hundred exact junctions matches almost nothing in a bulk repertoire, and the counts are then dominated by whether one public clonotype happened to be sequenced. Measured on the Pogorelyy replication the same correlation is 0.51–0.61 at r = 0 and 0.76–0.86 at r = 1. Both this estimator and the Pgen route (pgen(..., mismatches=1)) therefore take the closed one-substitution ball as their unit.

Counting objects, not recaptures. coverage_corrected_mass() takes a multiplicity that means “how many independent units re-reported this same junction” — recaptures of one object, used to estimate what was missed. Here distinct donors are distinct objects, each contributing its own event. The two are not the same quantity and the same donor counts must not be fed to both: a donor count passed as multiplicity says “seen again”, while the same number here says “happened again”.

Parameters:

  • r (int)

  • denominator (int | None)

  • max_members (int)

Return type:

dict

class vdjmatch.precursor.ClusterMotif(cid, epitope, gene, species, v, j, length, size, allowed)[source]#

Bases: object

One VDJdb cluster PWM as a degenerate motif ready for motif_mass().

allowed is one string of permitted residues per position ("" = wildcard). The cluster is V/J/length-pinned, so v/j are the conditioning to pass alongside it.

Parameters:
cid: str#
epitope: str#
gene: str#
species: str#
v: str#
j: str#
length: int#
size: int#
allowed: tuple[str, ...]#
vdjmatch.precursor.load_cluster_motifs(path, threshold=0.008, species=None, gene=None, epitope=None, min_size=2)[source]#

Read VDJdb’s motif_pwms.txt into per-position allowed-residue sets.

The file is one row per (cluster, position, residue) with freq the residue’s frequency within the cluster. Because a cluster is pinned to one V, one J and one length, thresholding freq per position yields the allowed argument motif_mass() wants — no alignment, no register search, no enumeration.

threshold is the per-position frequency cut, default MOTIF_FREQ_THRESHOLD. A position never comes back empty: if no residue clears the cut the modal residue(s) are kept, so raising the threshold shrinks the motif monotonically towards the consensus sequence rather than zeroing its mass.

The file is not a complete count table, and assuming it is silently drops real cluster members. Measured on the 2026-06 release: 2,326 of 24,036 listed positions (9.7%), spread over 1,042 of 1,791 clusters (58%), have listed frequencies summing to less than 1, and 172 positions are missing outright — residues the release filtered out. In one worked case (H.B.ATDALMTGY.5) the position-12 row lists only F at freq=0.4 while 6 of the 10 member junctions carry Y there. So: whenever the listed frequencies fall short of 1 by more than threshold, the unlisted residues could individually clear the cut and there is no way to know which they are, so the position becomes a wildcard rather than a set that excludes known members. Positions absent from the file are wildcards for the same reason.

species/gene/epitope filter exactly; min_size drops clusters below that many members (csz). Accepts a plain or gzipped path.

Parameters:

  • threshold (float)

  • species (str | None)

  • gene (str | None)

  • epitope (str | None)

  • min_size (int)

Return type:

list[ClusterMotif]

vdjmatch.precursor.motif_mass(model, allowed, v=None, j=None)[source]#

Pgen of every junction matching a degenerate motif.

allowed is one entry per position, each a string of permitted residues; "" or "X" means any residue. A VDJdb cluster PWM is V/J/length-pinned, so thresholding it per position gives exactly this — and one call returns the whole cluster’s mass with no enumeration and no inclusion–exclusion.

Unlike observed_mass() and union_mass() this scores a set, so it carries no observed-sample coverage bias. Pass the motif’s own v/j — cluster motifs are V/J-pinned and the conditioned quantity is the right one for them.

Return type:

float

vdjmatch.precursor.cross_check(model, junctions, allowed, r=1, alpha=0.1, threads=0)[source]#

Two independent estimates of the same F(e) — and their disagreement is the missing mass.

Route A (motif_mass()) scores a set: it asks the recombination model for the total mass of every junction the motif admits, so it never touches the observed sample and carries none of its coverage bias. Route B (observed_mass() / shell_profile()) starts from the junctions VDJdb actually recorded, so it is bounded below by how deeply the epitope was sampled. Both are estimates of F(e) for the same epitope.

Returns {"set_mass", "observed_mass", "ball_mass", "retained_mass", "ratio_observed", "ratio_retained", "missing_fraction", "n_seqs"}.

Interpretation. ratio_observed = set_mass / observed_mass is the factor by which the sample under-counts, and missing_fraction = 1 - observed_mass/set_mass the share of the cognate mass never observed. ratio_retained repeats it against the shell-weighted estimate: if the neighbourhood correction is doing its job, ratio_retained is closer to 1 than ratio_observed — that is the whole claim of the ball/shell route, tested rather than assumed. A ratio below 1 is informative in the other direction: the motif is tighter than the sample it was built from, i.e. the threshold is too strict or the cluster is a proper subset of the epitope’s cognate TCRs (which it usually is — one epitope has several clusters).

Both routes are marginalised over V/J so the two numbers are the same quantity. Cluster motifs are V/J-pinned and motif_mass() will condition on request, but a conditioned A against a marginal B is a category error and is not offered here.

Parameters:
Return type:

dict

vdjmatch.precursor.precursor_frequency(model, junctions, r=1, alpha=0.1, q=1.0, n_cells=100000000000.0, compartment=1.0, n_eff=None, threads=0)[source]#

F(e) as a frequency, with the cell count and precursor probabilities that follow.

F_hat(e) = q * sum_{k=0..r} alpha^k * mass(shell k)
cells    = n_cells * compartment * F_hat(e)
lambda   = n_eff * F_hat(e)        P(>=m precursors) = 1 - Poisson_cdf(m-1; lambda)

The mass term is shell_profile()’s retained: the union of the observed junctions’ Hamming-r balls, resolved into min-distance shells and down-weighted by the measured cognacy retention alpha per edit, so a neighbour k substitutions away contributes alpha^k of its mass rather than all of it. Shells partition the union, so this neither double-counts the overlap between two cognate junctions’ balls nor assumes every neighbour is still cognate.

q is the selection constant that carries a generation probability to a post-selection repertoire frequency. It defaults to 1.0, i.e. uncalibrated — the returned F is then a raw model mass and only its ranking is meaningful. Pass ALICE_Q for the published TRB value, or a constant fitted against your own cohort. Whatever you pass is echoed back in the result so a reported number always carries its calibration.

compartment is the fraction of n_cells the epitope’s restriction actually addresses — for an MHC-I epitope, the CD8 fraction, times the naive fraction if you are after a naive precursor count rather than a realised one. Left at 1.0 the cells field is “cells in the whole T-cell pool”, which is rarely the quantity a tetramer experiment measured.

n_eff is the number of independent rearrangements the probabilities are about, and it is not the same as n_cells: set it to a donor’s distinct-clonotype count to ask “does this person have a precursor at all”, or to a sample’s sequencing depth to ask “will I see one in this tube”. Left as None the Poisson fields are omitted rather than guessed.

Returns {"F", "q", "alpha", "r", "union", "retained", "overlap", "n_seqs", "shells", "cells", "n_cells", "compartment", "lambda", "p_ge_1", "p_ge_10"}, with the last three None when n_eff is not given.

What `F` is not. It is the probability that a precursor exists, not that a response is mounted: a detectable response needs more than one cell, which is what the p_ge_* fields are for, and it needs priming, help and an absence of tolerance, which nothing here models.

Parameters:
Return type:

dict

vdjmatch.precursor.occupancy(model, junctions, r=1, alpha=0.1, n_eff=None, selection=1.0, sample=3000, seed=0, threads=0)[source]#

How many cognate clonotypes exist, how many are visible at depth n_eff, and their mass.

S           = sum_k alpha^k n_k                          effective cognate-set size
F           = sum_k alpha^k m_k                          fraction of the naive repertoire
n_seen(N)   = sum_k alpha^k n_k E_k[1 - exp(-N Pgen)]
n_unseen(N) = S - n_seen(N)

with n_k and m_k the size and Pgen mass of shell k of the neighbourhood, and alpha the per-edit cognacy retention.

Why this rather than the summed mass. A specificity database samples cognate TCRs size-biased by repertoire frequency, p_i = pi_i / F. Under that sampling the arithmetic sum over a catalogued set has expectation

E[ sum_{i in sample} pi_i ] = n * (sum_i pi_i^2) / F

— proportional to how many TCRs were catalogued, proportional to how concentrated the spectrum is, and inversely proportional to the quantity it is meant to estimate. A set total is therefore a measurement of database attention as much as of biology, and correlations built on one vanish when the catalogued count is controlled for. The occupancy form has no such term: it is a count of sequence-space members weighted by how likely each is to be generated.

n_seen saturates in n_eff, which is the substantive difference. A cognate set concentrated on a few high-Pgen junctions is exhausted at shallow depth; a broad one keeps accumulating clonotypes as depth grows. Two epitopes with the same summed mass and different shapes therefore have different precursor counts, and it is the count that a response is built from.

selection multiplies the depth: a repertoire of n_eff observed rearrangements behaves like selection * n_eff draws from the generation distribution. It defaults to 1.0, i.e. uncalibrated. SELECTION_BY_CHAIN carries the measured per-chain values, which differ (TRB 4.62, TRA 1.07) and should not be pooled.

Returns {"S", "F", "n_seen", "n_unseen", "seen_fraction", "n_eff", "selection", "alpha", "r", "shells"}, with the n_seen fields None when n_eff is not given.

Counting distinct junctions is deliberate. Nothing here consumes a record or donor multiplicity: how often a clonotype was reported is expansion and curation attention, not how many clonotypes exist.

The shell sizes are exact; the Pgen spectrum within a shell is estimated from a uniform subsample of sample members, because a radius-1 shell around a few thousand junctions holds millions of sequences and only its distribution enters. seed fixes that subsample.

Parameters:
Return type:

dict

vdjmatch.precursor.expected_cells(f, n_cells=100000000000.0, compartment=1.0)[source]#

Expected number of epitope-specific cells: n_cells * compartment * f.

compartment restricts n_cells to the subset the epitope can address (the CD8 fraction for MHC-I, times the naive fraction for a naive precursor count). The result inherits the error of both f and the assumed pool size, so quote it as an order of magnitude.

Parameters:
Return type:

float

vdjmatch.precursor.p_at_least(f, n_eff, k=1)[source]#

P(X >= k) for X ~ Poisson(n_eff * f) — the probability that k precursors exist.

F(e) on its own answers “is there at least one?”, i.e. k = 1. A detectable assay or clinical response needs more, and once n_eff * F is of order one the two stop being a monotone reparametrisation of each other — which is exactly the regime a rare neoepitope sits in. That is why k is an argument rather than a fixed 1.

n_eff is a count of independent rearrangements; see precursor_frequency().

Parameters:
Return type:

float

vdjmatch.precursor.paired_frequency(f_alpha, f_beta, kappa=1.0)[source]#

Frequency of cells whose both chains are cognate: kappa * f_alpha * f_beta.

Alpha and beta pair essentially without constraint across the repertoire, so unconditionally P(A, B) = P(A) P(B). Given an epitope they do not: a cognate beta works with a restricted set of alphas, not with any alpha drawn from the cognate alpha set. kappa carries that departure — kappa = 1 asserts fully permissive pairing within the cognate sets and is an upper bound; kappa < 1 is restricted pairing. It is measurable on VDJdb’s paired records and is not assumed here, which is why it defaults to the bound and must be passed to claim anything tighter.

Both inputs must be the same kind of quantity — two calibrated F values, or two raw masses, never one of each.

Parameters:
Return type:

float

vdjmatch.precursor.summarise(frame, *, junction_col='cdr3', group_col=None, chain_col=None, capture_col=None, locus='TRB', source='olga', organism='human', r=1, alpha=0.1, q=1.0, n_cells=100000000000.0, compartment=1.0, n_eff=None, selection=1.0, min_junctions=1, threads=0, progress=False)[source]#

Run summarise_group() over every group in frame.

group_col names the grouping column (epitope on a normalised VDJdb table); without it the whole frame is one group. chain_col splits by locus and loads the matching recombination model per chain, so a mixed TRA/TRB table is scored correctly rather than with one model. capture_col (reference_id on VDJdb) supplies the recapture multiplicities the unseen-species fields need — the number of distinct units that re-reported each junction.

Parameters:
Return type:

polars.DataFrame

vdjmatch.precursor.summarise_group(model, junctions, *, group='', chain='', r=1, alpha=0.1, q=1.0, n_cells=100000000000.0, compartment=1.0, n_eff=None, selection=1.0, multiplicity=None, n_units=None, threads=0)[source]#

Every estimator for one set of junctions, as a flat row keyed by SCHEMA.

junctions are junctions (Cys104..Phe/Trp118 inclusive), not IMGT CDR3s; the ones that fail the anchor check are dropped and counted in n_dropped rather than scored as 0. multiplicity — one capture count per input junction, i.e. how many independent units re-reported it — switches on the unseen-species fields; without it they are null and unseen_status says why.

Parameters:
Return type:

dict

On this page
Show Source