API reference#

See Annotate your first receptors for a runnable example, Work with releases and repertoires for source conventions and offline inputs, and Understand matching evidence for score definitions.

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

Batch repertoire annotation#

annotate_repertoires accepts a table, directory, list of paths or a metadata TSV with path and optional sample_id. It preserves sample/source identities and batches single or linked populations. Compact kernel selection is the default; selection_only=False also exports full candidate/per-target evidence. threads also bounds the Polars pool. Set POLARS_MAX_THREADS before importing Polars if that library has already initialized its pool; a larger existing pool raises an error before input loading. Generated-field collisions retain input metadata under input_ prefixes.

vdjmatch.annotate_repertoires(inputs=None, *, metadata=None, reference=None, output_prefix='annotation', locus='auto', species='HomoSapiens', threads=1, selection_only=True, **options)[source]#

Annotate files, a directory, or path/sample_id metadata as population batches.

Source rows are read once through shared ingestion. Query IDs are global numeric offsets; source IDs, metadata and sample-scoped pairing remain available in the returned tables. Defaults and advanced scorer options come from the existing targeted annotation parser. No per-file search pool or computed cache is used. Outputs are <prefix>.queries.tsv, .selection.tsv, optional detail tables, and .manifest.json. The result contains frames, manifest and output paths.

vdjmatch.api#

The original Annotator API’s automatic calibration derives the background organism from the queried reference loci; mixed species or an explicit mismatch raises an error, and missing metadata requires species="human" or species="mouse". Supplied control mappings must cover every active locus (both TRA and TRB for paired queries); their organism is caller-owned.

Standalone selected-receptor annotation with explicit reference and candidate evidence.

class vdjmatch.api.Annotator(index, *, provenance=None)[source]#

Bases: object

An explicitly owned reference index, reusable across queries without result memoization.

Parameters:

index (VdjdbIndex)

classmethod latest(*, species='HomoSapiens', source='github', asset='default', cache=None, pin=None)[source]#

Resolve GitHub latest/pinned release; HF benchmark mirror requires source=’hf’.

classmethod version(tag, *, species='HomoSapiens', source='github', asset='default', cache=None)[source]#
classmethod from_path(path, *, species=None, asset='default', **filters)[source]#
classmethod from_frame(vdj, *, species=None)[source]#
property reference_report#

Retained reference observations versus searchable valid junctions.

property loci#
graded_candidates(data, *, cdr3='cdr3', v=None, j=None, locus=None, threads=1, control=None, calibrate=False, species=None, matrix=None, sequence_convention=None, exclude_exact=False, return_hits=False, match_v=False, match_j=False)[source]#

Five fixed edit balls with radius-constrained weighted rankings.

Unit-cost balls share one minimum-edit search. Weighted search may choose a cheaper alignment with more edits, so each weighted radius is searched in one complete query batch. No best-radius P-value is selected. Hard V/J predicates are supported for uncalibrated ranking; sequence-only controls cannot calibrate these predicates and are rejected.

unified_candidates(data, *, cdr3='cdr3', v=None, j=None, locus=None, threads=1, control=None, species=None, exclude_exact=False, sequence_convention=None, return_hits=False, distance='edit')[source]#

Experimental fixed-ball background-mass evidence for a selected distance.

Both edit and anchored BLOSUM62 gap-block distances cover all query lengths. No length switch, label-dependent route or top-K cap. Sequence-only controls do not provide V-conditioned or paired calibration.

ranked_candidates(data, *, cdr3='cdr3', v=None, j=None, locus=None, k=10, threads=1, matrix=None, control=None, sequence_convention=None, exclude_exact=False, paired=False)[source]#

Global gap-block top-K evidence and a separate global enrichment test.

K is predeclared; candidate label scores are not calibrated probabilities. Paired input uses cdr3a/cdr3b and genuine reference pairs.

hits(cdr3s, *, locus='TRB', scope='1,0,0,1', match_v=False, match_j=False, exclude_exact=False, threads=1, **kw)[source]#

Long reference observations; internal cdr3 names denote anchor-inclusive junctions.

candidates(data, *, cdr3='cdr3', v=None, j=None, locus=None, scope='1,0,0,1', threads=1, match_v=False, match_j=False, control=None, calibrate=False, species=None, score_scale=400.0, soft_v=True, sequence_convention=None, exclude_exact=False)[source]#

All per-pMHC candidates, competing evidence and optional fixed-ball calibration.

With no control/calibrate request, ranking is explicitly uncalibrated. Supply a seqtree Index or locus→Index mapping, or calibrate=True for the versioned standard background. NED-v1 uses +1 control counts, temperature 400 and existing soft-V weight 0.25. exclude_exact punctures zero-edit hits in both target and control searches.

annotate(data, *, cdr3='cdr3', v=None, j=None, locus=None, scope='1,0,0,1', prefix='vdjmatch_', threads=1, match_v=False, match_j=False, control=None, calibrate=False, species=None, score_scale=400.0, soft_v=True, sequence_convention=None, exclude_exact=False)[source]#

Append best candidate evidence/status, preserving every input row and its order.

Exact ranking ties abstain (ambiguous). Confidence is exposed as evidence; NED and the fixed-ball enrichment P-value are never a probability of the winning label.

paired_candidates(data, *, cdr3a='cdr3_alpha_aa', cdr3b='cdr3_beta_aa', scope='1,0,0,1', threads=1, control=None, calibrate=False, species=None, align=False, score_scale=400.0, return_hits=False, progress=False, exclude_exact=False)[source]#

Genuine same-complex evidence; optional fixed-ball chain-independent calibration.

control is a TRA/TRB → seqtree Index mapping. Paired reference similarity ranks unique junction pairs. Independent observation counts do not multiply ranking votes. Joint enrichment assumes independent background chains; it is not a posterior label confidence and does not model repertoire co-occurrence (deferred research). return_hits=True returns (joint_hits, candidates) from the same searches; chain metadata in the detailed hits uses alpha_ and beta_ prefixes. exclude_exact=True removes only the exact/exact pair. Its fixed-ball calibration uses simultaneous category bounds and a binomial tail with failure budget 1e-6, conditional on IID category/target sampling and independent background chains. Controls must index unique junctions; duplicate exact memberships fail. Deduplicated controls do not establish IID.

annotate_paired(data, *, cdr3a='cdr3_alpha_aa', cdr3b='cdr3_beta_aa', scope='1,0,0,1', prefix='vdjmatch_', threads=1, control=None, calibrate=False, species=None, align=False, score_scale=400.0, exclude_exact=False)[source]#

Append same-reference-complex evidence with explicit chain-independent calibration.

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

Convenience call; provide an Annotator to reuse reference inputs explicitly.

Detailed paired output#

paired_candidates(..., return_hits=True) returns (detail, candidates) from one reference search per chain. Detailed rows retain shared query/complex/pMHC keys and chain metadata under alpha_ and beta_ prefixes. align=True adds each chain’s alignment/CIGAR fields. score_scale must be finite and positive; its default is 400 for the bundled scoring convention. Both arguments are also supported by annotate_paired.

detail, candidates = annotator.paired_candidates(
    cells, cdr3a="cdr3a", cdr3b="cdr3b", align=True,
    score_scale=400.0, return_hits=True, threads=1,
)

vdjmatch.db#

VDJdb fetch, cache, and parse.

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

Fetch latest/pinned release; default selects primary, or historical long table.

primary/airr return a ZIP suitable for load(); historical and explicit legacy/slim/full requests return an extracted table. Downloads stream into unique staging files and are checksummed and published under an OS process lock.

Parameters:
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 gzipped legacy table, with cache identity scoped to HF repo/tag/table.

Parameters:
Return type:

Path

vdjmatch.db.load(source=None, *, asset='default', species=None, gene=None, mhc_class=None, epitope=None, mhc_a=None, mhc_b=None, mhc_match='exact', reference_id=None, exclude_reference_ids=None, evidence_type=None, min_score=0, paired_only=False, pin=None, cache=None)[source]#

Load a table, release ZIP or directory; local inputs never access the network.

Metadata and independent observations survive normalization. Primary tables use stable record IDs for pairing. AIRR requires Rearrangement + Reactivity linked by cell_id and junction_aa; species filtering errors when the export has no receptor species. Peptide, MHC and study selectors accept a string or list of strings (OR within each selector, AND across selectors). Evidence selects any matching observation attached to the record and retains all of that record’s chain/evidence metadata. A requested evidence predicate requires a supplied evidence_type field. mhc_match="compatible" matches family/allele restrictions with intersecting declared resolution, retaining original restriction labels. Exact matching is the default for backwards compatibility.

Parameters:
Return type:

polars.DataFrame

vdjmatch.db.provenance(source)[source]#

Describe a resolved local reference, including fetched release identity when available.

Parameters:

source (str | PathLike)

Return type:

dict

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

Curated clonotype/pMHC keys supported by distinct nonempty reference IDs.

Retain species and complete MHC fields when present. Repeated rows from one source do not establish replication; distinct IDs do not prove independent experiments or biological binding. Junction/V/J identities are exact source values; gene_level=True explicitly strips V/J allele suffixes. Missing receptor fields cannot establish replicated clonotype support.

Parameters:

  • df (polars.DataFrame)

  • min_refs (int)

  • gene_level (bool)

Return type:

polars.DataFrame

vdjmatch.match#

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

class vdjmatch.match.VdjdbIndex(by_gene, reference=None)[source]#

Bases: object

Explicitly owned reference index; unique junction search keys, independent records retained.

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

records_for(gene)[source]#

Read-only reference observations for candidate coverage/provenance.

Parameters:

gene (str)

Return type:

polars.DataFrame

empty_hits(query_id_type=polars.UInt32, *, align=False, region_aware=False)[source]#

Typed empty evidence with the same reference metadata columns as ordinary hits.

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

Long hit evidence, preserving query_id/locus and all reference metadata.

Search is one native batch of distinct junctions even with progress enabled. threads controls seqtree; Polars performs sequential post-search stages in its existing process-wide pool. chunk is retained for signature compatibility; it does not split native search. exclude_exact removes zero-edit hits before alignment and observation expansion.

Parameters:

  • queries (polars.DataFrame)

  • params (seqtree.SearchParams)

  • gene (str)

  • threads (int)

  • match_v (bool)

  • match_j (bool)

  • exclude_exact (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 indices linked to unique paired junctions and their source observations.

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

Validate paired cardinality; retain observations and index unique paired keys.

Parameters:

  • vdjdb (polars.DataFrame)

  • species (str | None)

Return type:

PairedVdjdbIndex

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

Fixed-ball pair/chain counts and calibrated E/p values; ties sort by epitope.

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]#

Unique junction/V/J/locus clonotypes, query rows and their abundances, without database-record duplication.

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

Compatibility reducer using NED-v1; optional global statistics stay separately named.

vdjmatch.io#

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

vdjmatch.io.read_rearrangement(path, dedup=False, valid_aa=True, *, source='auto', sequence_convention=None, return_report=False)[source]#

Read original query rows with stable IDs; explicitly request key deduplication.

dedup=True sums counts per junction/V/J/locus/pair and retains original identities as query_ids and (when supplied) sequence_ids lists. return_report=True returns (frame, ingestion_counts).

Parameters:
vdjmatch.io.read_cell(path, link=None, valid_aa=True, *, source='auto', sequence_convention=None, ambiguity='error', return_report=False)[source]#

Pair TRA/TRB rows in first-cell order, preserving chain IDs and counts.

link selects the original linkage column. Multiple rows for one cell/locus raise by default; ambiguity='first' explicitly selects the first. Missing chains have typed null fields and pair_status='missing_alpha'/'missing_beta'. Wide junction_aa_alpha/beta input also retains wholly missing pairs and exposes original wide IDs as source_query_ida/b/source_sequence_ida/b.

Parameters:
vdjmatch.io.read_tcrvdb(path, valid_aa=True, *, return_report=False)[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:
vdjmatch.io.normalize_annotation_input(path, *, source='auto', sequence_convention=None)[source]#

Return canonical raw long rows and ingestion diagnostics, without dropping rows.

AIRR and simple junction tables retain source columns, identities and metadata. Wide junction_aa_alpha/beta rows expand to TRA/TRB rows in input order; original wide IDs remain in source_query_id/source_sequence_id and computational chain IDs are unique row offsets. Missing chains remain present. Explicit pair/cell links take priority over clone IDs and identify one wide pair each; omitted links use the supplied unique query ID or source row offset. No anchors are fabricated. source_rows counts physical input rows; input_rows and sequence diagnostics count expanded chain rows, and output_rows counts raw rows returned, including missing/invalid sequences.

Parameters:

path (str | PathLike)

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

Precursor estimation is a separate API. Its optional dependencies are available with pip install 'vdjmatch[precursor]'; annotation does not run precursor estimation.

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, per_donor=False)[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_donors_hit", "prevalence"}}, "n_cognate", "n_donors", "r"}.

Per-donor, and why it is not the same number. The donor is already part of the event key, so the pooled f_hat above can be resolved per donor at no extra pass. Two quantities come out of that and they answer different questions. prevalence is the fraction of donors carrying at least one matching event – the readout a tetramer-enrichment study reports as “n of N donors responded”, and the one the saturating link prevalence = 1 - exp(-c F(e)) predicts. per_donor=True additionally returns each donor’s own f_hat under "donors", which is what a spread across donors has to be computed from: the pooled ratio is a ratio of sums and says nothing about how much donors differ. Deep donors dominate the pooled figure and are one point each in the spread, so the two disagree by construction and both are reported rather than one being derived from the other.

denominator overriding the count does not propagate to the per-donor denominators: those are always counted off the stream that was actually passed, because there is no way to split an externally supplied total across donors. A caller that pre-filters and overrides should read the per-donor numbers as conditional on that filter.

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

  • per_donor (bool)

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.load_compendium(source=None, *, measured_only=True, dedup_retabulations=False, species=None, cache=None)[source]#

Published measurements of naive precursor frequency — experimental, not computed.

147 rows from 14 studies tracing to 31 originating laboratories (one is a review, whose rows carry the laboratory that made them in lab). freq_per_1e6_naive is the harmonized column: cells per million of the naive compartment. The number as published is kept beside it in freq_per_1e6_published with the denominator it was published under, so every conversion is reversible.

Parameters:

  • source (str | PathLike | None) – path to an existing table; None fetches from the HuggingFace mirror.

  • measured_only (bool) – keep only rows that measured an endogenous naive frequency — drops the adoptive-transfer titrations, where the frequency is an experimental input rather than an observation, and the memory arms. Almost always what you want.

  • dedup_retabulations (bool) – drop the ten review rows that restate a primary already in the table (Obar 2008 on six mouse epitopes, Alanio 2010 on four human). The two copies differ only by the denominator the review renormalized to, so counting both turns one measurement into two. Off by default because a per-row median legitimately weights by how many donors stand behind a row; turn it on when you are counting measurements.

  • species (str | None) – "human" or "mouse".

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

Returns:

One row per published record. See the mirror’s README.md for the full schema.

Return type:

polars.DataFrame

Note

Two denominator traps are already applied and must not be applied twice. Adult blood is ~50% naive-phenotype, so a per-total-CD8 figure is doubled to reach the naive scale; cord blood is essentially the naive compartment in full and is not doubled.

vdjmatch.precursor.load_estimates(source=None, *, species=None, chain=None, model=None, cache=None)[source]#

What F(e) predicts for the epitopes the compendium measured — derived, not observed.

One row per (species, model set, chain, epitope), at radius 1 with alpha = 0.1. The estimator columns are repertoire fractions; multiply by 1e6 to compare against measured_per_1e6, which is carried alongside with its across-study range.

Parameters:

  • source (str | PathLike | None) – path to an existing table; None fetches from the HuggingFace mirror.

  • species (str | None) – "human" or "mouse".

  • chain (str | None) – "TRA" or "TRB".

  • model (str | None) – recombination model set — "olga" or "arda".

  • cache (str | PathLike | None) – cache directory.

Return type:

polars.DataFrame

Note

Human is scored under both model sets on purpose, because mouse exists only under arda and a human-olga against mouse-arda comparison would confound species with model family. The sets are not interchangeable: olga scores 5.7% of human TRA junctions at exactly zero, where they contribute nothing to any mass and raise no error.

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

Candidate evidence#

Per-pMHC evidence; NED ranks candidates, enrichment tests the stated reference set.

NED-v1 is the size-invariant +1-control neighbour-density formula in the existing holdout scorer. It is a ranking statistic, not a posterior probability of specificity.

vdjmatch.aggregate.candidates.candidate_schema(query_id_type=polars.UInt32)[source]#
vdjmatch.aggregate.candidates.candidates(hits, reference, controls=None, *, control_size=None, score_scale=400.0, soft_v=True, match_v=False, match_j=False, control_counts=None, distance_column='edits', background_mass=False, control_totals=None, exclude_exact=False, control_sizes=None)[source]#

Reduce independent observations to candidate evidence for each stable query_id.

controls: hit table query_id/edits, one row per unique control sequence match. Counts use unique (junction,V,J) keys; study support excludes empty identifiers. control_size=None requests explicitly uncalibrated ranking. Soft V uses existing germline-loop similarity (weight 0.25 across families); unknown/missing calls are neutral.

Query normalization#

Source-aware normalization of query junctions and their original identities.

vdjmatch.io.columns.validate_counts(df)[source]#

Reject missing, malformed, fractional, negative or overflowing supplied counts.

Parameters:

df (polars.DataFrame)

Return type:

None

vdjmatch.io.columns.normalize_query(df, valid_aa=True, *, source='auto', sequence_convention=None, return_report=False)[source]#

Normalize junctions, preserving supplied query/sequence identities and row order.

source='airr' requires junction_aa. source='custom' requires an explicit sequence_convention='junction' for ambiguous sequence columns. Auto accepts legacy short-form producers, but recognizes AIRR v_call or sequence_id headers. Bare IMGT CDR3s cannot be converted without anchors. return_report returns (frame, counts) with retained/dropped row counts.

Parameters:

  • df (polars.DataFrame)

  • valid_aa (bool)

  • source (str)

  • sequence_convention (str | None)

  • return_report (bool)

Legacy paired first-hit helpers#

These helpers retain their adaptive first-hit contract separately from the Annotator’s fixed-ball paired-reference interface. paired_scan returns (hits, alpha_control_costs, beta_control_costs). By default each hit is the historical two-tuple (radius, epitope) and pvalue counts observation associations; its target size N must use that same unit.

With include_identity=True, a hit is the three-tuple (radius, epitope, (alpha_junction, beta_junction)). pvalue then counts distinct junction pairs at the nearest paired radius; pass a distinct-pair target size N. When selecting an epitope, compute N for that epitope’s reference set. Mixing two-tuple and three-tuple hits is rejected.

Paired first-hit E-values under the existing independent-chain background model.

The nearest radius is max(alpha cost, beta cost). Historical scans count paired reference observations; include_identity=True additionally exposes paired junction identity so callers can count distinct pairs and use the same unit for target size N.

vdjmatch.evalue.paired.build_paired_ref(df)[source]#

Validate TRA/TRB linkage and return one row per complete paired observation.

Integer and rich string IDs work. Missing partners remain unpaired; multiple distinct chains at one locus or inconsistent peptide/MHC annotations raise. Exact duplicated chain rows do not create extra paired observations.

Parameters:

df (polars.DataFrame)

Return type:

polars.DataFrame

vdjmatch.evalue.paired.paired_scan(ref, control_a, control_b, pairs, *, params=None, threads=0, exclude_exact=False, include_identity=False)[source]#

One batched search per chain/reference/control.

Default hits are historical (radius, epitope) observation associations. include_identity emits (radius, epitope, (alpha, beta)) for distinct-pair counting by pvalue(); pass a distinct-pair target size N in that mode.

Parameters:

  • ref (pl.DataFrame)

  • control_a (Index)

  • control_b (Index)

  • params (SearchParams | None)

  • threads (int)

  • exclude_exact (bool)

  • include_identity (bool)

vdjmatch.evalue.paired.pvalue(paired_hits, ctrl_a_costs, ctrl_b_costs, N, Ma, Mb, epitope=None)[source]#

Conditional first-hit calibration at the nearest paired radius.

N uses observation associations for historical two-tuples, or distinct paired junctions for identity-bearing three-tuples; when restricting epitope, use its N. An empty control errors. Zero neighbours use 3/M separately for each chain.

Parameters:
Return type:

dict

Experimental local-rank evidence#

vdjmatch.evalue.local_rank.local_rank_evidence(reference_counts, control_counts, reference_population, control_population, *, prior=0.5)[source]#

Reduce batched cumulative distance balls to finite-control rank evidence.

Rows are queries; columns are identical increasing, fixed distance cutoffs. Distances between cutoffs are tied in the same bin on both sides. References beyond the last cutoff contribute their null-conditional mean rank evidence. No target density, minimum target size, random tie breaking or fitted scale is needed. Empty populations return neutral evidence and the declared prior; callers must retain an explicit unavailable-population disposition.

rank_bayes_factor is the mean normalized reciprocal background rank. p_rank_bound=min(1,1/BF) is a conservative bound, not an exact tail. posterior_rank_signal is the rank-tilted model posterior, not Prob(TP). The model conditions on fixed pooled distances and exchangeable labels; generated-draw/unique-reference differences do not establish that assumption.

Curated source replication#

vdjmatch.db.replicated() retains source provenance and complete pMHC identity when selecting corroborated records.

Experimental V-conditioned motif evidence#

vdjmatch.match.motifs.vdjbet_motif_scores(queries, reference, *, exclude_exact=False, target_reference=None)[source]#

Return supported V+masked-4-mer scores per query/full pMHC.

Inputs use full junction cdr3 and v; queries have unique query_id, reference has PMHC columns. Call separately for each species and locus. Follow vdjbet: IMGT classes, all four masks, distinct junction presence, x_length,V,target>=2, x_V,target>=3, x_target>=5, and sum(max(MF_V×IKF)) per motif divided by query length. Benchmark exclusion removes a full query junction from every target/V stratum and global counts before any gate. Missing/unsupported motifs emit no row, never a prior-only specificity. Optional target_reference defines compatible-restriction unions labelled with each query’s declared PMHC; global motif frequencies still use reference.

Targeted kernel evidence#

Uncalibrated kernel-ladder evidence over cumulative distance-ball counts.

vdjmatch.evalue.kernel.ladder_cutoffs(matrix_scale, *, radius_factor=22.5, reach_factor=100.0)[source]#

Recorded 90 inner/62 outer intervals, in ordinary matrix cost units.

vdjmatch.evalue.kernel.kernel_ladder(reference_counts, same_v_counts, control_counts, reference_sizes, control_sizes, cutoffs, *, kernel_scale, cross_v=0.25, floor=0.01, scorer='v-bounded', extras=True)[source]#

Recorded V-weighted bounded score plus five advanced ladder alternatives.

Each curve counts all eligible references through each ordered cutoff. Exact exclusions must also puncture population sizes. Scores are relative evidence, not probabilities, E-values or calibrated specificity confidence.

vdjmatch.evalue.kernel.select_kernel(queries, scores)[source]#

One ranked peptide per query, with ties/no-support retained explicitly.

Restrictions of the same peptide do not compete against one another. The support share is winner/(winner+strongest other peptide), a relative scale, not a binding probability. Original query order and metadata are retained.

vdjmatch.match.tcrdist.precompute_v_loop_distances(*, species='human', alphabet='full', matrix='tcrdist', v_regions='cdr12plus')[source]#

Exact aligned V-loop costs for reuse within one annotation invocation.

The small allele tables are immutable inputs to later population batches; no sample predictions, process-global state or persisted cache are retained. TRA/TRB have separate tables because their aligned widths can differ. v_regions="cdr12" omits CDR2.5; the default cdr12plus retains it.

Legacy calculation reproduction#

Accepted legacy matrix/LOO calculations, migrated from bench; not annotation confidence.

Legacy contracts deliberately retain unique-junction eligibility, composition caps, first-epitope assignment and balanced trapezoidal PR area rather than modern AP.

vdjmatch.match.scoring_audit.matrix_audit(frame, *, min_epi=30, vlen=4, jlen=6, legacy, threads=1)[source]#

Flat positional-matrix scores and offdiagonal Pearson correlations.

vdjmatch.match.scoring_audit.loo_audit(frame, *, locus, retention, significance, structural, min_epi=30, top=8, max_queries=600, subs=2, seed=0, threads=1)[source]#

Legacy epitope-withheld matrix training and balanced PR candidate retrieval.

On this page
Show Source