API reference #

vdjtools.model.io.from_germline(germline, *, locus, organism='custom', palindrome_max=None, ins_max=40, strict=True, source=None)[source]#

Build a recombination Model scaffold from any V(D)J germline library.

This is the entry point for a custom reference: supply your own alleles (from vdjtools.model.reference.read_germline_fasta(), from arda, or assembled by hand) and get back a model with the right gene set, germline geometry and event graph, whose marginal tables are placeholders — uniform gene usage, small-trim / small-insert biased deletions and insertions with wide support — meant to be refit by vdjtools.model.infer.infer_native(). The placeholders’ ndel / length support ranges bound what EM can learn, so they are sized to each segment’s full cut-segment length and to ins_max.

Parameters:

germline (DataFrame) – Germline frame in vdjtools.model.reference.load_germline()’s schema — required allele, segment, sequence, optional gene, functional, cdr3_anchor, full_germline. sequence is the CDR3-region germline for V/J and the full germline for D. Alleles with non-ACGT bases are dropped (the native DP cannot encode them); strict decides whether that and the other audit findings are fatal.
locus (str) – e.g. "TRB", "IGH".
organism (str) – Free-text organism tag carried in the manifest.
palindrome_max (dict[str, int] | None) – Max palindromic nt per trimmable end (default matches OLGA human).
ins_max (int) – Maximum N-region insertion length in the placeholder tables.
strict (bool) – Raise if vdjtools.model.reference.validate_germline() reports an error. Warnings (a misplaced-looking anchor, ambiguous bases, a gene-level allele name) are never fatal — call validate_germline yourself to inspect them.
source (str | None) – Manifest provenance string; defaults to "germline:<locus>".

Returns:

A validated Model — VDJ if the frame carries any D allele, else VJ.

Raises:

ValueError – If the germline frame is unusable (missing columns, no V or no J, duplicate or empty alleles) and strict is set.

Return type:

vdjtools.model.io.from_arda(locus, organism='human', *, palindrome_max=None, ins_max=40)[source]#

Build a recombination Model whose gene set + germline come from arda.

A thin wrapper over from_germline() that sources the germline library from arda — the single source of germline truth — so generated sequences carry arda’s IMGT allele names and germline, and stitching gets the full-length V/J germline + anchor from arda’s scaffolds. The marginal tables are placeholders meant to be refit; see from_germline().

Parameters:

locus (str) – e.g. "TRB", "IGH".
organism (str) – e.g. "human", "mouse".
palindrome_max (dict[str, int] | None) – Max palindromic nt per trimmable end (default matches OLGA human).
ins_max (int) – Maximum N-region insertion length in the placeholder tables.

Returns:

A validated Model (VDJ if arda has D germline for the locus, else VJ).

Raises:

ImportError – If arda is not importable (it is a base dependency; a plain pip install vdjtools ships it).
ValueError – If arda has no germline for locus / organism.

Return type:

vdjtools.model.io.marginals_frame(model)[source]#

Flatten every marginal into ONE long, self-describing table.

The per-event tables are already tidy polars, but they have different schemas, so inspecting a whole model means juggling a dozen frames. This stacks them into a single frame with the union of all realization columns (null where a column does not apply to that event), which is what you want for eyeballing, diffing, spreadsheeting or shipping a model as one file.

Parameters:: model (Model) – The model to flatten.
Returns:: event, kind, given + _REALIZATION_COLS + p, one row per probability.
Return type:: DataFrame

Example

>>> marginals_frame(m).filter(pl.col("event") == "v_choice").head()

vdjtools.model.io.set_marginals(model, frame, *, validate=True)[source]#

Rebuild a model’s marginals from a marginals_frame()-shaped table.

The inverse of marginals_frame(), and the reason a hand-edited TSV is a first-class model input: edit the probabilities in whatever tool you like, read the table back, and get a model. The germline frames and the event graph come from model — this replaces probabilities only.

Parameters:

model (Model) – Supplies the manifest and germline; its own tables are replaced.
frame (DataFrame) – Long-format marginals with at least event and p, plus whichever realization columns the events need. Extra columns are ignored; rows for unknown events are an error.
validate (bool) – Run validate() on the result.

Returns:

A new Model; the input is not modified.

Raises:

ValueError – If an event has no rows, an event name is unknown, or a required realization column is missing for some event.

Return type:

vdjtools.model.io.save_model(model, path, *, fmt='parquet')[source]#

Write a model to path/ as manifest.json + one file per event/germline table.

Parameters:

model (Model) – The model to write.
path (str | Path) – Destination directory (created if absent).
fmt (str) – "parquet" (default, exact and compact), or "tsv" / "csv" for a hand-editable directory. Text formats lose dtypes on the way out; load_model() restores them from the manifest, so the round-trip is still exact.

Raises:

ValueError – On an unknown fmt.

Return type:

None

vdjtools.model.io.load_model(path, *, validate=False)[source]#

Load a native model directory written by save_model().

The on-disk format is detected per table, so parquet, TSV and CSV directories all load, and text-format tables are cast back to the schema’s dtypes. training.json is read when present; a model written before training logs existed simply has training is None.

Parameters:

path (str | Path) – The model directory.
validate (bool) – Also run validate(). Off by default so a deliberately-broken model can still be loaded for check_model() to diagnose.

Returns:

The loaded Model.

Return type:

Precomputed recombination models shipped with vdjtools.

Three model sets live under vdjtools/model/_bundled/:

olga — imported from OLGA’s default models (the exact-Pgen bootstrap; single-D). Seven human loci, OLGA’s germline namespace.
learned — EM-inferred from real non-functional reads (out-of-frame + stop-codon, HuggingFace), tandem-D on the D-bearing loci (IGH/TRD/TRB). These carry a learned P(n_D=2) and broader trim/insertion distributions than the synthetic OLGA models. Seven human loci, OLGA’s germline namespace.
arda — the same EM inference on the arda IMGT allele namespace (from_arda() scaffold refit by infer_native()). Nine models: the seven human loci plus mouse TRA/TRB — the only bundled set with a non-human organism. Use this one when the rest of your pipeline is arda-annotated, so generated sequences share one allele namespace with your query data instead of needing a name fallback.

arda models are keyed by {organism}_{LOCUS} on disk; the other two by {LOCUS} alone. load_bundled() hides that — pass organism= and it picks the right key.

Each model is a directory of parquet marginal tables + manifest.json (see vdjtools.model). Provenance and the build command are recorded in SOURCES.md.

vdjtools.model.bundled.SOURCES = ('olga', 'learned', 'arda')#: The three shipped model sets.

vdjtools.model.bundled.LOCI = ('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL')#: The seven human loci with bundled models (arda adds mouse TRA/TRB).

vdjtools.model.bundled.load_bundled(locus, source='olga', *, organism='human', collapse=True)[source]#

Load a precomputed model shipped with the package (no OLGA/HuggingFace at runtime).

Parameters:

locus (str) – One of TRA TRB TRG TRD IGH IGK IGL (case-insensitive).
source (str) – "olga" (OLGA bootstrap, exact Pgen), "learned" (EM-inferred from real non-functional reads; tandem-D on IGH/TRD/TRB), or "arda" (the same EM inference on the arda IMGT allele namespace, and the only set covering mouse).
organism (str) – "human" (default) or "mouse". Only the arda set ships a non-human organism; asking the others for one is an error rather than a silent human model.
collapse (bool) – If True (default), collapse each gene to a single *01 allele via collapse_alleles() — the working gene-level resolution, in which Pgen also collapses a clonotype’s allele to *01 (short-read aligners cannot resolve alleles reliably, so the suffix is noise). Pass collapse=False for full allele resolution, e.g. the exact-OLGA-Pgen fidelity check.

Returns:

The Model.

Raises:

ValueError – If source is not one of SOURCES, or a non-human organism is asked of a human-only set.
FileNotFoundError – If no bundled model exists for (source, organism, locus).

Return type:

Example

>>> m = load_bundled("TRB")                                  # olga, human
>>> m_arda = load_bundled("TRB", "arda")                     # arda namespace, human
>>> m_mouse = load_bundled("TRB", "arda", organism="mouse")  # the only mouse models

vdjtools.model.bundled.list_bundled()[source]#

Return the available bundled models as {source: [key, ...]}.

Keys are bare loci for olga/learned and {organism}_{LOCUS} for arda, matching what load_bundled() accepts for each set.

Return type:: dict[str, list[str]]

Collapse a V(D)J model to one allele per gene — the default working resolution.

IMGT alleles of one gene (TRBV12-4*01, TRBV12-4*02) differ by a handful of SNPs, usually in framework rather than the CDR3 region, and short-read aligners cannot reliably tell them apart (arda splits usage across *01/*03 as mapping noise). So the default working model sums each gene’s sub-variants into a single *01 representative, and Pgen collapses a clonotype’s allele to that representative too — the allele suffix stops carrying (spurious) information. Pass collapse=False (or keep the un-collapsed model) when you genuinely want allele resolution, e.g. the exact-OLGA-Pgen fidelity check.

The collapse is a proper marginalisation, not a truncation:

a choice probability (P(V=a)) becomes the plain sum over the gene’s alleles — exact;
a conditional (P(delV | V=a), P(D | J=ja)) becomes the usage-weighted average over the gene’s alleles, P(x | gene) = Σ_a [P(a)/P(gene)] · P(x | a) — the correct allele-marginal;
the germline of one representative allele is kept and relabelled gene*01. This is the one lossy step: where alleles differ inside the CDR3 region the mixture cannot be reproduced by a single germline, so collapsed Pgen is approximate there (exact wherever the CDR3-region germline is allele-invariant, which is the common case).

The representative is ranked, in order, by CDR3-region germline length, then IMGT functionality (F > ORF > P, read from arda’s per-allele CDR3-anchor table), then usage, then a preference for *01, then the allele name as a last resort.

Length leads because IMGT ships some alleles with a truncated — sometimes empty — CDR3-region germline, and a truncated allele must never define the gene’s trim range. Human IGKV3-20 is the usage version of that trap: *02 is 11 nt against *01’s 30 and had the higher learned usage, so the collapsed gene inherited an 11-nt germline (labelled *01, which was doubly misleading) and stranded 25% of its own deletion distribution on trims it could no longer reach. Functionality cannot lead for the same reason: over every bundled model the two orders differ on exactly one gene, human TRBV23/OR9-2, where the only non-pseudogene allele (*02, ORF) has an empty CDR3 germline and *01 (P) has 21 nt — leading with functionality would install the empty one and make Pgen through that gene exactly 0. Alleles of one gene are near-identical through the CDR3 region, so a large length gap means an incomplete database entry, not biology.

Functionality is second because ties on length are common and the name tie-break underneath them is meaningless — it just takes the lexicographically last allele. That is how TRBJ2-7*02, an ORF allele whose germline templates SYEQYV instead of the conserved SYEQYF, became the representative of TRBJ2-7 and was relabelled *01: both alleles are 19 nt, so length did not separate them, no usage was passed for J at all, and *02 won on the name. Every real …YEQYF junction then scored exactly zero against the collapsed model, with no error raised, and the TRBJ2-7 marginal fell 5–6 orders of magnitude.

Having fixed the representative, the collapsed deletion conditionals are then projected onto what that germline supports and renormalized: the collapsed model is a single-germline model, so any residual mass on unreachable trims would leak straight out of Pgen.

vdjtools.model.collapse.collapse_alleles(model)[source]#

Return a copy of model with every gene reduced to a single *01 allele.

Parameters:: model (Model) – A V(D)J Model (OLGA-derived or EM-learned).
Returns:: one allele (gene*01) per gene in every table and in the germline. Marginal usage is preserved exactly; conditionals are the correct usage-weighted allele averages; germline is the representative allele’s (ranked by length, then IMGT functionality, then usage, then *01), relabelled *01.
Return type:: A new Model at gene resolution

Example

>>> m = collapse_alleles(load_bundled("TRB", "learned"))
>>> m.tables["v_choice"].height        # one row per V gene, not per allele

Germline reference — arda is the single source of germline truth.

Every V/D/J germline sequence and CDR3 anchor used anywhere in vdjtools (annotation, scenario enumeration, contig stitching, generation, and arda-native models) resolves from arda’s germline library by allele name, so the whole pipeline speaks one coordinate frame.

arda’s anchor convention is byte-identical to OLGA’s: anchor_nt is a 0-based offset into the full germline marking the conserved Cys104 codon (V) or [FW]118 codon (J); the CDR3-region germline is full[anchor:] for V and full[:anchor+3] for J. So no coordinate conversion is needed between the two — this module documents and enforces that shared frame (reconcile_olga() catalogs any residual sequence differences).

Full-length V/J germline for contig stitching (Phase 1c) is recovered from arda’s bundled scaffold reference by load_full_vj_germline() / arda_full_germline(): arda’s alleles.fasta scaffolds are full_V + N-pad + full_J, sliced per allele at the v_sequence_end / j_sequence_start boundaries from arda.annotate.reference. The sliced full germline is anchor-consistent with load_germline() (verified across all functional V/J alleles: full_V[anchor:] == the CDR3-region germline, and the CDR3-region J germline is a prefix of full_J). Only functional/ORF alleles with complete markup have a full germline in arda; pseudogenes (present in the CDR3 anchors) may be absent.

vdjtools.model.reference.reverse_complement(seq)[source]#

Reverse complement of a nucleotide string (ACGT).

Parameters:: seq (str)
Return type:: str

vdjtools.model.reference.translate(seq)[source]#

Translate a nucleotide string to amino acids (standard code; trailing partial codon dropped).

Parameters:: seq (str)
Return type:: str

vdjtools.model.reference.cut_segment(seq, segment, max_pal)[source]#

Palindrome-extend a CDR3-region germline for the Pgen DP (mirrors OLGA’s cutR/cutL_seq).

Appends up to max_pal reverse-complement (P-nucleotide) bases on the trimmable end so a deletion index into the result directly counts nt removed:

V: append at the 3’ end. - J: prepend at the 5’ end.

Parameters:

seq (str) – CDR3-region germline (V: Cys104→3’ end; J: 5’→[FW]118 codon end).
segment (str) – "V" or "J".
max_pal (int) – Maximum palindromic nt for this end.

Return type:

vdjtools.model.reference.cut_segment_d(seq, max_pal5, max_pal3)[source]#

Palindrome-extend a D germline on both ends (5’ then 3’).

Parameters:

seq (str)
max_pal5 (int)
max_pal3 (int)

Return type:

vdjtools.model.reference.load_germline(locus, organism='human')[source]#

V/D/J germline + CDR3 anchors for a locus, from arda (the source of truth).

Parameters:

locus (str) – e.g. "TRB", "TRA", "IGH".
organism (str) – e.g. "human", "mouse".

Returns:

One row per allele with columns allele, gene, segment, sequence, cdr3_anchor, functionality, functional, status. For V/J, sequence is the CDR3-region germline and cdr3_anchor the 0-based anchor codon offset in the full germline; for D, sequence is the full D germline, cdr3_anchor = -1.

Raises:

ImportError – If arda is not importable (it is a base dependency; a plain pip install vdjtools ships it).
ValueError – If no germline is found for locus / organism.

Return type:

DataFrame

vdjtools.model.reference.load_full_vj_germline(organism='human')[source]#

Full-length V and J germline nucleotide sequences from arda, by (segment, allele).

arda ships full V/J germline only inside deduplicated V–J scaffolds (database/vdj/<organism>/alleles.fasta, keyed by opaque scaffold id); this slices each scaffold at the v_sequence_end / j_sequence_start boundaries from arda.annotate.reference.load_reference to recover the per-allele full germline (V: FR1 → 3’ end of V-REGION; J: 5’ J → end of FR4). The first scaffold carrying an allele wins (all are byte-identical for that allele’s segment).

Parameters:: organism (str) – e.g. "human", "mouse".
Returns:: {("V"|"J", allele): full_germline_nt}. Only functional/ORF alleles with complete arda markup are present; pseudogenes (in the CDR3 anchors) may be missing.
Raises:: ImportError – If arda is not importable (it is a base dependency; a plain pip install vdjtools ships it).
Return type:: dict[tuple[str, str], str]

vdjtools.model.reference.arda_full_germline(locus, organism='human')[source]#

Stitch-ready full V/J germline + anchor for a locus, entirely from arda.

Combines load_full_vj_germline() (full-length germline) with load_germline() (the CDR3-region germline) so the anchor is derived self-consistently by length — no reliance on a cross-source coordinate assumption:

V: anchor = len(full) - len(cdr3_region) — full[:anchor] is the framework 5’ of the conserved Cys104.
J: anchor = len(cdr3_region) - 3 — full[anchor + 3:] is the framework 3’ of the conserved Phe/Trp118 codon.

This is exactly the (full_germline, anchor) pair vdjtools.model.stitch.stitch_contig() consumes, so an arda-native model (no OLGA germline) can stitch full contigs. Alleles whose full germline is absent from arda (pseudogenes / incomplete markup) are skipped.

Parameters:

locus (str) – e.g. "TRB", "IGH".
organism (str) – e.g. "human".

Returns:

{("V"|"J", allele): (full_germline_nt, anchor)}.

Return type:

dict[tuple[str, str], tuple[str, int]]

vdjtools.model.reference.GERMLINE_REQUIRED = ('allele', 'segment', 'sequence')#: Columns a germline frame must carry, and the defaults filled in for the optional ones.

vdjtools.model.reference.normalize_germline(germline)[source]#

Fill a germline frame’s optional columns with their defaults; returns a new frame.

gene defaults to allele.split("*")[0]. See GERMLINE_OPTIONAL for the rest.

Parameters:: germline (DataFrame)
Return type:: DataFrame

vdjtools.model.reference.validate_germline(germline)[source]#

Audit a germline frame destined for vdjtools.model.io.from_germline().

Catches the mistakes that otherwise produce a model that builds cleanly and scores wrongly — above all a misplaced CDR3 anchor, which shifts every deletion profile by a constant and is invisible downstream.

Parameters:: germline (DataFrame) – A frame with at least allele, segment, sequence (see GERMLINE_REQUIRED); optional columns are described by GERMLINE_OPTIONAL. sequence is the CDR3-region germline (V: Cys104 codon → 3’ end; J: 5’ end → through the [FW]118 codon), or the full germline for D.
Returns:: Tidy issue frame severity, check, event, segment, allele, detail, value; empty when the library is clean. severity == "error" means the frame cannot build a model.
Return type:: DataFrame

vdjtools.model.reference.forbidden_dj_pairs(d_alleles, j_alleles, locus)[source]#

(D, J) pairs that deletional V(D)J recombination cannot produce, for this locus.

Recombination excises the DNA between the two segments it joins, so a D can only be joined to a J that lies 3’ of it. In almost every locus that is vacuous — all D genes sit 5’ of all J genes — but TRB interleaves:

5’ ─ TRBV… ─ TRBD1 ─ TRBJ1 cluster ─ TRBD2 ─ TRBJ2 cluster ─ TRBC2 ─ 3’

so TRBD2 can only reach the TRBJ2 cluster; a TRBD2–TRBJ1 join is physically impossible. TRBD1, sitting 5’ of both clusters, can reach either.

This is a hard genomic fact, not a probability the data gets to inform, and EM will happily learn the impossible pair from noise if nothing forbids it: a D match is only 10–18 nt, which is weak evidence on short reads. Before this constraint existed the retrained human TRB model assigned 0.091 to that pair, and OLGA’s own model assigns 0.333. Ambiguous d_call values in the training data concentrate exactly there, so the error is systematic, not random.

Parameters:

d_alleles – D allele names, e.g. ["TRBD1*01", "TRBD2*01"].
j_alleles – J allele names.
locus (str) – e.g. "TRB". Loci with no interleaving return an empty set.

Returns:

The set of impossible (d_allele, j_allele) pairs.

Return type:

set[tuple[str, str]]

Example

>>> forbidden_dj_pairs(["TRBD2*01"], ["TRBJ1-1*01", "TRBJ2-1*01"], "TRB")
{('TRBD2*01', 'TRBJ1-1*01')}

vdjtools.model.reference.read_fasta(path)[source]#

Read a FASTA into (header, sequence) pairs, transparently handling .gz.

Parsing itself is delegated to arda (one FASTA parser in the stack, not two), but arda’s reader takes a path and plain open(), so a gzipped file reaches it as mojibake and dies on the first byte. Gzipped input is therefore decompressed to a temporary file first.

Return type:: list[tuple[str, str]]

vdjtools.model.reference.read_germline_fasta(v, j, d=None, *, anchors=None)[source]#

Build a germline frame from your own FASTA files — the entry point for a custom library.

The segment comes from which argument a file was passed as, so no header convention is assumed: a header is either >ALLELE or >ANYTHING|ALLELE (arda’s D convention), and everything after the second | is ignored.

Parameters:

v – FASTA path for the V alleles.
j – FASTA path for the J alleles.
d – Optional FASTA path for the D alleles. Supplying it makes the model VDJ; omitting it makes it VJ.
anchors – Optional CDR3-anchor CSV in OLGA’s *_gene_CDR3_anchors.csv format (gene,anchor_index,function). When given, the FASTAs are taken to hold full-length germline and are sliced to the CDR3 region (full[anchor:] for V, full[:anchor + 3] for J). When omitted, the V/J sequences are taken to be CDR3-region germline already — validate_germline() flags it if they are not.

Returns:

A germline frame in load_germline()’s schema, ready for vdjtools.model.io.from_germline().

Return type:

DataFrame

vdjtools.model.reference.reconcile_olga(model)[source]#

Catalog how an OLGA-loaded model’s germline relates to arda’s (the shared-frame audit).

OLGA bootstrap models keep OLGA’s germline geometry for exact-Pgen fidelity; this reports, per V/J allele, whether it resolves in arda and whether the CDR3-region germline matches — so divergences (IMGT-version drift, V 3’ extent) are flagged, never silent.

Parameters:: model – A Model loaded via from_olga.
Returns:: allele, segment, in_arda, germline_equal, olga_len, arda_len.
Return type:: Per-allele report

Reconstruct full-length nt contigs from (V, J, CDR3) and run them through arda.

OLGA generation (and OLGA-synthetic bootstrap data) emits only (V call, J call, CDR3 nt), but arda — the aligner that drives scenario enumeration for real reads — consumes full nt reads. stitch_contig rebuilds the read by prepending the V germline 5’ of the CDR3 anchor and appending the J germline 3’ of it, so synthetic and real reads flow through the identical arda.annotate_sequences → scenarios → EM path (with known ground truth on the synthetic side).

vdjtools.model.stitch.stitch_contig(model, v, j, cdr3_nt)[source]#

Rebuild a full nt contig: V framework 5’ of Cys104 + CDR3 + J framework 3’ of [FW]118.

Parameters:

model (Model) – Model providing full germline + anchors (genes_v/genes_j).
v (str) – V and J allele names.
j (str) – V and J allele names.
cdr3_nt (str) – The junction-space CDR3 nt (Cys104 → [FW]118 inclusive, i.e. AIRR junction).

Returns:

The contig, or None if either gene lacks a usable anchor / full germline.

Return type:

str | None

vdjtools.model.stitch.stitch_frame(model, gen, *, cdr3_col='junction_nt')[source]#

Add a contig column to a generated frame (rows with no anchor drop to null).

Parameters:

model (Model)
gen (DataFrame)
cdr3_col (str)

Return type:

DataFrame

vdjtools.model.stitch.annotate(contigs, *, organism='human')[source]#

Annotate nt contigs with arda → a polars frame of the calls the scenario/EM path needs.

Returns columns v_call, d_call, j_call, junction, junction_aa, productive (arda’s best alignment). arda is a base dependency (ships with vdjtools).

Parameters:

contigs (list[str])
organism (str)

Return type:

DataFrame

Most-likely nucleotide CDR3 for an amino-acid CDR3, and the recombination scenario behind it.

VDJdb-style records carry (V, J, CDR3aa) and no nucleotides. infer_nt() returns the nucleotide CDR3 that maximises the model’s generation probability, together with the argmax scenario — which is the V/D/J boundary markup (v_end, d_start, d_end, j_start).

This is the argmax counterpart of a sum that already exists. pgen.pgen_aa() marginalises over every nucleotide sequence encoding a given amino-acid CDR3; the most likely such sequence is the maximum of the same set, and pgen.pgen_nt()’s own loops — (V, len_v), (J, len_j), then the D/insertion middle — already enumerate the scenarios. Nothing here re-derives the recombination model: every probability comes from pgen.prepare()’s tables.

Two objectives are deliberately separated, because they are not the same number:

infer_nt() maximises the marginal P(nt) = pgen_nt(nt), which sums over every scenario that could produce that sequence. That is the right target for “the most likely nucleotide sequence”.
best_scenario() then maximises over scenarios for that fixed sequence, which is what the boundary markup has to come from. The two maxima are over different things and taking the second as an answer to the first would be wrong.

STATUS.

best_scenario() is complete and validated. It needs no enumeration: it is a max-product walk of the (V, len_v) x (J, len_j) x middle loops pgen.pgen_nt() already sums over. Checked on 500 generated draws per locus: the reported V span is the V germline 500/500, the J span the J germline 500/500, the D span the D germline 500/500 (TRB), and scenario_p <= pgen_nt 500/500 — a max cannot exceed the sum it is taken over.
infer_nt() is the production path: a codon-constrained max-product DP over the same scenario enumeration pgen.pgen_aa() sums over, with exact branch-and-bound, followed by an exact pgen.pgen_nt() re-score of the top candidates.
infer_nt_bruteforce() is exact but exponential, and exists as the TEST ORACLE.

Why brute force cannot be the answer, measured on the real VDJdb (79,997 records). The codon search space — the product of synonymous-codon counts over the CDR3 — has a median of 5.3e6 for TRA and 1.9e7 for TRB; only 8.9 % / 1.6 % of records are at or under 1e5 candidates. Each candidate costs a full pgen.pgen_nt(), so enumeration is out by orders of magnitude.

WARNING: And it cannot be rescued by pinning the germline-templated flanks, which was tried and is unsound: pinning a codon to the V/J germline excludes every sequence in which that germline was trimmed, and the true maximum can be one of them. Caught against the brute-force oracle — CAVSDMRF under TRAV21*01 gives …GTGAGTGAC… where the pinned version returns …GTGAGCGAT…. The pin is gone; do not reintroduce it as an optimisation.

NOTE: Both functions assume an in-frame CDR3, i.e. len(nt) == 3 * len(aa). That is the VDJdb population (real, productive receptors) but NOT every draw from generate() — an out-of-frame rearrangement has e.g. 8 aa and 25 nt, and no in-frame nt can explain it. Test with productive_only=True or the comparison is meaningless.

class vdjtools.model.viterbi.Scenario(cdr3_nt, v_call, j_call, v_end, j_start, d_call=None, d_start=None, d_end=None, pgen=0.0, scenario_p=0.0, n_candidates=0, truncated=False, runner_up_pgen=0.0)[source]#

Bases: object

A nucleotide CDR3 and the recombination that most likely produced it.

Coordinates are 0-based, half-open, in CDR3 nucleotide space (the conserved Cys through the conserved Phe/Trp, both included) — the same space as cdr3_nt itself. d_start/d_end are None when the model has no D (a VJ chain) or no D was placed.

Parameters:

cdr3_nt (str)
v_call (str)
j_call (str)
v_end (int)
j_start (int)
d_call (str | None)
d_start (int | None)
d_end (int | None)
pgen (float)
scenario_p (float)
n_candidates (int)
truncated (bool)
runner_up_pgen (float)

cdr3_nt: str#

v_call: str#

j_call: str#

v_end: int#

j_start: int#

d_call: str | None = None#

d_start: int | None = None#

d_end: int | None = None#

pgen: float = 0.0#

scenario_p: float = 0.0#

n_candidates: int = 0#

truncated: bool = False#

runner_up_pgen: float = 0.0#

property margin: float#

pgen / runner_up_pgen. 1.0 means a tie; large means the argmax is unambiguous.

NOTE: Report this. With a long non-templated core many sequences are near-equally likely and a bare “most likely” claim is then close to meaningless — the margin is what says so.

vdjtools.model.viterbi.best_scenario(model_or_prep, cdr3_nt, v=None, j=None)[source]#

The single most likely recombination for a KNOWN nucleotide CDR3.

A max-product walk of exactly the loops pgen.pgen_nt() sums over — (V, P(V), len_v) × (J, P(J), len_j) × the middle — so the scenario space, the pruning and every probability are the shipped ones. Returns None when no (V, J) pair can explain the sequence.

NOTE: The D placement is taken from the middle by the same gapless scan the model uses; d_call is None for a VJ chain.

Parameters:

cdr3_nt (str)
v (str | None)
j (str | None)

Return type:

Scenario | None

vdjtools.model.viterbi.infer_nt(model_or_prep, cdr3_aa, v=None, j=None, *, n_best=8, keep=0.01, top_vj=8)[source]#

The most likely nucleotide CDR3 for an amino-acid CDR3, and its recombination markup.

Two stages, and the second is what makes the first honest:

Max-product over (scenario, nucleotides) jointly. Every scenario pgen.pgen_aa() sums over is enumerated — (V, delV) x (J, delJ) x (D, delD, position) — and for each, _aa_dp_max() picks the single best codon assignment: germline positions are pinned to their segment, and each free N-region position takes the nucleotide maximising P(nt_1)·prod P(nt_k | nt_{k-1}) under the VD / DJ / VJ dinucleotide model covering it.
The surviving candidates are re-scored with the exact marginal pgen.pgen_nt(), and the winner is the best of those. Stage 1 maximises the joint P(nt, scenario), which is not the marginal P(nt) this function’s contract asks for; the re-scoring is what turns a Viterbi path into an answer about the marginal, and it is why pgen on the result is a real pgen_nt rather than a path weight.

Both stages earn their place, measured against infer_nt_bruteforce() on generated productive draws with V and J pinned (25 TRG, 19 TRA, 10 residues each).

variant	TRG	TRA
one scenario, best codon per residue (no search)	9/25	4/19
the joint argmax alone, no re-score	21/25	15/19
argmax + marginal re-score (the defaults)	25/25	19/19

The cheap shortcut fails because a trim chosen before the codons pins a codon the true optimum would have trimmed away — the same unsoundness as pinning the germline flanks.

Two implementations. Pass a Model and stage 1 runs natively (vdjtools.model.native.best_aa_scenarios()) — 2.5 ms per human TRB CDR3, 0.5 ms per TRA, i.e. all 80k VDJdb records in about 3 minutes. Pass an already-pgen.prepare() -d model and it runs the pure-Python scenario search instead, which is the reference implementation the native one is validated against and is ~600x slower on TRB. keep and top_vj apply only to that reference path; the native DP needs neither, because it marginalizes over V and J at essentially no cost (0.26 vs 0.23 ms with both free).

Three input modes for the calls, because annotation tables have all three:

v="TRBV5-1*01", j="TRBJ2-3*01" — the normal mode. Exhaustive over the trims of that pair.
v="TRBV6-2*01,TRBV6-3*01" (or a list) — the aligner offered several and could not choose. Every listed allele is searched and the model picks the most plausible.
v=None — nothing known; the DP marginalizes over every V (or J).

Parameters:

model_or_prep – A Model (native, the normal call) or a pgen.prepare() -d one (the pure-Python reference search).
cdr3_aa (str) – The CDR3 amino-acid sequence (conserved Cys -> conserved Phe/Trp inclusive).
v – V allele, comma-separated alleles, a list, or None — see above. Allele-keyed, never gene-level: see the _gene_idx trap.
j – J allele, comma-separated alleles, a list, or None.
n_best (int) – How many distinct candidates to re-score with the exact marginal Pgen.
keep (float) – Reference path only. Candidates within this factor of the running best survive the prune. 1.0 is the pure argmax and the fastest.
top_vj (int) – Reference path only. Shortlist size when a call is None, ranked by the model’s own P(gene)·P(del); 0 searches everything. It exists because the Python search costs ~23 s per unconstrained TRB CDR3. The native path ignores it.

Returns:

A Scenario with pgen the marginal of the winner and runner_up_pgen the marginal of the next distinct candidate — report Scenario.margin, because with a long non-templated core many sequences are near-equally likely. None when no nucleotide CDR3 encoding cdr3_aa is explicable under the model.

Return type:

Scenario | None

Note

Tandem-D (n_D = 2) scenarios are not enumerated in stage 1. They cannot add candidates — a single D trimmed to zero length already reaches every middle — so they can only reorder them, and the stage-2 pgen_nt re-scoring includes the tandem-D mass in full. Raise n_best if a D-D model’s ordering matters to you.

vdjtools.model.viterbi.infer_nt_bruteforce(model, cdr3_aa, v=None, j=None, *, max_candidates=200000)[source]#

Exact most-likely nucleotide CDR3 by enumerating every codon assignment. Oracle only.

Maximises the marginal pgen_nt over every nucleotide sequence encoding cdr3_aa, then reports the argmax scenario for the winner via best_scenario(). Exact, and exponential in the number of residues.

WARNING: Not usable on real data — see the module header: the median VDJdb record has 5.3e6 (TRA) to 1.9e7 (TRB) candidates. This exists so the production DP, when written, has an exact reference to be tested against on the small cases where both can run. max_candidates bounds the work and sets truncated; a truncated result is the best of a PARTIAL set and is not the maximum, which is why it is flagged rather than returned quietly.

Returns None when no candidate is explicable under the model — including the out-of-frame case, where no in-frame sequence of length 3 * len(aa) exists.

Parameters:

model (Model)
cdr3_aa (str)
v (str | None)
j (str | None)
max_candidates (int)

Return type:

Scenario | None

vdjtools.model.viterbi.codon_options(aa)[source]#

Per residue, the codons encoding it. Unknown residues get no options (an empty list).

Parameters:: aa (str)
Return type:: list[list[str]]

Bootstrap the recombination model from real AIRR reads (private HF isalgo/airr_model_read).

The dataset ships raw 5’-RACE FASTQ reads per organism-group (human, human_fetal — TdT-low, mouse) and chain, split into functional (productive) and non-functional (out-of-frame or stop) buckets. Training / benchmarking the engine means:

fetch_fastq() — pull one {group}/{CHAIN}.{label}.fq.gz from the hub.
annotate_reads() — map the reads with arda (arda rnaseq map) to a per-read AIRR table with V/D/J (and D2 for D-D joins) calls, junction, CIGARs, isotype, and productivity.
unique_clonotypes() — collapse to the model’s clonotype identity: same V allele, J allele, and junction nucleotides. Reads differing only in alignment (CIGAR) or isotype (IGH c_call) are one clonotype; isotype is dropped from the key (isotype switching is the same clonotype) and the collapse can be restricted to naive IgM.

The non-functional clonotypes are the unbiased EM training set; the functional ones are the selection-shaped test set. huggingface_hub (fetch) and arda[rnaseq] (annotate; needs arda + mmseqs2 + seqtree) are lazy, tool-only imports — never runtime dependencies of the model math.

Example:

from vdjtools.model import from_olga, data
from vdjtools.model.infer import infer_native
clones = data.prepare("human", "TRB", "nonfunctional", out_dir="/tmp/arda", cap=200_000)
template = from_olga(olga_dir, locus="TRB")
fit, rep = infer_native(template, clones["junction"].to_list())

vdjtools.model.data.MODEL_READS_REPO = 'isalgo/airr_model_read'#: Private HuggingFace dataset of real AIRR reads (owner isalgo, cc-by-nc-nd); fetched, not vendored.

vdjtools.model.data.ORGANISM = {'human': 'human', 'human_fetal': 'human', 'mouse': 'mouse'}#: arda reference organism per group (fetal T cells use the human reference).

vdjtools.model.data.fetch_fastq(group, chain, label, *, repo='isalgo/airr_model_read')[source]#

Download one {group}/{CHAIN}.{label}.fq.gz from the dataset; return its local path.

Parameters:

group (str)
chain (str)
label (str)
repo (str)

Return type:

vdjtools.model.data.annotate_reads(fq_path, *, out_dir, prefix, organism='human', cap=None, threads=0, reconstruct=False, verbose=False)[source]#

Map a FASTQ to a per-read AIRR table with arda (arda rnaseq map → <prefix>.airr.tsv).

Returns the raw per-read annotation (V/D/J and D2 calls, junction, CIGARs, isotype, productivity); collapse it to clonotypes with unique_clonotypes(). We deliberately stop at map — the clonotype identity “same V/J allele + junction, up to CIGAR” is exactly that dedup, so arda’s heavier error-model correct stage is not needed (and, on these single-end 5’-RACE reads, it discards reads that map but whose mate-spanned junction it cannot reassemble).

Parameters:

fq_path (str | Path) – Input FASTQ (.fq/.fq.gz), e.g. from fetch_fastq().
out_dir (str | Path) – Directory for arda’s output.
prefix (str) – Output basename.
organism (str) – arda reference organism ("human" / "mouse").
cap (int | None) – If set, annotate only the first cap reads (bounded benchmark scale).
threads (int) – mmseqs threads (0 = all cores).
reconstruct (bool) – Merge overlapping paired mates into one fragment before mapping (recovers longer junctions single reads don’t span; needs the FASTQ to carry both mates).
verbose (bool)

Returns:

The per-read AIRR polars.DataFrame. Reads that do not span a full junction have an empty junction and are dropped by unique_clonotypes().

Return type:

DataFrame

Requires the arda CLI with the rnaseq extra (arda + mmseqs2) on PATH.

vdjtools.model.data.unique_clonotypes(clones, *, naive_igm_only=False)[source]#

Collapse an arda clonotype table to unique clonotypes = (v_call, j_call, junction).

Isotype (c_call) is dropped from the key — isotype switching is the same clonotype. Rows are keyed to allele resolution (V/J allele) and full junction nt, so reads that differ only in alignment (CIGAR) or isotype collapse together; read support sums into count.

Parameters:

clones (DataFrame) – An arda clonotype frame (from annotate_reads()).
naive_igm_only (bool) – Keep only IgM (naive B) clonotypes before collapsing (IGH only; no-op elsewhere) — for a naive-repertoire (pre-selection-independent) subset.

Returns:

v_call, j_call, junction, junction_aa, locus, d_call, d2_call, count (one row per unique clonotype), sorted by descending count.

Return type:

Deduplicated clonotype frame

vdjtools.model.data.prepare(group, chain, label, *, out_dir, cap=None, reconstruct=False, naive_igm_only=False, verbose=False)[source]#

Fetch → map (arda) → unique clonotypes for one (group, chain, label) bucket.

Parameters:

group (str)
chain (str)
label (str)
out_dir (str | Path)
cap (int | None)
reconstruct (bool)
naive_igm_only (bool)
verbose (bool)

Return type:

DataFrame

vdjtools.model.data.PREPARED_DIR = PosixPath('/home/runner/work/vdjtools/vdjtools/tests/python/fixtures/model_reads')#: Directory of the arda-mapped clonotype examples shipped with the source tree.

vdjtools.model.data.write_prepared(clones, path)[source]#

Write a clonotype frame as the gzipped FASTA load_prepared() reads.

FASTA rather than a table because the payload is a sequence set — it stays greppable, readable by any bioinformatics tool, and about half the size of the equivalent TSV. The V/J/D calls EM needs for its per-read masks ride in the header, pipe-separated:

>{id}|{v_call}|{j_call}|{d_call}|{d2_call}|{count} then the junction nucleotides.

Empty fields are written as empty strings and read back as nulls.

Parameters:

clones (DataFrame)
path (str | Path)

Return type:

vdjtools.model.data.load_prepared(group='human', chain='TRB', label='nonfunctional', *, path=None)[source]#

Load a small arda-mapped clonotype example — no arda, no mmseqs2, no network, seconds.

prepare() runs the full pipeline on the raw FASTQ, which needs arda + mmseqs2 and takes minutes per chain. For examples, notebooks, tests and a quick look at real junctions that is overkill, so a few already-mapped subsets ship with the source tree as gzipped FASTA (see write_prepared() for the header format).

These live under tests/ and are not packaged into the wheel — they exist for working in a checkout. From an installed vdjtools, use prepare() or pass path=.

Parameters:

group (str) – human, human_fetal or mouse.
chain (str) – e.g. "TRB".
label (str) – "nonfunctional" (the EM training bucket) or "functional".
path (str | Path | None) – Read this file instead of looking one up.

Returns:

junction, v_call, j_call, d_call, d2_call, count — the columns unique_clonotypes() produces, minus the ones a FASTA cannot carry.

Raises:

FileNotFoundError – If no example ships for that combination. Only a couple do; build the rest with prepare() (or the whole corpus with build_all()).

Return type:

DataFrame

Example

>>> clones = load_prepared("human", "TRB")
>>> model, report = infer_frame("TRB", clones, max_iter=10)

vdjtools.model.data.BUILD_DEFAULTS = {'gene_prior': 1.0, 'iters': 15, 'nd_prior': 0.0, 'tol': 0.0001}#

P(V) = 0 is an absorbing state of this EM, so without a Dirichlet pseudocount over the germline’s functional alleles one unlucky iteration deletes a real gene permanently (human TRB kept 30 of 57 V genes unregularized, having seen 54 in the data).

Type:: EM defaults for a corpus build. gene_prior is not optional in spirit

vdjtools.model.data.build_model(chain, *, group='human', template=None, clones=None, work_dir='/tmp/vdjtools_build', cap=None, iters=15, tol=0.0001, single_d=False, nd_prior=0.0, gene_prior=1.0, threads=0, ambiguous='A', verbose=False, checkpoint=None, checkpoint_every=1)[source]#

Fetch, annotate and EM-fit a model for one chain — the whole corpus pipeline, end to end.

Parameters:

chain (str) – e.g. "TRB", "IGH".
group (str) – human, human_fetal or mouse. Sets the arda reference organism.
template – Model supplying the gene set, germline and event graph. Defaults to from_arda() for the chain and the group’s organism.
clones (DataFrame | None) – Skip fetch+annotate and fit these clonotypes instead (an EM_DATA_DIR parquet, a load_prepared() subset, or your own reads). Must carry junction, v_call, j_call and, for a D locus, d_call/d2_call.
work_dir (str | Path) – Scratch directory for arda’s per-read output (it caches there).
cap (int | None) – Read cap for annotation. None uses every read — a subsample is a silent claim that the tail does not matter, and the tail is where the rare V genes are.
iters (int) – EM iteration cap. Convergence is on relative log-likelihood, so this is a safety net.
tol (float) – Stop when the relative log-likelihood improvement falls below this.
single_d (bool) – Force a strict single-D model on a D-bearing locus.
nd_prior (float) – Dirichlet pseudocount pushing P(n_D=2) toward 0.
gene_prior (float) – Dirichlet pseudocount over the germline’s functional V/J alleles — see BUILD_DEFAULTS.
threads (int) – E-step worker threads (0 = auto).
ambiguous (str | None) – Substitute this base for any non-ACGT character in a junction (default "A"), or None to drop those clonotypes. See sanitize_junctions().
verbose (bool)
checkpoint_every (int)

Returns:

(model, report, stats) — stats records chain, group, chain_type, n_clonotypes, n_used, p_nd2, loglik_first, loglik_last, iters, seconds.

Example

>>> model, rep, stats = build_model("TRB", work_dir="/tmp/em")
>>> model.save("models/TRB")

vdjtools.model.data.build_all(chains=('IGH', 'IGK', 'IGL', 'TRA', 'TRB', 'TRG', 'TRD'), *, groups=('human',), workers=None, out_dir=None, **kw)[source]#

Build models for several chains concurrently — the full-corpus entry point.

Each (group, chain) runs the whole build_model() pipeline. Concurrency is threads, not processes: arda annotation is a subprocess (so it blocks on I/O, not the GIL) and the native E-step releases the GIL, so a thread pool gets real parallelism without the memory cost of forking the germline and count arrays per chain.

Parameters:

chains – Chains to build. Defaults to all seven.
groups – Organism groups to build each chain for.
workers (int | None) – Concurrent builds. None = min(len(jobs), cpu_count // 2), leaving cores for each build’s own E-step threads.
out_dir – If given, each model is saved to out_dir/{group}_{chain}/.
**kw – Passed to build_model() (iters, tol, cap, gene_prior, …).

Returns:

{f"{group}_{chain}": {"model": ..., "report": ..., "stats": ...}} for the builds that succeeded, plus {"error": "..."} entries for those that did not — one failing chain never aborts the rest of a 30-minute corpus run.

Return type:

Example

>>> results = build_all(["TRB", "TRA"], workers=2, out_dir="models")
>>> pl.DataFrame([r["stats"] for r in results.values() if "stats" in r])

Generation probability, sampling and inference#

Generation probability (Pgen) of a CDR3 from a Model.

This is the reference nucleotide implementation: a transparent sum over every recombination scenario (V/J/D choice × deletions × insertion partition) that could produce the observed CDR3 nt, reading straight from the model’s polars tables (no OLGA at runtime). It reproduces OLGA’s compute_nt_CDR3_pgen to numerical tolerance (see the tests).

It is correctness-first, not speed-first — the fast amino-acid transfer-matrix DP is the job of the native _core port (Phase 1f). The nt scenario sum here is exactly the quantity the EM E-step needs (the bootstrap training data is out-of-frame nucleotide reads).

vdjtools.model.pgen.prepare(model)[source]#

Unpack a model’s polars tables into the lookup structures the Pgen sum uses.

Parameters:: model (Model)
Return type:: _Prepared

vdjtools.model.pgen.pgen_aa(prep, cdr3_aa, v=None, j=None)[source]#

Generation probability of an amino-acid CDR3, marginalizing over synonymous codons.

Parameters:

prep (_Prepared) – A prepare() -d model.
cdr3_aa (str) – The CDR3 amino-acid sequence (conserved Cys → conserved Phe/Trp inclusive).
v (str | None) – Optional gene names to restrict the sum to.
j (str | None) – Optional gene names to restrict the sum to.

Returns:

Pgen as a float. Sums tandem-D (n_D=2) scenarios for a D-D model, mixed with the single-D term by P(n_D) — consistent with pgen_nt().

Return type:

vdjtools.model.pgen.pgen_nt(prep, cdr3_nt, v=None, j=None)[source]#

Generation probability of a nucleotide CDR3, optionally restricted to a V and/or J.

Parameters:

prep (_Prepared) – A prepare() -d model.
cdr3_nt (str) – The CDR3 nucleotide sequence (conserved-Cys → conserved-Phe/Trp inclusive).
v (str | None) – Optional gene names to restrict the sum to (an OLGA-style usage mask of size 1).
j (str | None) – Optional gene names to restrict the sum to (an OLGA-style usage mask of size 1).

Returns:

Pgen as a float.

Return type:

Bridge from the polars Model to the native _core hot loops.

pack reconstructs the model’s dense arrays (gene choice, deletion, insertion, dinucleotide) and germline cut segments into a C++ PackedModel; pgen_nt calls the native Pgen. The result matches the pure-Python reference (and OLGA) exactly — the native path is just faster.

vdjtools.model.native.pack(model)[source]#

Build (and cache) the native PackedModel for model.

The native nt Pgen, aa Pgen (incl. Hamming-1 and v/j-agnostic), and the EM E-step (infer_native()) all support tandem-D (n_D=2).

The cache is keyed by id(model) but stores the model reference and verifies identity on hit: CPython recycles object ids after GC, so a bare-id cache could return a stale PackedModel for a different model that reused a freed id (e.g. running TRB then TRD EM in one process — the stale TRB pack has a different gene count and crashes the M-step). Keeping the ref also pins the id.

Parameters:: model (Model)

vdjtools.model.native.pgen_nt(model, cdr3_nt, v=None, j=None)[source]#

Native nucleotide Pgen — same result as vdjtools.model.pgen.pgen_nt(), faster.

Parameters:

model (Model)
cdr3_nt (str)
v (str | None)
j (str | None)

Return type:

vdjtools.model.native.pgen_aa(model, cdr3_aa, v=None, j=None, mismatches=0)[source]#

Native amino-acid Pgen — same result as vdjtools.model.pgen.pgen_aa(), much faster.

Parameters:

model (Model) – A recombination Model.
cdr3_aa (str) – The junction/CDR3 amino-acid sequence (Cys → Phe/Trp inclusive).
v (str | None) – V allele to condition on (e.g. "TRBV9*01"), or None to marginalize over all V (V-agnostic). A gene-level name ("TRBV9") or any unknown allele raises KeyError — it must not silently degrade to the V-agnostic value.
j (str | None) – J allele to condition on (e.g. "TRBJ2-3*01"), or None to marginalize (as v).
mismatches (int) – 0 for the exact sequence; 1 to also sum the Pgen of every amino-acid sequence within Hamming distance 1 (one substitution) — the total probability mass in the 1-mismatch ball, computed natively far faster than OLGA’s per-neighbour approach.

Returns:

Pgen as a float.

Return type:

vdjtools.model.native.pgen_aa_degenerate(model, allowed, v=None, j=None)[source]#

Total Pgen of every CDR3 matching a motif — a per-position set of permitted residues.

The same masked transfer-matrix DP pgen_aa() runs (that is the all-singletons case) and pgen_aa() with mismatches=1 runs per position (the one-wildcard case), with the residue set at each position free. Cost is one DP pass whatever the sets contain: the sequences the motif matches are summed over, never enumerated. This is the exact Pgen of a V/J/length- pinned motif such as a VDJdb cluster PWM thresholded per position.

Parameters:

model (Model) – A recombination Model.
allowed (list[str]) – One item per amino-acid position of the junction/CDR3 (so len(allowed) is the junction length). Each item is a string of the residues permitted there: "C" pins one, "ILVF" allows a subset, and "" or a string containing "X" is a wildcard (any of the 20 amino acids). list(seq) reproduces pgen_aa() exactly. A character the genetic code does not name raises ValueError — an empty codon mask would score a silent 0.0.
v (str | None) – V allele to condition on, or None to marginalize — as pgen_aa(), including the KeyError on a gene-level or unknown name.
j (str | None) – J allele to condition on, or None to marginalize (as v).

Returns:

Pgen as a float — the summed Pgen of all matching sequences.

Return type:

Note

X means wildcard here only. pgen_aa() matches residues by exact character against the genetic code, so pgen_aa(m, "CASSLAPGATNEKLXF") is 0.0, not a wildcard query. Pass list("CASSLAPGATNEKLXF") to this function to get the latter.

vdjtools.model.native.pgen_aa_degenerate_batch(model, allowed, v=None, j=None, threads=0)[source]#

Batch pgen_aa_degenerate() over many motifs, parallelized across queries.

Mirrors pgen_aa_batch(): the GIL is released and the queries are partitioned across worker threads, and the result is bitwise-identical to the serial per-motif computation for any threads.

Parameters:

model (Model) – A recombination Model.
allowed (list[list[str]]) – One per-position residue-set list per motif (see pgen_aa_degenerate()).
v (list[str | None] | None) – Optional per-motif V alleles (same length as allowed); None marginalizes over all V for every motif. Individual entries may be None.
j (list[str | None] | None) – Optional per-motif J alleles (as v).
threads (int) – Worker threads; 0 = auto (hardware_concurrency - 2). Batches under 64 motifs run single-threaded.

Returns:

Per-motif Pgen in input order.

Return type:

list[float]

vdjtools.model.native.pgen_aa_batch(model, cdr3_aas, v=None, j=None, mismatches=0, threads=0)[source]#

Batch amino-acid Pgen over many CDR3s, parallelized across sequences in native code.

Computes the same value as calling pgen_aa() per sequence, but releases the GIL and partitions the sequences across worker threads — the clean, exact speedup for the real workload (Pgen / 1-mismatch matching over many clonotypes). The result is bitwise-identical to the serial per-sequence computation for any threads.

Parameters:

model (Model) – A recombination Model.
cdr3_aas (list[str]) – Junction/CDR3 amino-acid sequences.
v (list[str | None] | None) – Optional per-sequence V alleles to condition on (same length as cdr3_aas); None marginalises over all V for every sequence. Individual entries may be None. An unknown or gene-level name raises KeyError (see pgen_aa()).
j (list[str | None] | None) – Optional per-sequence J alleles (as v).
mismatches (int) – 0 for exact Pgen, 1 for the Hamming-1 ball (as pgen_aa()).
threads (int) – Worker threads; 0 = auto (hardware_concurrency - 2). Batches under 64 sequences run single-threaded.

Returns:

Per-sequence Pgen in input order.

Return type:

list[float]

vdjtools.model.native.best_aa_scenarios(model, cdr3_aa, v=None, j=None, k=8)[source]#

Top-k recombination scenarios for an amino-acid CDR3, best first.

The argmax counterpart of pgen_aa(), over the same Pi_L*Pi_R transfer matrix: max in place of the sums, with the winning (V, delV) and (J, delJ) carried through the state. Because the DP marginalizes over V and J at no extra cost, leaving both unspecified is barely slower than pinning them (0.26 vs 0.23 ms on human TRB).

Parameters:

model (Model) – The recombination model.
cdr3_aa (str) – CDR3 amino-acid sequence (conserved Cys -> conserved Phe/Trp inclusive).
v (str | None) – Optional allele names to condition on. None marginalizes.
j (str | None) – Optional allele names to condition on. None marginalizes.
k (int) – How many scenarios to return.

Returns:

[(w, v_allele, len_v, j_allele, len_j, d_allele | None, idx5, idx3, pos)], descending by w. len_v/len_j are the nucleotides each germline contributes to the CDR3; idx5/idx3 are the D’s 5’/3’ trims and pos where its contribution starts.

Return type:

list[tuple]

Sample recombined sequences from a Model (OLGA-style generation).

Ancestral sampling over the model’s Bayes net: pick genes, deletions, insertion lengths and non-templated nt (via the dinucleotide Markov chain), assemble the CDR3, translate. Returns a polars DataFrame. This is the reference sampler; the native port is Phase 1f.

vdjtools.model.generate.prepare_generation(model)[source]#

Precompute cumulative distributions for fast ancestral sampling.

Parameters:: model (Model)
Return type:: _GenPrep

vdjtools.model.generate.generate(model, n, *, seed=None, productive_only=False)[source]#

Sample n recombined CDR3s from the model.

Parameters:

model (Model) – The generative model.
n (int) – Number of sequences to return.
seed (int | None) – RNG seed for reproducibility.
productive_only (bool) – If True, reject out-of-frame / stop-codon draws and keep sampling.

Returns:

DataFrame with junction_nt, junction_aa, v_call, d_call, d2_call, j_call, productive. d2_call is the second D of a tandem (n_D=2) draw, else null; for a single-D model it is all-null.

Return type:

DataFrame

EM inference of model marginals from nucleotide CDR3 sequences.

Expectation–Maximization over the recombination Bayes net: the E-step enumerates every scenario that could produce each observed nt CDR3 (the same enumeration as pgen), weights them by the current model, and accumulates soft counts per event realization; the M-step re-normalizes those counts in polars to get the next marginals. Trained on out-of-frame reads, it recovers the raw generation model (no productivity conditioning, so no selection bias).

Closed-loop oracle: generate synthetic sequences from a known model, then infer must recover that model’s marginals (see the tests). This is the reference driver; the E-step hot loop is a Phase 1f native-port candidate.

class vdjtools.model.infer.InferenceReport(loglik=<factory>, n_scoreable=<factory>, gene_tv=<factory>, n_iter=0, converged=False, n_sequences=0, max_iter=0, tol=0.0, init='', native=False, elapsed_s=0.0, finished_at='', template_source='')[source]#

Bases: object

Per-iteration diagnostics from infer() — the training log.

Every field after converged is metadata recorded so a saved model can say what it was fitted on and how; all are defaulted, so constructing a bare report still works.

Parameters:

loglik (list[float])
n_scoreable (list[int])
gene_tv (list[float])
n_iter (int)
converged (bool)
n_sequences (int)
max_iter (int)
tol (float)
init (str)
native (bool)
elapsed_s (float)
finished_at (str)
template_source (str)

loglik: list[float]#

n_scoreable: list[int]#

gene_tv: list[float]#

n_iter: int#

converged: bool#

n_sequences: int#

max_iter: int#

tol: float#

init: str#

native: bool#

elapsed_s: float#

finished_at: str#

template_source: str#

to_dict()[source]#

The report as a JSON-serializable dict (one entry of a model’s training["runs"]).

Return type:: dict

classmethod from_dict(obj)[source]#

Rebuild a report from to_dict(), ignoring keys this version does not know.

Parameters:: obj (dict)
Return type:: InferenceReport

to_frame()[source]#

Per-iteration training log: iter, loglik, n_scoreable, rel_change.

Return type:: DataFrame

vdjtools.model.infer.resume(path, sequences, **kw)[source]#

Continue EM from a checkpoint written by infer_native(checkpoint=...).

A warm start from the saved marginals, so the fit picks up where it stopped rather than realigning from scratch. The checkpoint carries its own training log and this run appends to it, so the full history survives across as many interruptions as it takes — which is what makes a multi-hour fit on a time-limited queue practical.

Parameters:

path – The checkpoint directory (or an already-loaded Model).
sequences – The same sequences the interrupted run was fitting.
**kw – Passed to infer_native(). init is forced to "template", and checkpoint defaults to path so the continued run keeps saving to the same place — pass checkpoint=None to stop checkpointing.

Returns:

(model, report) — report covers this run only; training_frame(model) shows every run.

Example

>>> model, rep = resume("ckpt/IGH", seqs, max_iter=10)

vdjtools.model.infer.print_progress(stream=None, prefix='')[source]#

A ready-made progress= callback that reports each EM iteration as it happens.

EM on a large D-bearing locus runs for tens of minutes with nothing to show for it, and the training log only becomes readable once the fit returns — so a long run is indistinguishable from a hung one. This prints the log-likelihood and its relative change per iteration, which is exactly the quantity the convergence test uses, so you can watch it approach tol.

Parameters:

stream – Where to write; defaults to sys.stderr so it never pollutes piped output.
prefix (str) – Prepended to each line, e.g. the locus being built.

Returns:

A callable suitable for infer(progress=...) / infer_native(progress=...).

Example

>>> infer_native(template, seqs, progress=print_progress(prefix="[TRB] "))
[TRB] iter  1  loglik -37.0067  rel      inf  n=122703
[TRB] iter  2  loglik -33.9738  rel 8.20e-02  n=122703

vdjtools.model.infer.training_frame(obj)[source]#

The training log of a model (or a single report) as one tidy frame across all runs.

Parameters:: obj – A Model carrying a training log, or an InferenceReport.
Returns:: run, iter, loglik, n_scoreable, rel_change — empty when the model was never fitted here (every bundled model, and anything imported straight from OLGA).
Return type:: DataFrame

Example

>>> m, rep = infer_native(template, seqs, max_iter=10)
>>> training_frame(m)     # loglik per iteration, ready to plot

vdjtools.model.infer.infer(template, sequences, *, max_iter=30, tol=0.001, init='align', masks=None, single_d=False, p_nd2_init=0.02, dd_allowed=None, nd_prior=0.0, progress=None, checkpoint=None, checkpoint_every=1)[source]#

Re-estimate a model’s marginals from nucleotide CDR3s by EM.

Parameters:

template (Model) – A model supplying the gene set, germline, and event graph (its marginals are replaced). Use one built by from_olga (or any Model).
sequences (list[str]) – Observed CDR3 nucleotide strings (typically out-of-frame reads).
max_iter (int) – Maximum EM iterations.
tol (float) – Stop when the V-usage total-variation between iterations falls below this.
init (str) – "align" (seed gene usage from a best-match vote — the default and fastest), "uniform" (each event uniform on its support), or "template" (warm start).
masks (list | None) – Optional per-sequence (v_genes, j_genes, d_genes) name lists (e.g. from arda_masks()) restricting each read’s scenario enumeration to its aligned genes. Strongly recommended for VDJ — without it the E-step enumerates every Cys-sharing V × the full D grid per read (tens of s/seq); with it, VDJ inference is tractable.
single_d (bool) – By default a tandem-D (D-D) model is learned for the D-bearing loci (IGH/TRD/TRB): a single-D template is promoted with to_dd() (seeding P(n_D=2)=p_nd2_init) and EM learns the true P(n_D=2). Set True to keep a strict single-D model. No effect on VJ loci or an already-tandem template.
p_nd2_init (float) – Initial P(n_D=2) seed when promoting to D-D (ignored if single_d).
dd_allowed (list | None) – Optional per-read booleans gating the tandem (n_D=2) E-step — a read may be tandem only where dd_allowed[i] is true (e.g. reads arda flags with a d2_call). Anchors D-D learning to alignment-detected tandems, countering the tandem-vs-long-insertion identifiability that inflates unregularized P(n_D=2) on real data. None = all reads.
nd_prior (float) – Dirichlet/Beta pseudocount added to the single-D (n_D=1) soft count each M-step, regularizing P(n_D=2) toward 0. Both anchors combine.
progress – Optional callable(iteration, loglik, rel_change, n_scoreable) invoked after every iteration — use print_progress() to watch a long fit converge live.
checkpoint – Directory to save the model into after each iteration, so a long fit survives being interrupted. Continue it with resume(). The checkpoint carries the training log so far, and the resumed run appends to it.
checkpoint_every (int) – Write a checkpoint every N iterations (default every one).

Returns:

(fitted_model, report). For a tandem-D template the E-step enumerates n_D=2 scenarios and the M-step learns P(n_D=2) along with the d2_gene / d2_del / dd events.

Return type:

tuple[Model, InferenceReport]

vdjtools.model.infer.call_alleles(index, call)[source]#

All model alleles compatible with one AIRR gene call, ambiguity included.

Two kinds of ambiguity, and both must widen the mask rather than narrow it:

allele — a call of TRBV20-1*03 where the truth is *01. Expanding to every model allele of the gene keeps the right scenario reachable.
comma-separated genes — AIRR writes an aligner’s tie as IGHV3-23*01,IGHV3-23D*01, which means the aligner could not tell these apart. Splitting on * alone keeps only IGHV3-23 and silently DROPS IGHV3-23D — a different gene on a duplicated locus. If the truth is the dropped one, its scenario is unreachable and EM misattributes the read. Measured on human IGH: 23,176 of 160,324 non-functional clonotypes (14.5%) carry an ambiguous V call; TRB 2.0%; TRA/TRD 0%.

Returns the union over every gene named, deduplicated and order-stable. Unknown genes contribute nothing; a call naming no known gene yields [], which the E-step reads as “unrestricted” — the honest degradation, since we know nothing about that read’s gene.

Parameters:

index (dict[str, list[str]])
call (str | None)

Return type:

list[str]

vdjtools.model.infer.gene_masks(model, v_calls, j_calls)[source]#

Build per-read (v_genes, j_genes, d_genes) E-step masks from V/J gene calls.

Each call is expanded to every model allele of every gene it names — see call_alleles() for why both allele- and comma-ambiguity must widen the mask. D is left unrestricted (few D genes, and D calls on the short D germline are unreliable).

Parameters:

model (Model)
v_calls (list[str])
j_calls (list[str])

Return type:

list[tuple]

vdjtools.model.infer.sanitize_junctions(df, col, *, ambiguous='A', where='infer_frame')[source]#

Make a junction column safe for the native encoder, which knows only A/C/G/T.

Real annotated reads carry the occasional ambiguous base (N, or an IUPAC code), and an unhandled one surfaces as a KeyError from deep inside the E-step.

Parameters:

df (DataFrame) – Clonotype frame.
col (str) – Junction column.
ambiguous (str | None) – A single base to substitute for every non-ACGT character (default "A"), or None to drop those clonotypes instead. Substituting is the default because it keeps the clonotype: an ambiguous base is one uncertain position in a junction that is otherwise perfectly good evidence, and on these reads it affects ~0.01% of rows, so dropping costs sample size for no gain in correctness. It is a substitution, not a marginalization — the base is treated as read, so a run with many ambiguous positions will bias the insertion model toward the substituted base and should use None.
where (str) – Caller name, used in the warning.

Returns:

A frame with the column uppercased and cleaned; rows are dropped only when ambiguous is None.

Return type:

DataFrame

vdjtools.model.infer.infer_frame(template, clones, *, seq_col=None, v_col='v_call', j_col='j_call', use_calls=True, native=True, ambiguous='A', **kw)[source]#

Fit a model from a clonotype frame — the ergonomic entry point to EM.

Wraps infer_native() with the two steps every caller otherwise repeats: find the nucleotide junction column, and turn the frame’s V/J calls into per-read E-step masks with gene_masks(). The masks matter enormously on a D-bearing locus — without them the E-step enumerates every Cys-sharing V against the full D grid for every read.

Parameters:

template – A Model supplying the gene set, germline and event graph, or a locus string (e.g. "TRB") to build one with from_arda().
clones (DataFrame) – Clonotype frame. Needs a nucleotide junction column and, for use_calls, v_call/j_call.
seq_col (str | None) – Explicit junction column; auto-detected from _JUNCTION_COLS otherwise.
v_col (str) – V-call column.
j_col (str) – J-call column.
use_calls (bool) – Build per-read masks from the V/J calls. Turn off only if the frame’s calls are untrustworthy — inference then enumerates every gene and gets much slower.
native (bool) – Use infer_native() (default). False runs the pure-Python infer().
ambiguous (str | None) – What to do with a junction holding a non-ACGT base — substitute this base (default "A"), or None to drop the clonotype. See sanitize_junctions().
**kw – Passed through (max_iter, tol, init, single_d, nd_prior, gene_prior, …).

Returns:

(fitted_model, report) — the model carries the run in its training log.

Raises:

ValueError – If no junction column is found, or it holds no usable sequences.

Example

>>> m, rep = infer_frame("TRB", clones, max_iter=10)
>>> training_frame(m)

vdjtools.model.infer.arda_masks(contigs, model, *, organism='human')[source]#

Annotate nt contigs with arda and build (junctions, masks) for masked infer().

The production path for real reads: junctions, masks = arda_masks(contigs, template); infer_native(template, junctions, masks=masks). arda is a base dependency (ships with vdjtools).

Parameters:

contigs (list[str])
model (Model)
organism (str)

Return type:

tuple[list[str], list[tuple]]

vdjtools.model.infer.augment_from_oracle(learned, oracle)[source]#

Fill functional genes the learned model left at P=0 with the ORACLE’s own usage and conditionals.

A learned model carries only the genes arda saw producibly in the training repertoire; a user’s library (different protocol/tissue) can be full of genes 5’RACE or this cohort never amplified, or that arda called but the germline can’t emit (a source mismatch, or a hard-call tie it lost). The OLGA oracle models every functional gene it knows with a real per-gene usage and deletion profile, and differs from the learned model ONLY in D/D-D handling — orthogonal to V/J gene identity — so its V/J genes transplant cleanly. For each functional (scoreable) gene present in the oracle but absent from the learned model, this copies the oracle’s choice mass AND every child table the gene parents (its deletion profile, P(J|V) on a VJ locus, P(D|J) for a J), then renormalizes — so no functional gene is silently missing. Usage is the oracle’s here; rescale_usage() adapts it to the user’s actual library (its cross-protocol job). Idempotent.

Parameters:

learned (Model)
oracle (Model)

Return type:

vdjtools.model.infer.enforce_dj_order(model)[source]#

Zero the genomically impossible (D, J) pairs in an existing model and renormalize.

The repair for a model fitted before this constraint existed — see forbidden_dj_pairs() for why TRBD2 cannot reach the TRBJ1 cluster. New fits get it from the M-step and need no repair.

Parameters:: model (Model) – A VDJ Model.
Returns:: A new model with P(D|J) (and P(D2|D)’s parent table, where applicable) masked and renormalized within each J. Unchanged for a VJ model or a locus with no interleaving.
Return type:: Model

Example

>>> fixed = enforce_dj_order(load_bundled("TRB", "learned", collapse=False))

vdjtools.model.infer.extend_alleles(model, germline, *, weight=1.0)[source]#

Add alleles from a larger germline library to an existing model, seeded from what it knows.

The use case is a model fitted against one reference meeting a richer one — a newer IMGT release, a population-specific library, your own genotyped alleles. Every new allele needs a germline row, a choice probability and a full set of child conditionals (its deletion profile, and P(J|V) / P(D|J) where it is a parent), none of which the library supplies.

Seeding uses the strongest evidence available for each case:

A new allele of a gene the model already has. Its choice mass is weight × the mean mass of that gene’s existing alleles, and its child tables are copied from a gene-mate. A new IMGT allele of a known gene is a polymorphism whose carriers use it about as often as the *01, so the gene’s own level is the right prior.
A brand-new gene. Child tables come from the germline-nearest existing allele, and the choice mass is a floor of half the smallest non-zero mass in the model. Plausible shape, deliberately tiny mass — there is no evidence at all for how often it is used.

Deletion rows copied from a donor are clipped to the new allele’s own germline length, so an extension can never introduce the unreachable mass check_model() flags.

Existing alleles are never modified, including their germline: silently swapping the sequence under an allele the model was fitted on would invalidate every conditional that references it. A library that disagrees about an existing allele is reported by check_model’s germline_source check, not fixed here.

Parameters:

model (Model) – The model to extend.
germline (DataFrame) – A germline frame (see from_germline() for the schema), typically a superset of the model’s own.
weight (float) – Scales the seeded mass for new alleles of known genes. 1.0 gives a new allele the gene’s average; 0.5 is a more conservative half of it.

Returns:

A new, validated and renormalized Model. Idempotent — extending twice with the same library changes nothing the second time.

Return type:

Note

This seeds, it does not estimate. Follow it with infer_native(extended, seqs, init="template") to let data set the new probabilities.

Example

>>> bigger = extend_alleles(m, load_germline("TRB", "human"))
>>> bigger, rep = infer_native(bigger, seqs, init="template", max_iter=5)

vdjtools.model.infer.infer_native(template, sequences, *, max_iter=30, tol=0.001, init='align', masks=None, single_d=False, p_nd2_init=0.02, dd_allowed=None, nd_prior=0.0, gene_prior=0.0, progress=None, checkpoint=None, checkpoint_every=1)[source]#

EM inference with the native C++ E-step — same result as infer(), much faster.

Requires the compiled _core extension. See infer() for the arguments (including single_d / p_nd2_init / dd_allowed / nd_prior / gene_prior).

gene_prior is a Dirichlet pseudocount spread over the germline’s functional V/J alleles in each M-step. P(V)=0 is an absorbing state of this EM — the E-step weights scenarios by P(V), so a zeroed allele can never be re-attributed — and on real data that silently deletes real genes for good (human TRB: 30 of 57 V genes survived unregularized). The germline says which alleles exist; the prior keeps all of them reachable and lets the data set the usage. 0.0 (default) is byte-identical to plain MLE. Learns tandem-D (n_D=2) by default on the D-bearing loci: the native E-step accumulates the second-D soft counts via a factorized forward/backward pass, read-parallelized across cores.

Parameters:

template (Model)
sequences (list[str])
max_iter (int)
tol (float)
init (str)
masks (list | None)
single_d (bool)
p_nd2_init (float)
dd_allowed (list | None)
nd_prior (float)
gene_prior (float)
checkpoint_every (int)

Return type:

tuple[Model, InferenceReport]

Tandem-D (D-D) extension#

Tandem-D (D-D) model extension: upgrade a single-D VDJ model to n_D ∈ {1, 2}.

OLGA (and every bootstrap model) assumes exactly one D per rearrangement. Real IGH and TRD repertoires contain tandem D-D joins — two D segments in the junction (V-D1-D2-J) — which no OLGA/IGoR model can represent. to_dd() adds the four events that a tandem rearrangement needs on top of an existing single-D Model:

n_d — P(n_D) over {1, 2} (the D-count prior; 2 = a tandem join).
d2_gene — P(D2 | D1) (the genomic-order mask lives in this table’s zeros).
d2_del — P(delD2_5', delD2_3' | D2) joint trimming of the second D.
dd_ins / dd_dinucl — the N-region between the two Ds.

The upgrade is an initialisation: it seeds P(n_D=2) = p_nd2 and copies the single-D deletion/insertion profiles for the new events. EM on real tandem-D reads then learns the true values (the enumeration in pgen_nt() already sums the n_D=2 scenarios). With p_nd2 = 0 the model is generatively identical to its single-D input.

vdjtools.model.dd.has_tandem(model)[source]#

True if the model places positive probability on a tandem (n_D=2) rearrangement.

Guards the paths that do not yet sum/sample the n_D=2 scenarios (amino-acid Pgen, the native _core Pgen/EM, generation, EM inference) so that a tandem model raises rather than silently returning a single-D result. The predicate is positive ``P(n_D=2)`` alone — so a to_dd(..., p_nd2=0) model (generatively single-D, d2_gene present but unused) correctly returns False and flows through the single-D fast paths byte-identically.

Parameters:: model (Model)
Return type:: bool

vdjtools.model.dd.to_dd(model, *, p_nd2=0.05)[source]#

Return a copy of model upgraded to a tandem-D (n_D ∈ {1, 2}) model.

Parameters:

model (Model) – A single-D VDJ Model (raises for VJ or already-DD).
p_nd2 (float) – Initial P(n_D = 2) prior mass (0 reproduces the single-D model exactly).

Returns:

A validated D-D Model with the four tandem events added.

Return type:

Model diagnostics, checking and scoring#

Entropy, mutual information and model-vs-model comparison live in analyze; the consistency audit in check; and everything that needs sequences — likelihood, BIC, Pgen distributions and diversity estimates — in score. See Recombination model workshop for the workflow these fit into.

Information-theoretic diagnostics for a recombination Model.

Turns a model’s declared Bayes net (events) and its marginal tables into three views used for validation and for the appendix figures:

entropy_table() — per-event Shannon entropy of that part of the rearrangement (bits): the marginal entropy H(X) of the event’s realization and, where the event is conditioned, the expected conditional entropy H(X | parents).
mutual_information() — the information each declared edge carries, I(child; parent) = H(child) − H(child | parent) (bits), plus I(V; J) and the within-D I(delD5; delD3).
bayes_net_dot() / render_bayes_net() — a graphviz DAG (bnlearn style) with nodes annotated by H(X) and edges by I; rendered to PDF/PNG via the dot CLI.

Everything is read straight from the polars tables, so it works identically on a legacy OLGA bootstrap model and an EM-inferred native model — the two are directly comparable.

vdjtools.model.analyze.gene_marginal(model, seg)[source]#

Marginal P(seg) as {allele: prob}, forward-propagated over the Bayes net.

V is always a root. VDJ J is a root (P(J)); VJ J is P(J|V) marginalized over V. D is P(D|J) marginalized over J.

Parameters:

model (Model)
seg (str)

Return type:

dict[str, float]

vdjtools.model.analyze.entropy_table(model)[source]#

Per-event entropy (bits): marginal H(X) and conditional H(X | parents).

Returns a tidy frame (event, kind, given, n_states, H_bits, H_cond_bits) — one row per event of the model’s declared graph, in graph order.

Parameters:: model (Model)
Return type:: DataFrame

vdjtools.model.analyze.mutual_information(model)[source]#

Mutual information (bits) carried by informative pairs of the model.

One row per declared parent→child edge (I(child; parent) = H(child) − H(child|parent)), plus I(V; J) (0 by construction for a VDJ model — V, J are independent roots) and the within-D deletion coupling I(delD5; delD3 | D) — for the second D too on a tandem model.

Parameters:: model (Model)
Return type:: DataFrame

vdjtools.model.analyze.total_entropy(model)[source]#

Per-event contribution to the scenario entropy H(recombination event), in bits.

For a Bayes net the joint entropy is the sum of each node’s conditional entropy given its parents, so this frame’s contribution_bits sums to the entropy of one whole recombination scenario — the information content of the process, and the basis of the 2^H diversity figure in vdjtools.model.score.diversity().

The insertion regions need care, and this is the only place the accounting is non-obvious. entropy_table() reports a dinucleotide event’s per-step conditional entropy (a composition summary, independent of how long the N-region is). Its contribution to the scenario is that per-step entropy times the expected number of steps, so the dinucleotide row here is E[length] · H_step using its paired *_ins event’s mean length. The first inserted base is treated as drawn from the chain’s stationary distribution rather than the model’s separate first-base bias — the standard decomposition, and worth a fraction of a bit at most.

Parameters:: model (Model) – The model to measure.
Returns:: event, kind, contribution_bits — one row per event, in graph order.
Return type:: DataFrame

Example

>>> total_entropy(m)["contribution_bits"].sum()   # bits per rearrangement

vdjtools.model.analyze.bayes_net_dot(model, *, title=None)[source]#

Graphviz DOT for the model’s Bayes net: nodes labelled with H(X), edges with I.

Parameters:

model (Model)
title (str | None)

Return type:

vdjtools.model.analyze.render_dot(dot, path, *, fmt='pdf')[source]#

Render any DOT source to path via the graphviz dot CLI; returns the output path.

Parameters:

dot (str) – DOT source, e.g. from bayes_net_dot() or compare_net_dot().
path (str | Path) – Output path; the suffix is replaced with fmt.
fmt (str) – Any format dot supports (pdf, png, svg, …).

Raises:

RuntimeError – If the dot CLI is not on PATH.

Return type:

vdjtools.model.analyze.render_bayes_net(model, path, *, fmt='pdf')[source]#

Render bayes_net_dot() to path via the dot CLI; returns the output path.

Parameters:

model (Model)
path (str | Path)
fmt (str)

Return type:

vdjtools.model.analyze.compare_entropy(models)[source]#

Stack entropy_table() across models into a wide event × model H(X) matrix.

Parameters:: models (dict[str, Model])
Return type:: DataFrame

vdjtools.model.analyze.compare_models(a, b, *, labels=('a', 'b'), by='allele')[source]#

Per-event distance between two models — the parameter-level compare_networks.

The two models’ tables are aligned on the union of their realization keys with zero fill, so a gene one model knows and the other does not contributes to the distance instead of being dropped. For a conditioned event the distance is computed per parent group and averaged weighted by the parent’s marginal (the same weighting the conditional-entropy code uses), so a rarely-used V’s deletion profile cannot dominate the number; tv_max reports the worst single group, which is what finds the one broken gene an average hides.

Parameters:

a (Model) – First model.
b (Model) – Second model.
labels (tuple[str, str]) – Names for the two models (used in error messages and the DOT title).
by (str) – "allele" (default) or "gene". Use "gene" to compare models built on different germline vintages or sources — an OLGA-namespace model against an arda-namespace one only lines up at gene level.

Returns:

event, kind, given, status, n_groups, support_a, support_b, support_shared, support_only_a, support_only_b, tv, tv_max, jsd_bits. status is shared, only_a, only_b, or schema_differs (the event exists in both but is factorized differently, e.g. j_choice is P(J|V) on a VJ locus and a root P(J) on a VDJ one). Distances are null unless the status is shared.

Return type:

One row per event of the union of both graphs

Note

Jensen-Shannon is the primary metric: it is symmetric, bounded by 1 bit, and finite when the supports differ, which is exactly the case here. KL is deliberately not reported — it is infinite whenever one model assigns zero to something the other does not.

Example

>>> compare_models(load_bundled("TRB", "olga"), load_bundled("TRB", "learned"), by="gene")

vdjtools.model.analyze.compare_usage(a, b, seg='v', *, by='gene')[source]#

Side-by-side gene (or allele) usage of two models — the protocol-bias view.

Parameters:

a (Model) – First model.
b (Model) – Second model.
seg (str) – "v", "j", "d" or "d2".
by (str) – "gene" (default) or "allele". Gene level is the meaningful comparison: allele calls on short reads are mismapping-prone, so allele-resolution usage is noise.

Returns:

name, p_a, p_b, log2_ratio sorted by descending p_a, over the union of both models’ genes. log2_ratio is null where either side is zero.

Return type:

DataFrame

vdjtools.model.analyze.compare_net_dot(a, b, *, labels=('a', 'b'), title=None)[source]#

Graphviz DOT contrasting two models’ Bayes nets — the compare_networks picture.

One DAG over the union of both graphs. Edges present in both are solid black; an edge only in a is blue and dashed, only in b red and dotted. Node fill intensity scales with that event’s Jensen-Shannon divergence between the two models, and node labels carry ΔH = H_a − H_b so the structural and the quantitative differences are visible at once.

Parameters:

a (Model) – First model.
b (Model) – Second model.
labels (tuple[str, str]) – Names for the two models, shown in the title and legend.
title (str | None) – Graph title; defaults to a summary of the two models.

Returns:

DOT source — render it with render_dot().

Return type:

Consistency audit of a Model against its own germline.

validate_tables() answers a narrow question — does every marginal normalize? — and raises on the first offender. This module answers the wider one a model builder actually has: is this model internally coherent, and does it agree with the germline it claims to be built on? It returns a tidy issue frame instead of raising, so every problem in a model is visible at once and the result is sortable, filterable and writable like any other table.

The checks exist because each one has silently produced a wrong answer at some point:

a functional gene left at P=0 makes Pgen exactly 0 for every clonotype using it, with no error;
deletion mass past a germline’s length is unreachable, so its probability is quietly lost;
an allele in a marginal but not in the germline frame crashes the native packer, far from the cause;
a model whose germline drifted from arda’s scores a different sequence than the one you annotated.

Severity is the contract: error means the model is broken (it will crash, or score wrongly), warn means it is suspicious but usable, info is a note.

vdjtools.model.check.DELETION_ENDS = {'d2_del': ('d_5', 'd_3'), 'd_del': ('d_5', 'd_3'), 'j_5_del': ('j_5',), 'v_3_del': ('v_3',)}#: Deletion event -> the palindrome_max key(s) bounding its most-negative (palindromic) ndel.

vdjtools.model.check.max_reachable_trim(event_name, kind, cut_len, palindrome_max)[source]#

Largest trim the Pgen DP can actually reach for a segment of cut_len nt.

The single definition of deletion reachability, shared by the checker, the collapse and the allele extension so the three cannot drift apart. Taken from pgen._v_options / pgen._d_middle rather than assumed: a trim consumes the palindrome-extended cut_segment and ndel = len(cut) - contributed - max_palindrome, so

V/J: ndel <= cut_len - Σ palindrome_max - 1 — the -1 is the invariant that V and J each contribute at least one nt to the CDR3;
D: ndel5 + ndel3 <= cut_len - Σ palindrome_max — a D may legally be deleted away.

Returns None for an event that is not a deletion, or one whose ends are not declared.

Parameters:

event_name (str)
cut_len (int)
palindrome_max (dict)

Return type:

int | None

vdjtools.model.check.check_model(model, *, germline='auto', tol=1e-05, raise_on=None)[source]#

Audit a model’s marginals against its manifest, its germline frames, and a reference library.

Parameters:

model (Model) – The model to check.
germline (str | DataFrame) – External germline to reconcile against. "auto" (default) tries load_germline() for the model’s own locus/organism and skips silently if arda has none (a custom library, or a non-arda organism); "none" disables the reconciliation; a pl.DataFrame uses that library directly.
tol (float) – Absolute tolerance for the “sums to 1 (or 0)” normalization check.
raise_on (str | None) – If "error" (or "warn"), raise ValueError when any issue at that severity or worse is found, instead of returning it.

Returns:

Tidy issue frame severity, check, event, segment, allele, detail, value — empty when the model is clean. severity is one of error, warn, info.

Raises:

ValueError – Only when raise_on is set and a matching issue is present.

Return type:

DataFrame

Example

>>> issues = check_model(load_bundled("TRB", "learned"))
>>> issues.filter(pl.col("severity") == "error").height
0

Score sequences under a model: likelihood, BIC, Pgen distributions, entropy and diversity.

Everything here is data-conditioned — it needs sequences, or sequences drawn from the model — which is what separates it from vdjtools.model.analyze (model-only information theory).

Three questions this answers:

How well does this model explain these sequences? model_fit() — log-likelihood, AIC, BIC.
Do two models score the same repertoire the same way? compare_pgen() + pgen_summary().
How much diversity does this model actually generate? diversity() — the entropy of the generated sequence distribution and the ~10^x effective-diversity number that follows from it.

Two conventions run through the module and are worth stating once.

Likelihoods use nucleotide Pgen. Σ Pgen_nt over all nt CDR3s is 1, so log Pgen_nt is a proper log-likelihood and BIC is meaningful. Pgen_aa sums only the in-frame, stop-free nucleotide fiber of a translation, so Σ Pgen_aa < 1: an amino-acid log-likelihood is unnormalized, and the missing constant differs between models. Amino-acid scoring is supported (real clonotype tables are often aa-only) but is a relative score on one fixed sequence set, never an absolute one.

Pgen = 0 never becomes ``-inf``. A sequence the model cannot generate gets pgen = 0.0 and a null log-probability, and aggregates are taken over the scoreable subset — the same convention InferenceReport already uses. Every aggregate reports n beside n_scoreable so a flattering log-likelihood earned on 10% of the data is visible rather than hidden.

vdjtools.model.score.pgen_frame(model, sequences, *, v=None, j=None, kind='auto', use_calls=True, on_unknown='error', threads=0, seq_col=None, v_col=None, j_col=None)[source]#

Generation probability of each sequence under model, as a tidy frame.

Parameters:

model (Model) – The recombination model to score with.
sequences – A list of junction/CDR3 strings, or a clonotype pl.DataFrame (the sequence column and v_call/j_call are auto-detected).
v – Optional per-sequence V calls, when sequences is a list.
j – Optional per-sequence J calls, when sequences is a list.
kind (str) – "nt", "aa", or "auto" (default) to detect from the sequences. "auto" requires a homogeneous set.
use_calls (bool) – Condition each sequence’s Pgen on its own V/J call. False marginalizes over all V/J — note this changes the quantity from P(junction, V, J) to P(junction).
on_unknown (str) – What to do with a call that does not resolve to exactly one model allele — "error" (default) or "marginalize".
threads (int) – Worker threads for amino-acid scoring (0 = auto). Nucleotide scoring is serial.
seq_col (str | None) – Explicit sequence column, when sequences is a frame.
v_col (str | None) – Explicit V-call column, when sequences is a frame.
j_col (str | None) – Explicit J-call column, when sequences is a frame.

Returns:

sequence, v_call, j_call, kind, pgen, log_pgen, log10_pgen, scoreable — one row per input sequence, in input order. log_pgen is null where pgen is 0.

Raises:

KeyError – If a V/J call does not resolve and on_unknown="error".
ValueError – On a mixed nt/aa set under kind="auto", or a missing sequence column.

Return type:

DataFrame

vdjtools.model.score.free_params(model, *, by_event=False, eps=0.0, tol=1e-09, reachable_only=True)[source]#

Number of free parameters k in a model — the penalty term of AIC / BIC.

k = Σ_events Σ_groups max(support − 1, 0), where a group is one normalization group (normalization_keys(), which already splits a dinucleotide table by from_nt, giving the right 4 × (4 − 1) = 12), the −1 is the simplex constraint, and support is the number of cells carrying probability, not the number of rows.

Counting rows instead of support would be badly wrong: deletion bins past a germline’s length, alleles pinned to zero and insertion tail bins are structural zeros that parameterize nothing. On human TRB that is the difference between ~3,600 and ~700 parameters for v_3_del alone.

Two group kinds contribute nothing and are dropped:

undefined conditionals — the all-zero groups the schema explicitly permits, kept only for gene-index alignment;
unreachable conditionals — a group whose parent allele has zero marginal probability (from_olga fills P(D|J) uniformly even where P(J) = 0). These sum to 1 but are not estimable, and counting them inflates k by n_D − 1 for every dead J.

Caveat, and it is a real one: a support-based count cannot tell a structural zero from a parameter EM happened to drive to exactly zero, so k is a lower bound. BIC is therefore only comparable between models counted the same way — which is the case for any two models compared through this function.

Parameters:

model (Model) – The model to count.
by_event (bool) – Return a per-event breakdown instead of the total.
eps (float) – A cell counts toward the support when p > eps.
tol (float) – A group whose probabilities sum to at most this is treated as an undefined conditional.
reachable_only (bool) – Drop groups whose parent allele has zero marginal probability.

Returns:

k as an int, or (with by_event) a frame event, n_groups, k.

Return type:

int | DataFrame

vdjtools.model.score.model_fit(model, sequences, *, weights=None, k=None, **kw)[source]#

Log-likelihood, AIC and BIC of a sequence set under a model.

Parameters:

model (Model) – The model to score with.
sequences – Sequences or a clonotype frame — see pgen_frame().
weights – Optional per-sequence weights (e.g. duplicate_count), or the name of a column when sequences is a frame. The log-likelihood becomes Σ w·log p and n becomes Σ w, i.e. exactly as if each clonotype were repeated w times.
k (int | None) – Override the free-parameter count; defaults to free_params().
**kw – Passed to pgen_frame() (kind, use_calls, on_unknown, threads, seq_col, v_col, j_col, v, j).

Returns:

kind, conditioned, n, n_scoreable, frac_scoreable, loglik_sum, loglik_mean, k, aic, bic. conditioned records whether V/J were conditioned on, because that changes the sample space (P(junction, V, J) vs P(junction)) and two fits are only comparable when it matches.

Return type:

A one-row frame

Example

>>> model_fit(load_bundled("TRB", "learned"), held_out_junctions)

vdjtools.model.score.compare_pgen(a, b, sequences, *, labels=('a', 'b'), **kw)[source]#

Score one sequence set under two models and pair the results up.

Parameters:

a (Model) – First model.
b (Model) – Second model.
sequences – Sequences or a clonotype frame — see pgen_frame().
labels (tuple[str, str]) – Names for the two models, used in the output column suffixes.
**kw – Passed to pgen_frame() for both models, so the two are scored identically.

Returns:

sequence, v_call_<a>, v_call_<b>, j_call_<a>, j_call_<b>, kind, pgen_<a>, pgen_<b>, log10_<a>, log10_<b>, delta_log10 — one row per sequence. The two V/J columns differ when the models resolve a gene-level call to different alleles. delta_log10 is null unless both models scored the sequence.

Return type:

DataFrame

Example

>>> summary = pgen_summary(compare_pgen(olga, learned, seqs, labels=("olga", "learned")))

vdjtools.model.score.pgen_summary(cmp, *, labels=('a', 'b'))[source]#

Summarise a compare_pgen() frame: agreement, offset, and coverage.

Parameters:

cmp (DataFrame) – The output of compare_pgen().
labels (tuple[str, str]) – The same labels that call used.

Returns:

n, n_scoreable_a, n_scoreable_b, n_scoreable_both, only_a_scoreable, only_b_scoreable, mean_log10_a, mean_log10_b, median_log10_a, median_log10_b, mean_delta, median_delta, sd_delta, q05_delta, q95_delta, pearson_log10, spearman_log10, ks_stat, ks_p.

Return type:

A one-row frame

Note

only_a_scoreable / only_b_scoreable are the headline numbers, not the correlations: one model assigning Pgen 0 to thousands of sequences the other scores fine is the finding, and a mean-delta-only report would hide it entirely. The KS statistic compares the two marginal log10 distributions; the samples are paired, so its p-value is anticonservative — lean on the statistic.

vdjtools.model.score.pgen_spectrum(model, *, n=10000, seed=0, bins=40, productive_only=False, sequences=None, **kw)[source]#

The model’s Pgen distribution, as a histogram table ready to plot.

Parameters:

model (Model) – The model whose Pgen spectrum to take.
n (int) – Sequences to generate when sequences is not given.
seed (int) – Generation seed.
bins (int) – Number of equal-width log10 bins.
productive_only (bool) – Restrict generation to productive rearrangements.
sequences – Score these instead of generating (e.g. a real repertoire), so an observed spectrum can be overlaid on the model’s own.
**kw – Passed to pgen_frame().

Returns:

bin_left, bin_right, bin_mid, count, frac over log10(Pgen).

Return type:

DataFrame

vdjtools.model.score.diversity(model, *, n=5000, seed=0, productive_only=False)[source]#

How much diversity this model generates — entropy in bits, and effective diversity.

Two independent readings of “total diversity”, reported side by side because they answer different questions and differ by orders of magnitude:

Scenario entropy H_scenario — the information in one recombination event, summed over the Bayes net (total_entropy()). This is an upper bound on the sequence entropy, because different scenarios can produce the same junction.
Sequence entropy H_sequence — the entropy of the junction distribution itself, estimated by Monte Carlo. Sequences drawn from the model are distributed as Pgen, so E[−log₂ Pgen] over generated sequences is an unbiased estimator of H and the sample standard error comes free.

From those, two Hill numbers:

diversity_shannon = 2^H_sequence (Hill q = 1) — the classic “effective number of distinct sequences” figure, the one usually quoted as ~10^x for a locus.
diversity_simpson = 1 / E[Pgen] (Hill q = 2) — the inverse coincidence probability, i.e. how many sequences you would need for two independent draws to collide. Exact, because Σ Pgen² = E_{s∼Pgen}[Pgen], so the generated sample estimates it directly. It is always the smaller number: it weights the common sequences more heavily.

Parameters:

model (Model) – The model to characterise.
n (int) – Sequences to generate. The Shannon estimate converges quickly; diversity_simpson is driven by the few highest-Pgen draws and needs more samples for a tight interval — check pgen_mean_se against pgen_mean.
seed (int) – Generation seed, so the estimate is reproducible.
productive_only (bool) – Estimate over productive rearrangements only. Off by default: the unrestricted distribution is the one Pgen normalizes over, so the estimator is unbiased there and merely descriptive here.

Returns:

n, scenario_entropy_bits, scenario_diversity, sequence_entropy_bits, sequence_entropy_se_bits, diversity_shannon, diversity_simpson, pgen_mean, pgen_mean_se, median_log10_pgen, q05_log10_pgen, q95_log10_pgen.

Return type:

A one-row frame

Example

>>> diversity(load_bundled("TRB", "olga"), n=20_000)

Input/output and schema (`vdjtools.io`)#

The canonical clonotype frame (AIRR Rearrangement column names + polars dtypes) and every reader that emits it: native vdjtools tables, AIRR Rearrangement TSV, and Parquet, plus format auto-detection, metadata-driven batch loading, and hive-partitioned cohort scans.

`vdjtools.io.schema`#

Canonical clonotype schema and coercion helpers.

Canonical clonotype-frame schema (AIRR-aligned) and coercion helpers.

The basic-analytics layer speaks a single, flat clonotype frame — one row per clonotype, AIRR Rearrangement column names, polars dtypes. Every reader emits it and every analysis function consumes it. Kept deliberately minimal (free functions, no classes) to mirror the vdjmatch / arda convention.

Columns:

v_call, d_call, j_call, c_call (Utf8, nullable) — IMGT segment calls; c_call is frequently absent in native vdjtools data.
junction_aa (Utf8) — the junction amino-acid sequence (conserved anchors Cys104 … Phe/Trp118 INCLUDED), per the AIRR junction_aa convention (equivalently the legacy vdjtools cdr3aa). This is two residues longer than the IMGT cdr3_aa (anchors excluded); readers prefer the junction form.
junction_nt (Utf8, nullable) — the junction nucleotide sequence (anchors included), matching junction_aa above. AIRR spells the nucleotide junction junction (no _nt suffix); readers accept that (and legacy cdr3_nt) as input aliases.
duplicate_count (Int64) — read/UMI count for the clonotype.
frequency (Float64) — duplicate_count normalised within the sample.
locus (Utf8, derived) — first three characters of v_call (TRB, IGH …).

vdjtools.io.schema.PRODUCTIVE = 'productive'#

OPTIONAL AIRR Rearrangement annotation columns. Not in SCHEMA – they are not required and most bulk pipelines do not emit them – but when a file DOES carry them they are authoritative and vdjtools.preprocess.filter_productive() reads them in preference to re-deriving the same fact from junction_aa. See https://docs.airr-community.org/en/latest/datarep/rearrangements.html

productive is the composite: an open reading frame, no defect in the start codon, splicing sites or regulatory elements, no internal stop codon, and an in-frame junction. The other two are components of it, useful when the composite is absent.

vdjtools.io.schema.AIRR_FUNCTIONAL_COLUMNS: tuple[str, ...] = ('productive', 'stop_codon', 'vj_in_frame')#: The optional AIRR annotation columns, in the order they are preferred as evidence.

vdjtools.io.schema.SCHEMA: dict[str, DataType] = {'c_call': String, 'd_call': String, 'duplicate_count': Int64, 'frequency': Float64, 'j_call': String, 'junction_aa': String, 'junction_nt': String, 'v_call': String}#: Canonical columns in canonical order, mapped to their polars dtype.

vdjtools.io.schema.COLUMNS: list[str] = ['v_call', 'd_call', 'j_call', 'c_call', 'junction_aa', 'junction_nt', 'duplicate_count', 'frequency']#: Column names in canonical order.

vdjtools.io.schema.column_names(df)[source]#

Column names of an eager or lazy frame.

Uses polars.LazyFrame.collect_schema() for a LazyFrame so the check does not emit polars’ “resolving schema” performance warning; falls back to .columns for an eager DataFrame.

Parameters:: df (DataFrame | LazyFrame) – A pl.DataFrame or pl.LazyFrame.
Returns:: The list of column names.
Return type:: list[str]

vdjtools.io.schema.locus_of(v_call)[source]#

Return the locus (first three characters) of an IMGT V-gene call.

Parameters:: v_call (str | None) – An IMGT V-gene call such as "TRBV12-3*01", or None.
Returns:: The three-letter locus ("TRB"), or None if v_call is None or shorter than three characters.
Return type:: str | None

Example

>>> locus_of("TRBV12-3*01")
'TRB'

vdjtools.io.schema.add_locus(df)[source]#

Add (or overwrite) the derived locus column from v_call.

Parameters:: df (DataFrame) – A clonotype frame carrying a v_call column.
Returns:: The frame with a locus column (null where v_call is null).
Return type:: DataFrame

vdjtools.io.schema.recompute_frequency(df)[source]#

Recompute frequency as duplicate_count / sum(duplicate_count).

Parameters:: df (DataFrame) – A clonotype frame with a duplicate_count column.
Returns:: The frame with frequency overwritten. If the total count is zero the frequency is set to 0.0 for every row.
Return type:: DataFrame

vdjtools.io.schema.normalize(df, *, recompute_freq=False, keep=())[source]#

Coerce an arbitrary frame to the canonical clonotype schema.

Missing canonical columns are added as nulls, present ones are cast to their declared dtype (non-strict — unparseable values become null). The result is the canonical columns in canonical order, followed by any keep columns; every other non-canonical column (e.g. native vdjtools markup like VEnd/DStart) is dropped.

Parameters:

df (DataFrame) – A frame that already uses canonical column names for whatever columns it carries.
recompute_freq (bool) – If True, recompute frequency from duplicate_count after coercion (use when the source lacks a trustworthy frequency).
keep (tuple[str, ...]) – Non-canonical columns to preserve, e.g. ("v_identity",). Dtypes are left alone — the canonical schema has nothing to say about them.

Returns:

A frame with the canonical columns, correctly typed and ordered, plus keep.

Return type:

DataFrame

vdjtools.io.schema.weight_expr(weight)[source]#

Return the per-clonotype weight expression for an analysis mode.

Parameters:: weight (str) – One of "reads" (weight by duplicate_count), "unique" (one per clonotype), or "freq" / "frequency" (weight by frequency).
Returns:: A polars expression yielding the per-row weight.
Raises:: ValueError – If weight is not a recognised mode.
Return type:: Expr

vdjtools.io.schema.strip_allele(expr)[source]#

Reduce a segment-call expression to gene resolution, ambiguity-safe.

Strips the IMGT allele suffix from every gene an AIRR call names, not just the first. The old \*.*$ regex matched from the FIRST * to end of string, so a comma-ambiguous call like IGHV3-23*01,IGHV3-23D*01 collapsed to IGHV3-23 – silently dropping IGHV3-23D, which then reported zero usage across a whole cohort despite being named in tens of thousands of rows. Genes are de-duplicated after stripping, so an allele-level tie within one gene (IGHV1-2*02,IGHV1-2*04) correctly collapses to the single unambiguous gene IGHV1-2, while a genuine cross-gene tie stays IGHV3-23,IGHV3-23D.

Parameters:: expr (Expr) – A polars string expression over segment calls.
Returns:: Each call reduced to its distinct gene(s), sorted and comma-joined (TRBV12-3*01 → TRBV12-3; A*01,A*02 → A; A*01,B*01 → A,B); nulls pass through unchanged.
Return type:: Expr

vdjtools.io.schema.resolve_gene(expr)[source]#

Reduce a segment call to exactly ONE gene: allele stripped, ambiguity resolved to the first.

The companion to strip_allele(), and the distinction matters:

strip_allele() keeps a genuine cross-gene tie as IGHV3-23,IGHV3-23D, because when you are reporting usage you must not invent certainty the aligner did not have.
resolve_gene() collapses it to IGHV3-23, because when the gene is a feature axis every distinct ambiguity string otherwise becomes its own category.

That second failure is not hypothetical. Fitting a V+k-mer vocabulary on 200 HIP samples produced 1,296 V “genes”, 1,235 of them comma-strings such as TRBV1,TRBV23-1,TRBV4-1,TRBV4-2,TRBV4-3. The cost is not the 21x wider axis: it is that the real TRBV9 bucket gets drained, since every TRBV9 clone that happened to be called ambiguously was filed elsewhere. Its features then fall below any incidence floor and vanish, so a cohort with clean calls is scored against columns nobody populated.

Where the ambiguity comes from matters, because it is not a parsing artifact and will not go away: that cohort is Adaptive/immunoSEQ realigned with MiXCR against IMGT from the junction plus short flanks. The realignment is the better call – MiXCR/IMGT is a sounder reference than Adaptive’s own – and the ambiguity is what honest calling looks like when V genes differ only outside the sequenced window. So first-listed is a resolution, not a correction.

Two things it cannot fix, and which belong to the assay rather than the reference: the window still bounds what is resolvable, and Adaptive’s multiplex V primers distort V usage frequencies. A V-conditioned feature axis fitted on such a cohort inherits both, which makes it a poor donor for a 5’RACE cohort whatever this function does.

First-listed rather than dropped: an ambiguous call still carries a clonotype, and the aligner lists its best call first. Note this takes the first call as written, not strip_allele(...).list.first() – strip_allele sorts, so composing the two silently returns the alphabetically-first gene instead (TRBV5*01,TRBV19*03 -> TRBV19, not TRBV5).

Parameters:: expr (Expr) – A polars string expression over segment calls.
Returns:: One gene per row (TRBV12-3*01 -> TRBV12-3; A*01,B*01 -> A); nulls pass through unchanged.
Return type:: Expr

`vdjtools.io.read`#

Single-file readers (native vdjtools, AIRR, Parquet).

Readers for the native vdjtools table format and AIRR Rearrangement TSV.

Both readers return the canonical clonotype frame described in vdjtools.io.schema. Gzip-compressed inputs (.gz) are handled transparently by polars.

vdjtools.io.read.read_vdjtools(path, n_rows=None, *, keep=(), recompute_frequencies=True)[source]#

Read a native vdjtools clonotype table into the canonical frame.

The native header is count freq cdr3nt cdr3aa v d j VEnd DStart DEnd JStart followed by optional annotation columns; . and empty strings mean missing. Ambiguous V/D/J calls ("TRBV12-4, TRBV12-3") are reduced to the first (highest-confidence) call. c_call is always null (the format has no C gene). frequency is recomputed exactly from duplicate_count (the native freq column is a rounded copy of the same ratio).

Parameters:

path (str | PathLike) – Path to a .txt or .txt.gz native vdjtools table.
n_rows (int | None) – If given, read at most this many data rows (preview huge files).
keep (tuple[str, ...]) – Annotation columns to preserve past the canonical set (this format’s trailing annotation fields are otherwise dropped). Same contract as the other two readers: absent names are skipped, not raised.
recompute_frequencies (bool)

Returns:

Canonical clonotype frame with a derived locus column.

Raises:

ValueError – If the required count and cdr3aa columns are absent.

Return type:

DataFrame

vdjtools.io.read.read_parquet(path, n_rows=None, *, keep=(), recompute_frequencies=True)[source]#

Read a Parquet clonotype table into the canonical frame.

Parquet is the at-scale storage format for repertoire cohorts: typed, columnar, compressed, and — laid out as a hive-partitioned directory — scannable as one polars.LazyFrame with predicate/projection pushdown (see vdjtools.io.cohort). Columns already using the canonical names (e.g. a file written by polars.DataFrame.write_parquet() from a canonical frame) are kept as-is; AIRR source names (junction_aa, duplicate_count …) are mapped to canonical only to fill a canonical column that is otherwise absent. Unlike the TSV path this reads native dtypes directly (no all-Utf8 pass), so it is both faster and lower-peak-memory on large files.

Parameters:

path (str | PathLike) – Path to a .parquet / .pq clonotype table.
n_rows (int | None) – If given, read at most this many data rows (preview huge files).
keep (tuple[str, ...]) – Extra source columns to carry through unrenamed, e.g. ("v_identity",). The canonical set is deliberately small, but a feature can need a field outside it – the signature’s SHM block is a weighted mean of v_identity, which no canonical column holds, so without this it is not merely missing but uncomputable, and ships as a permanently-nan column. Names absent from the file are skipped rather than raising: keep says what to preserve if present, not what to require.
recompute_frequencies (bool)

Returns:

Canonical clonotype frame with a derived locus column.

Raises:

ValueError – If no CDR3 amino-acid column (cdr3_aa / junction_aa) is present.

Return type:

DataFrame

vdjtools.io.read.read_airr(path, *, collapse=True, n_rows=None, keep=(), recompute_frequencies=True)[source]#

Read an AIRR Rearrangement TSV into the canonical frame.

Prefers the junction columns junction_aa / junction (conserved anchors INCLUDED — the canonical vdjtools / AIRR convention) over the IMGT cdr3_aa / cdr3 columns (anchors excluded), which are used only as a fallback. Known limitation: a file that provides only the IMGT cdr3_aa yields sequences two residues shorter than the junction, so downstream lengths / k-mers will differ. The count column may be duplicate_count or a vdjtools-style count / reads (AIRR-hybrid exports); absent, it defaults to 1. Per-read files (one row per rearrangement) collapse to unique clonotypes with summed counts; already-aggregated files pass through unchanged. frequency is always recomputed after collapsing.

Clonotype identity for collapsing is (v_call, j_call, junction_nt, junction_aa) — matching legacy Clonotype equality (V, J, CDR3nt); d_call and c_call are not part of the identity, so a representative (first non-null) value is attached to each collapsed clonotype.

Parameters:

path (str | PathLike) – Path to a .tsv / .tsv.gz AIRR Rearrangement file.
collapse (bool) – If True (default), sum the count over clonotypes identical on (v_call, j_call, junction_nt, junction_aa).
n_rows (int | None) – If given, read at most this many data rows (preview huge files).
keep (tuple[str, ...]) – Extra source columns to carry through unrenamed, e.g. ("v_identity",) for the signature’s SHM block, which no canonical column holds. Under collapse a numeric kept column is averaged over the clonotype’s reads and anything else takes a representative value. Names absent from the file are skipped rather than raising.
recompute_frequencies (bool)

Returns:

Canonical clonotype frame with a derived locus column.

Raises:

ValueError – If no CDR3 amino-acid column (junction_aa or cdr3_aa) is present.

Return type:

DataFrame

`vdjtools.io.batch`#

Format sniffing and metadata-driven batch / streaming reads.

Format auto-detection and metadata-driven batch reading.

Sits on top of the single-file readers in vdjtools.io.read: sniffs a file’s format, dispatches to the right reader, and loads a whole cohort described by a metadata table into one long clonotype frame (or a per-sample dict).

vdjtools.io.batch.sniff_format(path)[source]#

Detect a clonotype file’s format from its header.

Parameters:: path (str | PathLike) – Path to a clonotype table.
Returns:: "parquet" (.parquet / .pq extension), "vidjil" (.vidjil / .json or a leading {), "imgt", "migec", "mitcr", "rtcr", "mixcr", "immunoseq", "trust4" (each by its signature header columns), "airr_cell" (AIRR + a cell_id barcode — single-cell, see read()), "arda" (AIRR + arda’s d2_call), "vdjtools" (native table / MigMap — cdr3aa / count + cdr3nt), or "airr" (AIRR Rearrangement — v_call / junction_aa / junction_nt / cdr3_aa).
Return type:: The detected format string, one of
Raises:: ValueError – If no known format signature is recognised.

vdjtools.io.batch.read(path, fmt='auto', n_rows=None, *, keep=(), recompute_frequencies=True)[source]#

Read a clonotype table, auto-detecting the format by default.

Parameters:

path (str | PathLike) – Path to a native vdjtools, AIRR Rearrangement, or Parquet table, or a third-party tool export (MiXcr, MiGec, MiTCR/tcR, immunoSEQ, IMGT/HighV-QUEST, Vidjil, RTCR, TRUST4, arda) — see vdjtools.io.convert (.gz ok for the text formats).
recompute_frequencies (bool) – Derive frequency from duplicate_count (default, and the historical behaviour). Pass False to use the frequencies as they are in the file – which matters when the source carries a UMI-corrected or already-normalised frequency that count/total would silently discard. A file with no frequency column derives one either way; there is nothing to preserve. Only the vdjtools, airr and parquet readers honour it – the legacy converters build frequency themselves.
fmt (str) – "auto" (sniff the header / extension), "vdjtools", "airr", "parquet", or a legacy format string ("mixcr", "migec", "mitcr", "immunoseq", "imgt", "vidjil", "rtcr", "trust4", "arda").
n_rows (int | None) – If given, read at most this many data rows (preview huge files).
keep (tuple[str, ...]) – Non-canonical columns to preserve if present, e.g. ("v_identity",) for the signature’s SHM block. Honoured by the three native readers; the legacy converters narrow to the canonical schema and ignore it.

Returns:

Canonical clonotype frame.

Raises:

ValueError – If fmt is unknown, auto-detection fails, or the file is a single-cell AIRR table (a cell_id column) — those go through vdjtools.sc.read_airr_cell(), since collapsing them here would drop the cell barcode. Pass fmt="airr" to pool one into a bulk repertoire anyway.

Return type:

DataFrame

vdjtools.io.batch.read_metadata(path)[source]#

Read a sample metadata TSV.

The literal string "nan" (and empty strings) are treated as null. All columns are read as strings so metadata joins are stable. A leading # on the first column name (some metadata sheets comment out the header line, e.g. #file_name sample_id ...) is stripped.

Parameters:: path (str | PathLike) – Path to a metadata TSV (one row per sample).
Returns:: A pl.DataFrame of metadata, all-Utf8, with "nan" → null.
Return type:: DataFrame

vdjtools.io.batch.iter_samples(metadata, base_dir, sample_col='sample_name', file_template='{sample}.tsv.gz', fmt='auto', add_metadata=True)[source]#

Yield (sample_id, frame) one sample at a time — O(1-sample) RAM.

The streaming counterpart to read_samples(). Because it reads and yields each sample in turn (never accumulating), a caller can reduce-and-discard — per-sample diversity/usage, or sink each sample to a Parquet partition — and so process a 100k-sample cohort whose concatenation would not fit in memory. Each yielded frame carries the canonical columns + locus + the reserved sample_id / file_name, and (if add_metadata) that row’s metadata.

Parameters:

metadata (DataFrame) – Metadata frame (e.g. from read_metadata()).
base_dir (str | PathLike) – Directory holding the per-sample clonotype files.
sample_col (str) – Metadata column carrying the sample name.
file_template (str) – Filename template; {sample} is substituted with the sample name.
fmt (str) – Format for the readers ("auto" / "vdjtools" / "airr" / "parquet").
add_metadata (bool) – If True, attach the metadata columns to every clonotype row.

Yields:

(sample_id, frame) tuples in metadata order.

vdjtools.io.batch.read_samples(metadata, base_dir, sample_col='sample_name', file_template='{sample}.tsv.gz', fmt='auto', add_metadata=True, as_dict=False)[source]#

Read a batch of samples described by a metadata table into one frame.

Eagerly materialises the whole cohort — convenient for tens-to-hundreds of samples, but it holds every sample in RAM at once. For large cohorts (thousands+ of samples) use iter_samples() to stream, or vdjtools.io.cohort() to persist a hive-partitioned Parquet dataset and scan it lazily.

For each metadata row the file base_dir/file_template.format(sample=<sample>) is read into the canonical schema, tagged with the reserved columns sample_id (the sample_col value) and file_name, and — if add_metadata — joined with that row’s metadata columns.

Parameters:

metadata (DataFrame) – Metadata frame (e.g. from read_metadata()).
base_dir (str | PathLike) – Directory holding the per-sample clonotype files.
sample_col (str) – Metadata column carrying the sample name (default "sample_name").
file_template (str) – Filename template; {sample} is substituted with the sample name (default "{sample}.tsv.gz").
fmt (str) – Format for the readers ("auto" / "vdjtools" / "airr" / "parquet").
add_metadata (bool) – If True, attach the metadata columns to every clonotype row.
as_dict (bool) – If True, return {sample_id: frame} instead of one long frame.

Returns:

One concatenated long pl.DataFrame (canonical columns + locus + sample_id + file_name + metadata), or a dict[str, pl.DataFrame] if as_dict.

vdjtools.io.batch.map_samples(fn, items, *, fmt='auto', workers=None, keep=())[source]#

Read each sample and apply fn to it, in parallel, low peak memory.

The parallel, streaming replacement for “load the whole cohort, then loop over it”: each worker reads one sample, reduces it via fn (typically a per-sample summary statistic returning a small frame), and discards the raw clonotype frame — so peak memory is O(workers) samples, not the whole cohort. Reading (gzip TSV decode + polars parse) and fn (polars group_by / numpy) both release the GIL, so a thread pool runs them genuinely in parallel with no frame pickling.

Parameters:

fn – A callable pl.DataFrame -> T applied to each sample’s canonical clonotype frame (e.g. vdjtools.stats.diversity.diversity_stats()).
items – Iterable of (sample_id, path) pairs (e.g. from a metadata sheet).
fmt (str) – Reader format passed to read() ("auto" sniffs each file).
workers (int | None) – Max worker threads. None uses the pool default (min(32, os.cpu_count() + 4)); pass a smaller value if fn is compute-bound, to avoid oversubscribing polars’ own thread pool.
keep (tuple[str, ...]) – Non-canonical columns to preserve, passed to read(). Without it a reduction needing a field outside the canonical eight — v_identity for the SHM block — gets a frame that never carried it, and reports a hole rather than a number.

Returns:

A list of (sample_id, fn(frame)) in the same order as items (input/metadata order, regardless of which sample finishes first).

`vdjtools.io.cohort`#

Hive-partitioned Parquet cohort writer / lazy scanner.

Cohort-scale I/O: a hive-partitioned Parquet dataset scanned as one LazyFrame.

The single-file readers and vdjtools.io.read_samples() materialise a whole cohort in RAM — fine for tens-to-hundreds of samples, fatal for 100k. This module provides the at-scale path:

ingest_cohort() streams every sample through the readers one at a time (never holding two at once) and writes each to its own Parquet partition <out_dir>/sample_id=<id>/part.parquet, with the sample metadata sheet stored once as <out_dir>/metadata.parquet (not broadcast into every clonotype row).

scan_cohort() opens the whole cohort as a single polars.LazyFrame — sample_id recovered from the partition path, metadata joined lazily — so a cohort far larger than RAM is analysed with one streaming group_by:

>>> lf = scan_cohort("cohort/")
>>> usage = (lf.group_by(["sample_id", "v_call"])
...            .agg(pl.col("duplicate_count").sum())
...            .collect(engine="streaming"))

Every analysis in the package is already group_by(...).agg(...), so the cohort feature matrix is one streamed pass; only the final (bounded) sample×feature pivot is materialised.

vdjtools.io.cohort.SAMPLE_ID = 'sample_id'#: Reserved partition-key column recovered from the sample_id=<id> path.

vdjtools.io.cohort.ingest_cohort(metadata, base_dir, out_dir, sample_col='sample_name', file_template='{sample}.tsv.gz', fmt='auto')[source]#

Convert a sample cohort into a hive-partitioned Parquet dataset (streaming).

Reads each sample in turn via vdjtools.io.iter_samples() (O(1-sample) RAM) and writes the canonical clonotype frame to out_dir/sample_id=<id>/part.parquet. The metadata sheet is written once to out_dir/metadata.parquet with its sample column renamed to sample_id — it is not broadcast into the clonotype rows. Run this once; thereafter analyse the cohort with scan_cohort().

Parameters:

metadata (DataFrame) – Metadata frame (e.g. from vdjtools.io.read_metadata()).
base_dir (str | PathLike) – Directory holding the per-sample clonotype files.
out_dir (str | PathLike) – Destination directory for the partitioned dataset (created if absent).
sample_col (str) – Metadata column carrying the sample name.
file_template (str) – Filename template; {sample} is substituted per sample.
fmt (str) – Reader format ("auto" / "vdjtools" / "airr" / "parquet").

Returns:

The out_dir Path (feed straight into scan_cohort()).

Return type:

vdjtools.io.cohort.scan_cohort(out_dir, *, join_metadata=True)[source]#

Scan a hive-partitioned cohort (from ingest_cohort()) as one LazyFrame.

Nothing is read until .collect(); group_by / filter / column selection push down into the Parquet scan, so a cohort far larger than memory is reduced with .collect(engine="streaming"). sample_id is recovered from the partition path; the metadata columns are joined lazily on sample_id.

Parameters:

out_dir (str | PathLike) – The partitioned dataset directory written by ingest_cohort().
join_metadata (bool) – If True (default) and metadata.parquet is present, left-join the per-sample metadata onto the clonotype rows lazily.

Returns:

A polars.LazyFrame over the whole cohort (canonical clonotype columns + sample_id + metadata).

Return type:

LazyFrame

`vdjtools.io.convert`#

Converters for third-party repertoire formats (MiXcr, MiGec, Adaptive immunoSEQ, IMGT/HighV-QUEST, Vidjil, RTCR) to the canonical clonotype frame.

Converters for legacy repertoire input formats → the canonical clonotype frame.

Reimplements the battle-tested legacy vdjtools format parsers (their exact column mappings, gene-name normalisation and CDR3 handling were lifted verbatim from the Groovy com.antigenomics.vdjtools.io.parser classes) so third-party tool output can be read straight into the canonical AIRR-junction frame of vdjtools.io.schema:

MiXcr (v1/2 and v3/4 header dialects, incl. the C-gene / isotype hit) — read_mixcr()
MiGec — read_migec()
Adaptive immunoSEQ v1 and v2 — read_immunoseq()
IMGT/HighV-QUEST — read_imgt()
Vidjil (.vidjil JSON) — read_vidjil()
RTCR — read_rtcr()
TRUST4 (*_report.tsv) — read_trust4()
arda (AIRR annotation output) — read_arda()

Every reader returns the canonical frame: V/D/J IMGT calls, junction_nt / junction_aa (the AIRR junction — conserved Cys104 … Phe/Trp118 anchors included, matching the legacy vdjtools cdr3nt/cdr3aa), duplicate_count and recomputed frequency. Per-read formats are collapsed to unique clonotypes with summed counts.

vdjtools.io.convert.translate(seq)[source]#

Translate a (possibly out-of-frame) CDR3 nt sequence, V→J bidirectionally.

Port of the legacy CommonUtil.translate: in-frame sequences are a plain codon walk (stop → *); an out-of-frame sequence is padded in the middle with ? and translated inward from both ends, leaving the untranslatable middle codon(s) lower-cased (later collapsed to _ by to_unified_cdr3aa()).

Parameters:: seq (str)
Return type:: str

vdjtools.io.convert.to_unified_cdr3aa(seq)[source]#

Collapse each run of non-coding markers (lower-case nt / # ~ _ ?) to a single _.

Parameters:: seq (str | None)
Return type:: str | None

vdjtools.io.convert.extract_vdj(field)[source]#

First tie, allele stripped, quotes/space trimmed → an IMGT gene call (or None).

Port of CommonUtil.extractVDJ: "TRBV12-4,TRBV12-3" → "TRBV12-4", "TRBV13*00(401.5)" → "TRBV13". Empty / placeholder → None.

Parameters:: field (str | None)
Return type:: str | None

vdjtools.io.convert.read_mixcr(path, n_rows=None)[source]#

Read a MiXcr exportClones table — legacy (v1/2) or current (v3/4) dialect.

The two header dialects name every field differently, so both spellings are accepted for each column:

field	legacy (v1/2)	current (v3/4)
count	`Clone count`	`cloneCount` / `readCount` / `uniqueTagCountMolecule`
V/D/J/C hits	`All V hits` …	`allVHitsWithScore` …
CDR3 nt / aa	`N. Seq. CDR3` …	`nSeqCDR3` / `aaSeqCDR3`

MiXcr’s current field API exposes -readCount / -readFraction (there is no -cloneCount field any more), while the default preset still labels the column cloneCount — so a v4 export carries one spelling or the other depending on how it was produced, and both must parse.

Count precedence is read-based first (cloneCount → readCount), with the UMI molecule count (uniqueTagCountMolecule) used only when no read count column is present. On a UMI library the molecule count is the less PCR-biased abundance, so prefer exporting it alone if that is the quantity you want.

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

vdjtools.io.convert.read_migec(path, n_rows=None)[source]#

Read a MiGEC CdrBlast clonotype table.

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

vdjtools.io.convert.read_mitcr(path, n_rows=None)[source]#

Read a MiTCR / tcR (R package) clonotype table — the dot-separated dialect.

Header is Read.count Read.proportion CDR3.nucleotide.sequence CDR3.amino.acid.sequence V.gene J.gene D.gene V.end J.start D5.end D3.end VD.insertions DJ.insertions Total.insertions. Distinct from MiGEC’s CDR3 nucleotide sequence / V segments (spaces, not dots), so it needs its own picks. D.gene may carry an ambiguous call ("TRBD1, TRBD2"); extract_vdj() keeps the first.

Parameters:

path (str | PathLike) – Path to a MiTCR/tcR table.
n_rows (int | None) – Read only the first n_rows rows.

Returns:

Canonical clonotype frame.

Raises:

ValueError – If the signature columns are absent.

Return type:

DataFrame

vdjtools.io.convert.read_rtcr(path, n_rows=None)[source]#

Read an RTCR clonotype table (junction aa is re-translated from the nt junction).

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

vdjtools.io.convert.read_imgt(path, n_rows=None)[source]#

Read an IMGT/HighV-QUEST 1_Summary table (per-read → collapsed clonotypes).

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

vdjtools.io.convert.read_immunoseq(path, n_rows=None)[source]#

Read an Adaptive immunoSEQ export (v1 or v2 header dialect, auto-detected).

The Adaptive nomenclature (TCRBV29-01) is converted to IMGT (TRBV29-1) with the gene→family→family-ties fallback; the CDR3/junction nt is sliced out of the full rearrangement / nucleotide read via the vIndex + cdr3Length coordinates.

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

vdjtools.io.convert.read_vidjil(path, sample_id=0)[source]#

Read a Vidjil .vidjil JSON file.

Uses the anchor-inclusive seg.junction (never seg.cdr3, which excludes the anchors); the junction nt is sliced from the clone’s full sequence by the 1-based junction.start/junction.stop interval. sample_id selects the reads count for multi-sample files.

Parameters:

path (str | PathLike)
sample_id (int)

Return type:

DataFrame

vdjtools.io.convert.read_trust4(path, n_rows=None)[source]#

Read a TRUST4 clonotype report (*_report.tsv).

The TRUST4 report header is #count frequency CDR3nt CDR3aa V D J C cid cid_full_length. TRUST4’s CDR3 spans the conserved Cys104 … Phe/Trp118 anchors (anchors included), i.e. it is the AIRR junction, so CDR3nt/CDR3aa map straight to junction_nt/junction_aa (the _ stop / ? ambiguous-N markers are collapsed by to_unified_cdr3aa()). V/D/J/C keep the first allele’s gene (* → missing); the C column is the BCR isotype (or the TCR constant gene) when the constant region was captured. Rows whose CDR3 nt is not clean ACGT (TRUST4 partial / out_of_frame / N-containing) are dropped.

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

vdjtools.io.convert.read_arda(path, n_rows=None)[source]#

Read arda’s AIRR annotation output (per-sequence *.airr.tsv or clones.tsv).

arda (AIRR annotation + markup repair) writes standard AIRR Rearrangement column names, so this delegates to read_airr() — which maps v_call/d_call/j_call/c_call and junction/junction_aa and collapses reads to clonotypes — then nulls the literal "" arda emits for an empty gene call. arda’s extra columns (d2_call, c_class, mmseqs2_*) are ignored.

Parameters:

path (str | PathLike)
n_rows (int | None)

Return type:

DataFrame

Repertoire statistics (`vdjtools.stats`)#

Diversity estimators (observed, Chao1/ChaoE, Shannon/Simpson, d50, Efron-Thisted), iNEXT Hill-number rarefaction/extrapolation, similarity-aware (functional) diversity, spectratype, and V/J/VJ segment usage.

`vdjtools.stats.diversity`#

Point-estimate diversity indices.

Diversity estimators on a clonotype count vector.

Every estimator takes the per-clonotype duplicate_count vector of a single clonotype frame (one row per clonotype). Definitions follow legacy vdjtools (com.antigenomics.vdjtools.diversity); each deviation is noted in its docstring.

Notation: n = total reads, Sobs = observed richness (# clonotypes), f_i = count_i / n, F1 / F2 = number of singleton / doubleton clonotypes.

vdjtools.stats.diversity.observed_richness(counts)[source]#

Observed species richness Sobs — the number of clonotypes.

Parameters:: counts (ndarray) – Per-clonotype count vector.
Returns:: Number of clonotypes.
Return type:: int

vdjtools.stats.diversity.chao1(counts)[source]#

Chao1 lower-bound richness estimate.

Chao1 = Sobs + F1*(F1-1) / (2*(F2+1)) (bias-corrected form, legacy default).

Parameters:: counts (ndarray) – Per-clonotype count vector.
Returns:: The Chao1 richness estimate.
Return type:: float

vdjtools.stats.diversity.chao_e(counts, extrapolate_to=None)[source]#

Chao extrapolated richness at a larger sampling depth.

Follows legacy ChaoEstimator.chaoE: S = Sobs + F0 * [1 - (1 - F1/(n*F0))^(m* )] where F0 = F1*(F1-1)/(2*(F2+1)) and m* = extrapolate_to - n.

Parameters:

counts (ndarray) – Per-clonotype count vector.
extrapolate_to (int | None) – Target read depth (>= n). Defaults to 2*n — there is no cohort here to borrow “largest sample” from as legacy does, so the sample is extrapolated to double its size.

Returns:

The extrapolated richness. Falls back to Sobs when extrapolation is undefined (F0 == 0, i.e. too few singletons/doubletons).

Raises:

ValueError – If extrapolate_to < n.

Return type:

vdjtools.stats.diversity.efron_thisted(counts, max_depth=20, cv_threshold=0.05)[source]#

Efron–Thisted total-diversity lower-bound estimate.

Direct port of legacy ExactEstimator.getEfronThisted: increases the truncation depth (up to max_depth) until the coefficient of variation D/S reaches cv_threshold, accumulating alternating binomial coefficients over the count-frequency spectrum f_x (number of clonotypes seen exactly x times). The CV stopping rule is the only stopping rule — legacy has no “stop at the max observed count” cap.

Parameters:

counts (ndarray) – Per-clonotype count vector.
max_depth (int) – Maximum truncation depth.
cv_threshold (float) – Stop once std/estimate reaches this value.

Returns:

The Efron–Thisted richness estimate.

Return type:

vdjtools.stats.diversity.shannon_wiener(counts)[source]#

Shannon–Wiener diversity exp(H) (effective number of clonotypes).

Matches legacy shannonWienerIndex = exp(-Σ f_i ln f_i) (the Hill number of order 1), i.e. the exponential of Shannon entropy rather than the entropy itself.

Parameters:: counts (ndarray) – Per-clonotype count vector.
Returns:: exp(H).
Return type:: float

vdjtools.stats.diversity.normalized_shannon_wiener(counts)[source]#

Normalised Shannon–Wiener index H / ln(Sobs) (Pielou evenness).

Matches legacy normalizedShannonWienerIndex. This normalises the entropy H (not exp(H)) by ln(Sobs), so it lies in [0, 1].

Parameters:: counts (ndarray) – Per-clonotype count vector.
Returns:: H / ln(Sobs); 0.0 when Sobs <= 1 (no diversity to normalise).
Return type:: float

vdjtools.stats.diversity.inverse_simpson(counts)[source]#

Inverse Simpson index 1 / Σ f_i^2 (Hill number of order 2).

Parameters:: counts (ndarray) – Per-clonotype count vector.
Returns:: The inverse Simpson index.
Return type:: float

vdjtools.stats.diversity.d50(counts, fraction=0.5)[source]#

D50 — the legacy getDxxIndex dominance index 1 - k/Sobs.

Ranks clonotypes by descending count and finds the minimum number k whose cumulative read fraction reaches fraction; returns 1 - k/Sobs — the fraction of clonotypes not needed to cover fraction of the reads. This is the exact legacy ExactEstimator.getDxxIndex definition (1.0 - div / frequencyTable.diversity).

Parameters:

counts (ndarray) – Per-clonotype count vector.
fraction (float) – Cumulative-frequency target in [0, 1] (default 0.5).

Returns:

1 - k/Sobs in [0, 1); 0.0 when Sobs == 0.

Return type:

vdjtools.stats.diversity.diversity_stats(df, extrapolate_to=None)[source]#

Compute all diversity estimators for a clonotype frame as a one-row frame.

Parameters:

df (DataFrame) – A clonotype frame (one row per clonotype) with duplicate_count.
extrapolate_to (int | None) – Target depth for chaoE (see chao_e(); defaults to 2*n).

Returns:

A single-row pl.DataFrame with columns reads, observed_diversity, chao1, chaoE, efron_thisted, shannon_wiener, normalized_shannon_wiener, inverse_simpson, d50.

Return type:

DataFrame

vdjtools.stats.diversity.diversity_cohort(cohort, extrapolate_to=None, *, sample_col='sample_id')[source]#

Per-sample diversity_stats() for a whole cohort, low peak memory + exact.

Collapses the cohort to each sample’s count-frequency spectrum in one streamed pass — group_by([sample_id, duplicate_count]).agg(len) — which is bounded by the number of distinct counts per sample (tiny), never the number of clonotypes. Each sample’s count vector is then reconstructed from its spectrum and fed to the exact per-sample estimators, so every value is bit-identical to running diversity_stats() on the full sample (including Efron–Thisted / chaoE / d50, which need the spectrum/sorted form). Peak memory is O(Σ distinct counts), so a cohort far larger than RAM is handled with one collect().

Parameters:

cohort – A clonotype cohort with a sample_id column — a lazy vdjtools.io.scan_cohort() frame (streamed) or an eager pl.DataFrame.
extrapolate_to (int | None) – Target depth for chaoE (see diversity_stats()).
sample_col (str) – The per-sample id column (default "sample_id").

Returns:

sample_id followed by the diversity_stats() columns, in first-appearance order of sample_id.

Return type:

One row per sample

`vdjtools.stats.functional`#

Leinster-Cobbold similarity-aware diversity profile + Rao’s quadratic entropy.

Similarity-aware (functional) diversity of a single repertoire.

Classical diversity treats every clonotype as equally distinct: two repertoires of 1000 clones score the same whether those clones are near-identical variants of one motif or 1000 unrelated rearrangements. That is the wrong null for TCRs, where convergent recombination and clonal expansion both manufacture near-neighbours.

Leinster & Cobbold (2012, Ecology 93(3):477, doi:10.1890/10-2402.1) fix this by folding a similarity matrix Z into the Hill numbers. With p the clonotype abundance vector and (Zp)_i = Σ_j Z_ij p_j the ordinariness of clonotype i (how well-represented its neighbourhood is):

ᑫD^Z(p) = ( Σ_i p_i (Zp)_i^(q-1) )^(1/(1-q)) q ≠ 1 ¹D^Z(p) = exp( − Σ_i p_i ln (Zp)_i ) the q→1 limit Rao’s Q = Σ_ij p_i p_j (1 − Z_ij) = 1 − pᵀZp expected dissimilarity of two draws

Z = I (nothing resembles anything but itself) recovers the plain Hill numbers exactly — richness, exp(Shannon), inverse Simpson at q = 0, 1, 2. That identity is the defining special case and is what test_functional.py pins against vdjtools.stats.inext.

This is the single-community counterpart of vdjtools.overlap.similarity_overlap(), which applies the same kernel as a two-sample bilinear form pᵀZq. Both take Z from seqtree, so the kernel is defined in exactly one place.

q orders the profile by how much it cares about rare clonotypes: q=0 counts them fully (and is the most sensitive to sequencing depth), q=2 is dominated by the expanded ones. Read the profile, not a single number — that is the whole point of a Hill profile.

vdjtools.stats.functional.functional_diversity(df, *, q=(0, 1, 2), key=('junction_aa',), kernel='exp', tau=None, matrix=None, max_penalty=None, gap_prior='central', gap_open=None, dense=None, weight='freq', threads=0)[source]#

Leinster-Cobbold similarity-aware diversity profile of one repertoire.

Parameters:

df (DataFrame) – Clonotype frame (canonical schema), one row per clonotype.
q – Diversity orders to report. q=0 weights rare clonotypes fully, q=1 is the Shannon-equivalent, q=2 is dominated by expanded clones. Non-integer q works.
key (tuple[str, ...]) – Columns forming the clonotype identity; must include junction_aa.
kernel (str) – "exp" (Zᵢⱼ = exp(−Pᵢⱼ/τ)), "step" (1[P ≤ max_penalty]), or "identity" (Z = I — recovers the plain Hill numbers).
tau (float | None) – Kernel bandwidth for "exp".
matrix – Substitution matrix (default BLOSUM62, via seqtree).
max_penalty (int | None) – Alignment-penalty cutoff; the "step" kernel threshold.
gap_prior – Gap-placement prior passed to seqtree.
gap_open (int | None) – Gap-open penalty passed to seqtree.
dense (bool | None) – Force the dense (True) or sparse (False) kernel path.
weight (str) – "freq" (relative abundance) or "presence" (uniform — the unweighted, incidence-style profile).
threads (int) – Worker threads for the seqtree kernel build (0 = engine default).

Returns:

One row per requested q with columns q, diversity (ᑫD^Z) and rao (Rao’s quadratic entropy — constant across q, carried for convenience).

Raises:

ImportError – If seqtree (or, for the sparse path, scipy) is missing.
ValueError – On an unknown kernel/weight/gap_prior, or a key without junction_aa.

Return type:

DataFrame

Example

>>> functional_diversity(sample, q=(0, 1, 2))["diversity"].to_list()
[812.4, 402.1, 233.7]
>>> # Z = I is the sanity anchor: plain Hill numbers.
>>> functional_diversity(sample, q=(0,), kernel="identity")["diversity"][0]
4000.0

`vdjtools.stats.rarefaction`#

Rarefaction / extrapolation curves.

Rarefaction / extrapolation of repertoire diversity — the iNEXT framework.

rarefaction() is the single canonical entry point for rarefaction/extrapolation (R/E) of clonotype diversity. Its default (q=0, base="size") is the vdjtools-original richness R/E curve — the classic clonotype-accumulation curve, computed by the validated iNEXT q = 0 estimator. The same call also gives the full Hill-number profile (q) and coverage-based R/E (base="coverage").

The estimators themselves live in vdjtools.stats.inext (transcribed from the iNEXT R package and numerically validated against it). inext() and inext_coverage() are kept as aliases of the size- and coverage-based engines for discoverability; rarefaction() dispatches to them.

vdjtools.stats.rarefaction.rarefaction(data, q=0, base='size', *, sizes=None, coverages=None, endpoint=None, knots=40, se=True, nboot=50, conf=0.95, seed=0)[source]#

Rarefaction/extrapolation of Hill-number diversity (canonical entry point).

Single entry point for repertoire R/E, dispatching to the size- or coverage-based iNEXT engine (Chao et al. 2014; Hsieh et al. 2016). The default q=0, base="size" is the vdjtools-original richness rarefaction / extrapolation curve — the classic clonotype-accumulation curve, computed by the validated iNEXT q = 0 estimator. Pass a tuple of orders for the full Hill-number profile, or base="coverage" for coverage-based R/E.

Parameters:

data – A 1-D count vector (list/np.ndarray/pl.Series) of clonotype abundances, or a clonotype pl.DataFrame (duplicate_count used).
q – Hill order(s) — a scalar (default 0 = richness, the vdjtools-original curve) or a tuple such as (0, 1, 2) for the richness/Shannon/Simpson profile.
base – "size" for size-based R/E (uses sizes/endpoint/knots) or "coverage" for coverage-based R/E (uses coverages/knots).
sizes – Explicit sampling depths for base="size"; see inext().
coverages – Explicit target coverages for base="coverage"; see inext_coverage().
endpoint – Maximum depth for the default size grid (base="size").
knots – Number of grid points when sizes/coverages are None.
se – If True, compute bootstrap confidence intervals.
nboot – Number of bootstrap replicates.
conf – Confidence level for the intervals.
seed – Seed for the bootstrap RNG.

Returns:

columns order_q, m, method, sample_coverage, qD, qD_lo, qD_hi. For base="coverage" (see inext_coverage()): columns order_q, sample_coverage, m, method, qD, qD_lo, qD_hi.

Return type:

A tidy pl.DataFrame. For base="size" (see inext())

Raises:

ValueError – If base is not "size" or "coverage".

`vdjtools.stats.inext`#

iNEXT-style Hill-number rarefaction and extrapolation.

iNEXT interpolation/extrapolation of Hill-number diversity for a repertoire.

Size- and coverage-based rarefaction/extrapolation (R/E) plus asymptotic diversity estimators for Hill numbers of order q (0 = richness, 1 = Shannon, 2 = Simpson). The estimators are transcribed from the iNEXT R package (v3.0.2 internals TD.m.est, Chat.Ind, Diversity_profile, EstiBootComm.Ind, invChat.Ind) and the underlying papers, and are numerically validated against that package in tests/python/test_inext.py.

For extrapolation of orders q != 2 this uses iNEXT’s “beta” method (qD(n+m*) = D_obs + (D_asy - D_obs) * (1 - (1-beta)^m*)), which matches the R package rather than Eq. 10c of Chao et al. (2014). Order q = 2 uses the exact closed form for both interpolation and extrapolation. Only non-negative integer orders are supported (iNEXT’s default profile is q = (0, 1, 2)).

References

Chao, A., Gotelli, N. J., Hsieh, T. C., Sander, E. L., Ma, K. H., Colwell, R. K., & Ellison, A. M. (2014). Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecological Monographs, 84(1), 45-67. doi:10.1890/13-0133.1

Hsieh, T. C., Ma, K. H., & Chao, A. (2016). iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods in Ecology and Evolution, 7(12), 1451-1456. doi:10.1111/2041-210X.12613

vdjtools.stats.inext.inext(data, q=(0, 1, 2), *, sizes=None, endpoint=None, knots=40, se=True, nboot=50, conf=0.95, seed=0)[source]#

Size-based rarefaction/extrapolation of Hill-number diversity.

This is the size-based engine; the canonical entry point is vdjtools.stats.rarefaction() with base="size" — inext is kept as an alias for discoverability (people search for “iNEXT”).

Interpolates (m < n) and extrapolates (m > n) the Hill number of each order in q across a grid of sampling depths, following Chao et al. (2014) and the iNEXT R package (Hsieh et al. 2016). Point estimates are computed on the original data; confidence bands come from an augmented-assemblage bootstrap.

Parameters:

data – A 1-D count vector (list/np.ndarray/pl.Series) of clonotype abundances, or a clonotype pl.DataFrame (duplicate_count used).
q – Hill orders (non-negative integers); default (0, 1, 2).
sizes – Explicit sampling depths. If None, knots depths from 1 to endpoint are used, always including the observed depth n.
endpoint – Maximum depth when sizes is None. Defaults to 2*n.
knots – Number of depths when sizes is None.
se – If True, compute bootstrap confidence intervals.
nboot – Number of bootstrap replicates.
conf – Confidence level for the intervals.
seed – Seed for the bootstrap RNG.

Returns:

A tidy pl.DataFrame with columns order_q (Int64), m (Int64), method (rarefaction | observed | extrapolation), sample_coverage, qD, qD_lo, qD_hi (the CI columns are null when se=False).

Return type:

DataFrame

vdjtools.stats.inext.inext_coverage(data, q=(0, 1, 2), *, coverages=None, knots=40, se=True, nboot=50, conf=0.95, seed=0)[source]#

Coverage-based rarefaction/extrapolation of Hill-number diversity.

This is the coverage-based engine; the canonical entry point is vdjtools.stats.rarefaction() with base="coverage" — inext_coverage is kept as an alias for discoverability (people search for “iNEXT”).

Each target sample coverage is inverted to a sampling depth m (iNEXT invChat.Ind) and the Hill number is then evaluated at that depth.

Parameters:

data – Count vector or clonotype pl.DataFrame (see inext()).
q – Hill orders (non-negative integers); default (0, 1, 2).
coverages – Explicit target coverages in (0, 1). If None, a grid of the coverages attained by the default size grid (1..``2*n``) is used.
knots – Number of coverages when coverages is None.
se – If True, compute bootstrap confidence intervals.
nboot – Number of bootstrap replicates.
conf – Confidence level.
seed – Seed for the bootstrap RNG.

Returns:

A tidy pl.DataFrame with columns order_q, sample_coverage, m (Float64), method, qD, qD_lo, qD_hi.

Return type:

DataFrame

vdjtools.stats.inext.asymptotic_diversity(data, q=(0, 1, 2), *, se=True, nboot=50, conf=0.95, seed=0)[source]#

Asymptotic (estimated true) Hill-number diversity of each order.

Uses Chao1 (q = 0), the Chao-Wang-Jost Shannon Hill number (q = 1) and the MVUE Simpson Hill number (q = 2); see _asymptotic().

Parameters:

data – Count vector or clonotype pl.DataFrame (see inext()).
q – Hill orders (non-negative integers); default (0, 1, 2).
se – If True, compute bootstrap standard errors and intervals.
nboot – Number of bootstrap replicates.
conf – Confidence level.
seed – Seed for the bootstrap RNG.

Returns:

A pl.DataFrame with columns order_q, observed (plug-in), estimator (asymptotic), se, lo, hi.

Return type:

DataFrame

vdjtools.stats.inext.sample_coverage(data, m=None)[source]#

Estimated sample coverage Ĉ (Chao et al. 2014).

Parameters:

data – Count vector or clonotype pl.DataFrame (see inext()).
m – A single depth, an iterable of depths, or None (the observed depth n).

Returns:

Ĉ(n) (float) when m is None, Ĉ(m) (float) for a scalar m, or a pl.DataFrame with columns m and sample_coverage for an iterable m.

vdjtools.stats.inext.coverage(data, m=None)#: Clean top-level name for the sample-coverage curve Ĉ(n)/Ĉ(m).

vdjtools.stats.inext.estimate_d(data, base='size', level=None, q=(0, 1, 2), *, se=True, nboot=50, conf=0.95, seed=0)[source]#

Diversity at a target sample size or coverage (iNEXT estimateD analog).

Parameters:

data – Count vector or clonotype pl.DataFrame (see inext()).
base – "size" to fix a sampling depth, "coverage" to fix a coverage.
level – The depth (base="size"; defaults to the observed depth n) or the coverage (base="coverage"; required, in (0, 1)).
q – Hill orders (non-negative integers); default (0, 1, 2).
se – If True, compute bootstrap confidence intervals.
nboot – Number of bootstrap replicates.
conf – Confidence level.
seed – Seed for the bootstrap RNG.

Returns:

A pl.DataFrame with columns order_q, m (Float64), method, sample_coverage, qD, qD_lo, qD_hi.

Return type:

DataFrame

vdjtools.stats.inext.inext_batch(samples, q=(0, 1, 2), *, sizes=None, endpoint=None, knots=40, se=True, nboot=50, conf=0.95, seed=0, threads=0)[source]#

Size-based R/E of Hill-number diversity for many samples at once.

Computes the point curve and (optionally) bootstrap confidence intervals for every sample, parallelizing the per-sample work across a native thread pool (_core.inext_batch, GIL released). This is the “cohort of many repertoires, quickly” entry point; a single sample is better served by inext().

Parameters:

samples – A list of clonotype count vectors, or a clonotype pl.DataFrame with a sample_id column (grouped by it, weighted by duplicate_count).
q – Hill orders (non-negative integers); default (0, 1, 2).
sizes – Explicit sampling depths applied to every sample. If None, each sample gets knots depths from 1 to endpoint (its own default), always including its observed depth n.
endpoint – Maximum depth for the default per-sample grid. Defaults to 2*n.
knots – Number of depths when sizes is None.
se – If True, compute bootstrap confidence intervals.
nboot – Number of bootstrap replicates per sample.
conf – Confidence level for the intervals.
seed – Base RNG seed; sample i is seeded seed + i.
threads – Worker threads (0 = hardware_concurrency), capped at the number of samples.

Returns:

columns sample, order_q, m, method (rarefaction | observed | extrapolation), sample_coverage, qD, qD_lo, qD_hi (the CI columns are null when se=False).

Return type:

A tidy long pl.DataFrame with one block per sample

Raises:

ValueError – If a pl.DataFrame input lacks a sample_id column.

vdjtools.stats.inext.rarefaction_batch(samples, q=(0, 1, 2), *, sizes=None, endpoint=None, knots=40, se=True, nboot=50, conf=0.95, seed=0, threads=0)#

Alias of inext_batch() under the canonical rarefaction name.

Return type:: DataFrame

`vdjtools.stats.spectratype`#

CDR3 length spectratype (nt / aa, weighted / unweighted).

CDR3-length spectratypes (long-format polars).

The spectratype is the weighted distribution of CDR3 lengths. vj_spectratype additionally breaks the distribution down by V-J pairing (legacy SpectratypeV is the V-only special case, obtained by ignoring the j_call column).

vdjtools.stats.spectratype.spectratype(df, kind='aa', weight='reads', by_locus=True, by=())[source]#

CDR3-length distribution.

Parameters:

df – A clonotype frame (eager pl.DataFrame or a lazy pl.LazyFrame — e.g. a whole cohort from vdjtools.io.scan_cohort()); the result mirrors the input’s laziness.
kind (str) – "aa" (length of junction_aa) or "nt" (length of junction_nt).
weight (str) – "reads", "unique", or "freq".
by_locus (bool) – If True, break the distribution down per locus.
by – Extra column(s) to prepend to the group key (e.g. ["sample_id"] to compute the whole cohort in one grouped pass). Empty by default — output is byte-identical to the per-sample spectratype.

Returns:

Long-format frame with length, weight, the by columns and (if by_locus) locus, sorted by the group key; a pl.LazyFrame when df is lazy.

vdjtools.stats.spectratype.vj_spectratype(df, kind='aa', weight='reads', by_locus=True, keep_allele=False, by=())[source]#

CDR3-length distribution broken down by V-J pairing.

Parameters:

df – A clonotype frame (eager or lazy; see spectratype()).
kind (str) – "aa" or "nt" (see spectratype()).
weight (str) – "reads", "unique", or "freq".
by_locus (bool) – If True, break the distribution down per locus.
keep_allele (bool) – If True, keep allele suffixes; otherwise collapse V/J to gene level (default).
by – Extra column(s) prepended to the group key (e.g. ["sample_id"] for a one-pass cohort profile).

Returns:

Long-format frame with v_call, j_call, length, weight, the by columns and (if by_locus) locus (lazy when df is lazy). A V-only spectratype is obtained by summing weight over j_call.

`vdjtools.stats.usage`#

V / J / VJ segment-usage vectors and matrices.

V / D / J / C segment usage and V-J pairing profiles (long-format polars).

Usage is the summed weight per segment: reads (duplicate_count) or unique clonotypes. Allele suffixes (*01) are stripped to gene level by default. Results are long-format and normalisable to fractions by dividing within a locus.

vdjtools.stats.usage.segment_usage(df, segment, weight='reads', by_locus=True, keep_allele=False, by=())[source]#

Segment usage profile for one of the V, D, J, or C genes.

Parameters:

df – A clonotype frame (eager pl.DataFrame or a lazy pl.LazyFrame — e.g. a whole cohort from vdjtools.io.scan_cohort()); the result mirrors the input’s laziness.
segment (str) – One of "v", "d", "j", "c".
weight (str) – "reads" (sum duplicate_count), "unique" (count clonotypes), or "freq" (sum frequency).
by_locus (bool) – If True, break the profile down per locus.
keep_allele (bool) – If True, keep the IMGT allele suffix; otherwise collapse to gene level (default).
by – Extra column(s) to prepend to the group key (e.g. ["sample_id"] to compute the whole cohort in one grouped pass over a LazyFrame). Empty by default — output is byte-identical to the per-sample profile.

Returns:

Long-format frame with the segment call column (v_call …), a weight column, the by columns, and (if by_locus) a locus column; a pl.LazyFrame when df is lazy. Rows with a null segment call are dropped, so an all-null C gene yields an empty frame.

Raises:

ValueError – If segment is not one of v/d/j/c.

vdjtools.stats.usage.vj_usage(df, weight='reads', by_locus=True, keep_allele=False, by=())[source]#

V-J pairing usage profile.

Parameters:

df – A clonotype frame (eager or lazy; see segment_usage()).
weight (str) – "reads", "unique", or "freq" (see segment_usage()).
by_locus (bool) – If True, break the profile down per locus.
keep_allele (bool) – If True, keep allele suffixes; otherwise collapse to gene level (default).
by – Extra column(s) prepended to the group key (e.g. ["sample_id"] for a one-pass cohort profile).

Returns:

Long-format frame with v_call, j_call, weight, the by columns and (if by_locus) locus (lazy when df is lazy). Rows with a null V or J call are dropped.

CDR3 features (`vdjtools.features`)#

Per-clonotype and sample-level CDR3 sequence features: amino-acid physicochemical region profiles and k-mer / V+k-mer summaries.

`vdjtools.features.physchem`#

Physicochemical CDR3 region profiles (Kidera factors, charge, hydropathy, …).

CDR3 physicochemical-property profiles from the legacy amino-acid property table.

For each clonotype the mean of a property over the residues of a chosen CDR3 region is computed, then averaged (weighted by reads/frequency, or unweighted) within a group (e.g. per V-J pairing, or per locus).

Region definitions (over junction_aa, length L):

all — the entire CDR3.
trimmed — junction_aa[3:-3] (conserved-anchor-trimmed core); clonotypes with L <= 6 have an empty core and are skipped for this region.
center — the middle five residues junction_aa[L//2-2 : L//2+3]; clonotypes with L < 5 are skipped for this region.

vdjtools.features.physchem.DEFAULT_PROPERTIES = ('hydropathy', 'charge', 'polarity', 'volume', 'strength', 'kf1', 'kf2', 'kf3', 'kf4', 'kf5', 'kf6', 'kf7', 'kf8', 'kf9', 'kf10')#

Kidera-factor-free physicochemistry + the 10 Kidera factors.

Type:: Default property subset

vdjtools.features.physchem.load_property_table()[source]#

Load the legacy amino-acid property table (cached).

The shipped resources/aa_property_table.txt uses classic-Mac \r line endings and a leading ## reference comment; both are handled here.

Returns:: A pl.DataFrame with an amino_acid column and one column per property (all property columns cast to Float64).
Return type:: DataFrame

vdjtools.features.physchem.physchem_profile(df, group_by=('v_call', 'j_call'), region='all', weight='reads', properties=None, by=())[source]#

Group-wise weighted mean of CDR3 physicochemical properties.

For each clonotype the region residues are looked up in the property table and averaged per property (the per-clonotype property mean). These are then combined per group by a weighted mean Σ w_c m_c / Σ w_c (w_c = clonotype weight).

Parameters:

df – A clonotype frame (eager pl.DataFrame or lazy pl.LazyFrame — e.g. a whole cohort from vdjtools.io.scan_cohort()); result mirrors the input.
group_by – Grouping column(s). Either the string "locus" (derived if absent) or an iterable of column names present in df (default ("v_call", "j_call")).
region (str) – "all", "trimmed", or "center" (see module docstring).
weight (str) – "reads", "unique", or "freq".
properties (tuple[str, ...] | None) – Property names to compute. Defaults to DEFAULT_PROPERTIES.
by – Extra column(s) to prepend to the group key (e.g. ["sample_id"] to profile a whole cohort in one pass). Empty by default — byte-identical to the per-sample profile.

Returns:

Tidy frame with the by columns, the group columns, region, property and mean_value, sorted by group then property (lazy when df is lazy).

Raises:

ValueError – If region is unknown or a requested property is missing.

`vdjtools.features.kmer`#

k-mer and V+k-mer occurrence summaries.

CDR3 amino-acid k-mer profiles and joint V + k-mer + C feature summaries.

K-mers are overlapping sliding windows over junction_aa. Each k-mer occurrence carries its clonotype’s weight (reads, unique, or frequency); weights are summed.

Anchors carry no information. A junction begins with the conserved Cys104 and a few germline V residues and ends with germline J: on a 250k human-TRB control, an N-terminal 4-mer is shared by 31.0% of clonotypes (CASS alone by 56.5%) while a central 4-mer is shared by 0.080% — ~386x more selective (measured in seqtree.seeds). So a window that starts at index 0 spends most of its multiple-testing budget on germline. Pass flank to drop that many residues from each end first; flank=4 matches seqtree.seeds.core_kmers.

flank is a fixed trim, so it is only an approximation of the germline boundary — the real one is per-clonotype and locus-specific (TRBJ2-3’s germline is STDTQYF, seven residues, so flank=4 still leaves STD in the core). Where the data carries V/J markup, prefer it.

vdjtools.features.kmer.kmer_profile(df, k=3, weight='reads', by_locus=True, flank=0, by=())[source]#

CDR3 amino-acid k-mer spectrum.

Parameters:

df – A clonotype frame (eager pl.DataFrame or lazy pl.LazyFrame — e.g. a whole cohort from vdjtools.io.scan_cohort()); result mirrors the input.
k (int) – K-mer length (default 3).
weight (str) – "reads", "unique", or "freq".
by_locus (bool) – If True, break the spectrum down per locus.
flank (int) – Residues to drop from each end before windowing (see the module docstring); 0 (default) keeps the whole junction, 4 gives the seqtree core.
by – Extra column(s) to prepend to the group key (e.g. ["sample_id"] for a one-pass cohort spectrum). Empty by default — byte-identical to per-sample.

Returns:

Long-format frame with kmer, weight, the by columns and (if by_locus) locus, sorted by the group key (lazy when df is lazy).

vdjtools.features.kmer.v_kmer_c_profile(df, k=3, weight='reads', by_locus=True, keep_allele=False, flank=0, by=())[source]#

Joint (V gene, k-mer, C gene) profile — a tidy feature-matrix source.

Produces one aggregated row per (v_call, kmer, c_call) combination (plus locus if requested), suitable to pivot into a feature matrix. A null c_call (common in native vdjtools data) is retained as its own group.

Parameters:

df – A clonotype frame (eager or lazy; see kmer_profile()).
k (int) – K-mer length (default 3).
weight (str) – "reads", "unique", or "freq".
by_locus (bool) – If True, include locus in the grouping.
keep_allele (bool) – If True, keep V allele suffixes; otherwise collapse V to gene level (default).
flank (int) – Residues to drop from each end before windowing (see the module docstring).
by – Extra column(s) prepended to the group key (e.g. ["sample_id"] for a one-pass cohort feature source).

Returns:

Long-format frame with v_call, kmer, c_call, weight, the by columns and (if by_locus) locus (lazy when df is lazy).

vdjtools.features.kmer.kmer_cohort(cohort, k=4, flank=4, keep_allele=False)[source]#

Explode a cohort into a per-(sample_id, v_call, kmer) frame for an incidence test.

The bridge from k-mers to vdjtools.biomarker.association(): the result carries a kmer column alongside sample_id, so a V + k-mer phenotype test is just

from vdjtools.biomarker import association
from vdjtools.features.kmer import kmer_cohort
km = kmer_cohort(cohort, k=4, flank=4)
res = association(km, design, key=("v_call", "kmer"), match="exact")

match="exact" is the only valid mode here — fuzzy/1mm search on the CDR3, and a 1-mismatch ball around a 4-mer is most of k-mer space.

Pinning V matters. A short central k-mer is only selective given the germline context: on its own a 4-mer is shared by ~0.08% of a control repertoire, which over ~10⁵ features is still thousands of hits. (v_call, kmer) is the feature the V+k-mer search is named for.

Parameters:

cohort (DataFrame | LazyFrame) – A clonotype frame with sample_id (e.g. from vdjtools.io.scan_cohort()).
k (int) – K-mer length. Note seqtree.seeds’ measurement that a central k-mer’s median E-value crosses 1 at k=6; below that, prune with seqtree.seeds.SeedIndex against a real control rather than trusting a modelled background (its residual KL grows with k — D-gene germline runs correlate — so it must be counted, not fitted).
flank (int) – Residues dropped from each end (see the module docstring). 4 = seqtree’s core.
keep_allele (bool) – Keep V allele suffixes; default collapses to gene level.

Returns:

Unique (sample_id, v_call, kmer) rows.

Return type:

DataFrame

Repertoire overlap and TCRnet (`vdjtools.overlap`)#

Sample overlap and TCRnet built on vdjmatch / seqtree, including sequence-similarity-aware (TINA / Leinster-Cobbold) overlap, pairwise-distance matrices, clustering / MDS, and tracking.

`vdjtools.overlap.metrics`#

Exact-match overlap metrics (F, D, Jaccard, Morisita-Horn, …).

Exact-match pairwise repertoire-overlap metrics (pure polars + numpy).

Implements the four legacy vdjtools overlap metrics (OverlapEvaluator) on an exact clonotype match key. This is exact-match only; fuzzy / e-value overlap and TCRnet are delegated to vdjmatch (cluster.overlap / evalue.query_evalues).

vdjtools.overlap.metrics.DEFAULT_KEY = ('junction_aa', 'v_call', 'j_call')#

CDR3 aa + V + J).

Type:: Default exact match key (legacy “Strict” overlap

vdjtools.overlap.metrics.STANDARD_AA = '^[ACDEFGHIKLMNPQRSTVWY]+$'#: CDR3s outside this alphabet (*, _, X) are not scoreable by the similarity kernels.

vdjtools.overlap.metrics.overlap_pair(a, b, key=('junction_aa', 'v_call', 'j_call'))[source]#

Compute the shared-clonotype table and overlap metrics for two samples.

Both frames are first collapsed to unique clonotype keys (summing counts; frequencies are recomputed within each sample). The four legacy metrics are:

D (diversity): d12 / (d1 * d2).
F (frequency): sqrt(Σ_shared f_a · Σ_shared f_b).
F2: Σ_shared sqrt(f_a · f_b).
R: Pearson correlation of the raw shared-clonotype frequencies (legacy OverlapEvaluator: x[k] = it.getFreq(i) fed straight to PearsonsCorrelation, no log transform). Requires at least three shared clonotypes (legacy guard n > 2); with fewer, or when the correlation is undefined (a constant/degenerate vector yields NaN), R is None. Legacy coerced both of these cases to 0.0; None keeps “undefined” distinct from a genuine zero correlation.

Parameters:

a (DataFrame) – First clonotype frame.
b (DataFrame) – Second clonotype frame.
key (tuple[str, ...]) – Columns forming the exact match key (default ("junction_aa", "v_call", "j_call"); use ("junction_aa",) for CDR3-only or add "junction_nt" for nucleotide-level matching).

Returns:

A tuple (shared, metrics) where shared is a pl.DataFrame of the joined shared clonotypes (key columns plus count_a, count_b, freq_a, freq_b) and metrics is a dict with keys D, F, F2, R, d1, d2, d12 (R is None when undefined; see above).

Return type:

tuple[DataFrame, dict]

vdjtools.overlap.metrics.overlap_metrics(a, b, key=('junction_aa', 'v_call', 'j_call'))[source]#

Compute the four exact-match overlap metrics (D, F, F2, R) for two samples.

Parameters:

a (DataFrame) – First clonotype frame.
b (DataFrame) – Second clonotype frame.
key (tuple[str, ...]) – Exact match key (see overlap_pair()).

Returns:

Dict with keys D, F, F2, R, d1, d2, d12.

Return type:

`vdjtools.overlap.similarity`#

Sequence-similarity-weighted (TINA / Leinster-Cobbold) overlap.

Sequence-similarity-weighted repertoire overlap (TINA / Leinster-Cobbold).

Where vdjtools.overlap.metrics scores overlap on an exact clonotype match and vdjtools.overlap.fuzzy on a within-edit-scope match, this module scores it through a continuous CDR3 similarity kernel Z (Leinster-Cobbold / Nei; Schmidt et al. 2016 TINA). Two repertoires are relative-abundance vectors p, q over their clonotypes and the overlap is built on the bilinear form pᵀZq = Σᵢⱼ pᵢ Zᵢⱼ qⱼ:

cosine (TINA_w): S = pᵀZq / sqrt((pᵀZp)(qᵀZq))
Morisita-Horn: S = 2·pᵀZq / (pᵀZp + qᵀZq)
distance = 1 − S.

Z is a kernel of a seqtree gap-block alignment penalty P (≥0, 0 identical), symmetric, with Zᵢᵢ = 1:

kernel="exp" — Zᵢⱼ = exp(−Pᵢⱼ/τ) (Leinster-Cobbold / Nei); τ defaults to SubstitutionMatrix.blosum62().scale() (= 14).
kernel="step" — Zᵢⱼ = 1[Pᵢⱼ ≤ max_penalty]. On unit cost with indels prohibited this is exactly vdjmatch’s fuzzy edit-distance overlap.
kernel="identity" — Z = I (match on the clonotype key). This recovers the exact frequency overlap: cosine collapses to the classical cosine of the shared-frequency vectors and Morisita to the classical Morisita-Horn index.

The three kernels are the same code on the same spine — identity and step are the exact special cases (exact / fuzzy overlap) of the continuous exp form. The CDR3 kernel is block-diagonalised by the non-cdr3 key fields (V/J/locus): exp/step never connect clonotypes that differ outside the CDR3, so identity is the exact τ→0 limit of exp on the same key.

The penalty is built with seqtree (a base dependency): the dense path scores every clonotype pair via seqtree.score_matrix() (O(N²), for small N and the tests); the sparse path uses seqtree.pairwise_batch() only to find near candidates, then re-scores each with the same gap-block model (seqtree.gapblock_score(), same matrix / gap_open / gap_prior as dense) and keeps those with penalty ≤ the threshold — so dense and sparse agree on every retained pair. It assembles a scipy.sparse Z (the practical path at scale). Both within-sample blocks (Z_AA, Z_BB) are needed in full — a diagonal-only approximation of pᵀZp is wrong.

class vdjtools.overlap.similarity.SimilarityMatrices(z_ab, z_aa, z_bb, keys_a, keys_b, freq_a, freq_b, kernel, tau, sparse)[source]#

Bases: object

The three similarity blocks plus the aligned weight vectors for a sample pair.

Variables:

z_ab (Any) – Cross-sample kernel, shape (n_a, n_b) (dense np.ndarray or a scipy.sparse matrix).
z_aa (Any) – Within-a kernel, shape (n_a, n_a); diagonal 1, symmetric.
z_bb (Any) – Within-b kernel, shape (n_b, n_b).
keys_a (list) – Clonotype key tuples for a in row order.
keys_b (list) – Clonotype key tuples for b in row order.
freq_a (numpy.ndarray) – Within-sample relative abundance of each a clonotype (sums to 1).
freq_b (numpy.ndarray) – Within-sample relative abundance of each b clonotype.
kernel (str) – The kernel used ("exp", "step", "identity").
tau (float | None) – The kernel bandwidth τ (None for non-exp kernels).
sparse (bool) – Whether the blocks are scipy.sparse matrices.

Parameters:

z_ab (Any)
z_aa (Any)
z_bb (Any)
keys_a (list)
keys_b (list)
freq_a (ndarray)
freq_b (ndarray)
kernel (str)
tau (float | None)
sparse (bool)

z_ab: Any#

z_aa: Any#

z_bb: Any#

keys_a: list#

keys_b: list#

freq_a: ndarray#

freq_b: ndarray#

kernel: str#

tau: float | None#

sparse: bool#

vdjtools.overlap.similarity.similarity_matrix(a, b, *, key=('junction_aa',), kernel='exp', tau=None, matrix=None, max_penalty=None, gap_prior='central', gap_open=None, dense=None, threads=0)[source]#

Build the three CDR3-similarity blocks Z_AB, Z_AA, Z_BB for a sample pair.

Parameters:

a (DataFrame) – First clonotype frame (canonical schema).
b (DataFrame) – Second clonotype frame.
key (tuple[str, ...]) – Columns forming the clonotype identity; must include junction_aa (the similarity match unit). Default ("junction_aa",).
kernel (str) – "exp" (exp(−P/τ)), "step" (1[P ≤ max_penalty]), or "identity" (Z = I on the key — exact overlap, no seqtree).
tau (float | None) – Kernel bandwidth for "exp"; defaults to the matrix scale (blosum62().scale() == 14).
matrix – A seqtree.SubstitutionMatrix; defaults to BLOSUM62 for "exp" and to unit cost (edit-count penalty) for "step".
max_penalty (int | None) – Penalty ceiling. Required (and the step threshold) for "step" (default 1 = one edit); for the sparse "exp" path it is the neighbourhood cutoff (default from τ so exp(−P/τ) ≳ 1e-3).
gap_prior – Single-indel block placement rule: "central" (default — biases the indel to the loop centre, where CDR3 length variation sits; seqtree’s own gapblock recommendation for pairwise scoring), "frame" (left-anchored common frame), "none", or a seqtree.GapPrior.
gap_open (int | None) – Block-opening cost. Defaults to the seqtree score default for "exp" and to a gap-prohibiting value for "step" (substitution-only, matching vdjmatch’s default scope).
dense (bool | None) – Force the dense (True) or sparse (False) path; None auto-selects (dense while max(n_a, n_b) ≤ 1500).
threads (int) – Worker threads for the native search (0 = all cores).

Returns:

A SimilarityMatrices.

Raises:

ImportError – If seqtree (or, for the sparse path, scipy) is missing.
ValueError – On an unknown kernel/gap_prior or a key without junction_aa.

Return type:

SimilarityMatrices

vdjtools.overlap.similarity.similarity_overlap(a, b, *, key=('junction_aa',), metric='cosine', weight='freq', kernel='exp', tau=None, matrix=None, max_penalty=None, dense=None, threads=0)[source]#

Similarity-weighted overlap between two repertoires.

Builds Z via similarity_matrix(), forms the weight vectors p, q, and returns the cosine (TINA_w) or Morisita-Horn similarity on pᵀZq.

Parameters:

a (DataFrame) – First clonotype frame.
b (DataFrame) – Second clonotype frame.
key (tuple[str, ...]) – Clonotype identity key (must include junction_aa).
metric (str) – "cosine" (pᵀZq / sqrt(pᵀZp·qᵀZq)) or "morisita" (2·pᵀZq / (pᵀZp + qᵀZq)).
weight (str) – "freq" (relative abundance) or "presence" (uniform per clonotype, TINA-unweighted).
kernel (str) – "exp", "step", or "identity" (see similarity_matrix()). "identity" recovers the classical (exact) cosine / Morisita-Horn.
tau (float | None) – Kernel bandwidth for "exp" (default 14).
matrix – Substitution matrix (see similarity_matrix()).
max_penalty (int | None) – Penalty ceiling / step threshold (see similarity_matrix()).
dense (bool | None) – Force the dense (True) or sparse (False) kernel path; None auto-selects (see similarity_matrix()). Both agree on retained pairs.
threads (int) – Worker threads for the native search.

Returns:

Dict with keys similarity, distance, pTZq, pTZp, qTZq, metric, kernel.

Raises:

ValueError – On an unknown metric/weight/kernel.

Return type:

`vdjtools.overlap.fuzzy`#

Fuzzy (mismatch-tolerant) overlap via vdjmatch / seqtree.

Fuzzy (edit-distance) sample overlap — delegated to vdjmatch/seqtree.

Where vdjtools.overlap.metrics matches clonotypes exactly, this module matches CDR3s within an edit-distance scope (substitutions/indels), so a pair of repertoires that share only near-variants of a clonotype still register as overlapping. The fuzzy search itself is not reimplemented here — it is delegated to vdjmatch.cluster.overlap(), which runs on the native seqtree engine. These functions are thin polars wrappers that (a) collapse each sample to unique CDR3s, (b) hand the CDR3 lists to vdjmatch, and (c) join the matched pairs back to per-clonotype counts/frequencies.

Matching is on the amino-acid CDR3 (junction_aa): vdjmatch runs on the "aa" alphabet, so this module is aa-only. Scope syntax is vdjmatch’s "subs,ins,dels,total" (max substitutions, insertions, deletions, total edits); the default "1,0,0,1" is a single substitution.

vdjtools.overlap.fuzzy.fuzzy_overlap(a, b, scope='1,0,0,1', threads=0)[source]#

Fuzzy-matched clonotype pairs between two samples (within an edit scope).

Both samples are collapsed to unique junction_aa values (counts summed, frequencies recomputed within each sample); the two CDR3 lists are handed to vdjmatch.cluster.overlap(), and the returned within-scope pairs are joined back to per-clonotype counts and frequencies.

Parameters:

a (DataFrame) – First clonotype frame (canonical schema).
b (DataFrame) – Second clonotype frame.
scope (str) – vdjmatch edit-distance scope "subs,ins,dels,total" (default one substitution).
threads (int) – Worker threads for the native search (0 = all cores).

Returns:

A pl.DataFrame, one row per matched pair, with columns a_cdr3, b_cdr3, n_subs, score, count_a, freq_a, count_b, freq_b. Empty (with that schema) when nothing matches.

Raises:

ImportError – If vdjmatch is not importable (it is a base dependency).

Return type:

DataFrame

vdjtools.overlap.fuzzy.fuzzy_overlap_metrics(a, b, scope='1,0,0,1', threads=0)[source]#

Summary fuzzy-overlap metrics for two samples.

Computes fuzzy_overlap() once and derives:

pairs — number of within-scope matched clonotype pairs.
frac_a_matched / frac_b_matched — fraction of each sample’s unique clonotypes with at least one fuzzy match in the other.
fuzzy_F — the frequency-weighted fuzzy analogue of the exact F metric, sqrt(Σ_{a matched} freq_a · Σ_{b matched} freq_b), where a clonotype’s frequency is counted once if it has any fuzzy neighbour in the other sample (so a clonotype matching several near-variants is not double-weighted).

Parameters:

a (DataFrame) – First clonotype frame.
b (DataFrame) – Second clonotype frame.
scope (str) – vdjmatch edit-distance scope (see fuzzy_overlap()).
threads (int) – Worker threads for the native search.

Returns:

Dict with keys pairs, frac_a_matched, frac_b_matched, fuzzy_F.

Raises:

ImportError – If vdjmatch is not importable (it is a base dependency).

Return type:

`vdjtools.overlap.tcrnet`#

TCRnet neighbourhood-enrichment degree statistics.

TCRnet neighbourhood-enrichment (convergence) test — delegated to vdjmatch/seqtree.

The classic TCRnet / CalcDegreeStats analysis asks, for every clonotype in a sample, whether it has more close CDR3 neighbours than a background/generative process would produce — the signature of antigen-driven convergent selection. The legacy tool counted a clonotype’s within-sample degree and compared it to a control sample’s degree under a grouping (V/VJ/VJL) that bounded the comparison scope.

This reimplements it on a finite-sample, control-calibrated footing by delegating to vdjmatch.evalue.query_evalues() (the seqtree E-value engine):

target = a fuzzy seqtree.Index over the sample’s own unique CDR3s — so a clonotype’s neighbour count is its within-sample degree;
control = a matched background repertoire (vdjmatch.evalue.background());
for each query CDR3, E = (N/M)·n_control — where N is the target index size and M the control index size — is the neighbour count expected from the background, and p_enrichment is the Poisson tail probability of seeing at least the observed within-sample degree by chance. q_value is the Benjamini-Hochberg FDR over the clonotypes scored in one call.

Self-hit handling (this is subtle and was previously wrong). The query is a member of the target index but is not a member of the control, so the two sides need different treatment — and query_evalues(exclude_exact=True) punctures distance-0 hits on both. That silently dropped genuine public-clone neighbours from the background, understating E and inflating significance. We therefore query with exclude_exact=False and subtract the self-hit from the target side only: verified, n_target differs by exactly 1 between the two modes for every query (the query’s own copy), while n_control differs by 1 for the public clonotypes whose background neighbours were being discarded.

Locus handling. A TCRnet neighbourhood is only meaningful within one locus (a TRA sequence has no true neighbours in a TRB background). When neither control nor locus is given, the sample is split by the locus of its v_call and each locus is scored against its own matched background — mirroring the legacy V-grouping that kept loci disjoint. Clonotypes whose locus cannot be resolved (null v_call) are dropped with a warning. Passing control or locus explicitly overrides this and scores every clonotype against that single background.

Scope note: the legacy default was s,id,t = 1,0,1 (one substitution, no indels, one total edit), i.e. vdjmatch "1,0,0,1" — the default here.

vdjtools.overlap.tcrnet.tcrnet(sample, control=None, scope='1,0,0,1', locus=None, species='human', exclude_exact=True, threads=0)[source]#

Per-clonotype neighbourhood-enrichment (TCRnet) test for one sample.

Collapses the sample to unique CDR3s, builds a fuzzy seqtree.Index over them (the within-sample neighbourhood target), and queries each CDR3 against that target and a background control via vdjmatch.evalue.query_evalues(). When neither control nor locus is given the sample is scored per locus, each against its own matched background (see the module docstring).

Parameters:

sample (DataFrame) – Clonotype frame (canonical schema).
control – A prebuilt seqtree.Index background, or None to load matched background(s) via vdjmatch.evalue.background(). Passing one overrides the per-locus split (every clonotype is scored against it).
scope (str) – vdjmatch edit-distance scope "subs,ins,dels,total" defining the neighbourhood ball (default one substitution).
locus (str | None) – Force a single background locus (e.g. "TRB"); overrides the per-locus split. When None (and control is None) the sample is split by the locus of its v_call.
species (str) – Species for the background control (default "human").
exclude_exact (bool) – Whether a clonotype’s own copy counts toward its within-sample degree (default True = it does not). This applies to the target side only — the control side is never punctured, because the query is not a member of the control and a distance-0 control hit is a genuine public-clone background neighbour. Set False for the closed-ball convention ALICE uses.
threads (int) – Worker threads for the native search (0 = all cores).

Returns:

One row per unique clonotype (per locus) with columns junction_aa, v_call, j_call, duplicate_count, n_neighbors (within-sample degree), n_control (background degree), E (expected neighbours), p_enrichment (Poisson tail), q_value (Benjamini-Hochberg FDR over the clonotypes scored in this call), p_any, and locus — sorted by ascending p_enrichment.

Raises:

ImportError – If vdjmatch is not importable (it is a base dependency).
ValueError – If control and locus are None and no clonotype has a resolvable locus.

Return type:

DataFrame

`vdjtools.overlap.alice`#

The same neighbourhood-enrichment test against a V(D)J generation model (Pgen null) rather than a control repertoire — the complement of tcrnet.

ALICE neighbourhood enrichment — the generative null, against TCRnet’s empirical one.

ALICE (Pogorelyy et al., PLoS Biol 2019) asks the same question as vdjtools.overlap.tcrnet — does this clonotype have more close CDR3 neighbours than chance? — but answers it from a V(D)J generation model rather than from a control repertoire. For each amino-acid CDR3 σ within one V–J class:

λ_σ  =  n · Σ_{σ' : Hamming(σ,σ') ≤ 1}  Q · Pgen(σ')
p    =  P(Poisson(λ_σ) ≥ d(σ))

where n is the number of unique nucleotide clonotypes in that V–J class, d(σ) counts the nucleotide clonotypes whose CDR3aa lies within one mismatch of σ (σ included — different nucleotide variants of one amino-acid sequence are genuine neighbours), and Q rescales for thymic selection, which removes a fraction 1 − 1/Q of generated sequences.

The two nulls are complements, not rivals, and that is why both are here: ALICE’s generative null controls for the intrinsic biases of V(D)J recombination but knows nothing about selection or about which clonotypes are already common in people; TCRnet’s control repertoire absorbs thymic selection and endemic-pathogen expansions but needs a large, HLA-matched cohort to do it.

Both share a known blind spot worth stating plainly: neighbourhood enrichment cannot see a monoclonal expansion. A single hyperexpanded clone has no near neighbours by definition, so it scores as unremarkable no matter how dominant it is. That is precisely why this module and vdjtools.dynamics are complementary — enrichment measures breadth, the paired test measures magnitude.

The Hamming-1 ball sum is exact and already native: pgen_aa_batch(..., mismatches=1) computes Σ_k Pgen(a_{k→*}) − (L−1)Pgen(a), which is the closed ball with the centre counted exactly once — verified here against brute-force enumeration of all 19L neighbours to 2e-16.

vdjtools.overlap.alice.DEFAULT_Q = 9.41#: ALICE’s default thymic-selection factor — the paper’s average over V–J combinations.

vdjtools.overlap.alice.alice(sample, model=None, *, locus=None, source='olga', scope='1,0,0,1', selection_q=9.41, min_degree=3, min_count=2, threads=0)[source]#

Per-clonotype neighbourhood enrichment against a V(D)J generation model.

Parameters:

sample (DataFrame) – Clonotype frame (canonical schema). junction_nt is used when present so degrees are in nucleotide-clonotype units, as ALICE defines them; without it each amino-acid sequence counts once.
model – A Model. None loads the bundled model for locus.
locus (str | None) – Locus to score (e.g. "TRB"). None infers it from v_call; a sample spanning several loci is scored per locus.
source (str) – Bundled model source. Pinned to ``”olga”`` by default and you should leave it there: the "learned" models are EM-fit on ~2k clonotypes with no gene-usage pseudocount, so 68 of 89 bundled TRB V alleles have P(V) = 0 (vs 8 for OLGA) — those clonotypes get λ = 0 and come back infinitely significant. Their ball-Pgen scale is also ~16x below OLGA’s, which Q is calibrated against.
scope (str) – vdjmatch edit scope defining the neighbourhood (default one substitution).
selection_q (float) – Thymic-selection factor Q. The default is the paper’s average over V–J combinations; it is calibrated on OLGA’s TRB Pgen scale and there is no evidence it transfers to the other loci, so the benchmark titrates it per locus.
min_degree (int) – Only clonotypes with at least this many neighbours (self included) are tested; the rest are dropped. ALICE tests d(σ) > 2.
min_count (int) – Clonotypes below this count do not participate at all. Low-count variants of an abundant clonotype are usually sequencing error, and counting them as neighbours inflates every degree in their neighbourhood.
threads (int) – Worker threads (0 = all cores).

Returns:

the key columns, duplicate_count, n_neighbors (d(σ), self included), n_group (n, the V–J class’s nucleotide-clonotype count), pgen_ball (Σ Pgen over the closed ball), E (λ_σ), p_enrichment (Poisson tail), q_value (BH FDR over the clonotypes tested in this call) and locus — sorted by ascending p_enrichment. No threshold is applied: ALICE’s own papers use BH < 0.001 while the TCRnet framework paper used 0.05 for the same family of tests, and that was never reconciled — so the caller chooses.

Return type:

One row per tested clonotype

Raises:

ValueError – If no locus can be resolved, or the sample has no usable clonotypes.
KeyError – If a V/J call is neither a model allele nor a gene of one.

`vdjtools.overlap.cluster`#

Pairwise-distance matrices, clustering and MDS.

Pairwise sample distances and low-dimensional clustering.

Formalises the legacy CalcPairwiseDistances + ClusterSamples workflow (and the ad-hoc MDS the aging example notebook did inline): compute an all-pairs distance matrix from a repertoire-overlap metric, then embed it in 2-D (MDS) or build a hierarchy (hclust).

An overlap similarity is turned into a distance with the legacy per-metric normalisation (OverlapMetricNormalization):

F, F2, D (frequency/diversity overlaps, (0, 1]) → -log10(x + 1e-9);
R (correlation, [-1, 1]) → (1 - x) / 2;
jaccard (similarity index, [0, 1]) → 1 - x.

The diagonal is forced to 0 and the matrix is symmetric. When a scope is passed the (fuzzy) vdjtools.overlap.fuzzy.fuzzy_overlap_metrics() fuzzy_F is used instead of the exact metric.

vdjtools.overlap.cluster.DEFAULT_KEY = ('junction_aa', 'v_call', 'j_call')#: Default clonotype match key (CDR3 aa + V + J), matching the exact-overlap default.

vdjtools.overlap.cluster.pairwise_distances(samples, metric='F', key=('junction_aa', 'v_call', 'j_call'), scope=None, form='matrix')[source]#

All-pairs distance matrix over a collection of samples.

Parameters:

samples – A list of clonotype frames (named "0".."N-1") or a dict mapping sample name to frame.
metric (str) – Overlap similarity to base the distance on: "F", "F2", "D" (→ -log10), "R" (→ (1-x)/2), "jaccard" (→ 1-x), or the sequence-similarity-weighted "similarity_cosine" / "similarity_morisita" (TINA / Leinster-Cobbold, → 1-x). See the module docstring.
key – Exact-match clonotype key (default CDR3 aa + V + J); ignored when scope is given.
scope (str | None) – If set, use fuzzy overlap within this vdjmatch edit scope ("subs,ins,dels,total") and the fuzzy_F similarity instead of the exact metric (only metric="F" is valid then).
form (str) – "matrix" for a wide frame (a sample column plus one column per sample) or "long" for a sample_a, sample_b, distance frame.

Returns:

A symmetric distance matrix with a zero diagonal, in the requested form.

Return type:

DataFrame

vdjtools.overlap.cluster.cluster_samples(dist, method='mds', n_components=2, metadata=None)[source]#

Embed / cluster samples from a precomputed distance matrix.

Parameters:

dist (DataFrame) – A matrix-form distance frame from pairwise_distances() (a sample column plus one column per sample).
method (str) – "mds" — metric MDS (sklearn.manifold.MDS with dissimilarity="precomputed") → n_components coordinate columns mds1..mdsK; or "hclust" — average-linkage hierarchy (scipy.cluster.hierarchy) → a dendrogram leaf_order and a flat cluster label (fcluster into n_components clusters).
n_components (int) – MDS output dimensionality (method="mds") or the number of flat clusters (method="hclust").
metadata (DataFrame | None) – Optional frame carrying a sample column plus per-sample columns (e.g. age, group) to left-join onto the result for colouring.

Returns:

A pl.DataFrame with one row per sample and the embedding / cluster columns, plus any joined metadata.

Raises:

ImportError – If method="mds" and scikit-learn is not installed.
ValueError – If method is not "mds" or "hclust".

Return type:

DataFrame

`vdjtools.overlap.track`#

Clonotype tracking across a sample series.

Track clonotype frequencies across an ordered set of samples.

Reimplements the legacy TrackClonotypes: given several samples (e.g. a time course or an age series) it builds one row per clonotype and one frequency column per sample, so a clonotype’s trajectory can be read across the columns. Clonotypes present in at least one sample are kept, sorted by their summed frequency (most persistent/abundant first), optionally truncated to the top N.

vdjtools.overlap.track.DEFAULT_KEY = ('junction_aa', 'v_call', 'j_call')#: Default clonotype match key (CDR3 aa + V + J).

vdjtools.overlap.track.track_clonotypes(samples, order=None, top=None, key=('junction_aa', 'v_call', 'j_call'))[source]#

Pivot per-sample clonotype frequency into one column per sample.

Parameters:

samples – A dict mapping sample name to clonotype frame, or a list of frames (named "0".."N-1").
order – Sample names giving the left-to-right column order. Defaults to the samples’ natural order (dict insertion / list order). Names not present in samples are ignored.
top (int | None) – If given, keep only the top clonotypes by summed frequency.
key – Clonotype match key (default CDR3 aa + V + J).

Returns:

A pl.DataFrame with the key columns, one freq_<sample> column per sample in order (0.0 where the clonotype is absent), and a freq_sum column, sorted by freq_sum descending. Clonotypes present in at least one sample are included.

Return type:

DataFrame

Preprocessing (`vdjtools.preprocess`)#

Sample preprocessing and operations: downsampling, frequency error-correction, decontamination, segment/frequency/functional filters, VJ-usage batch-effect correction, and pooling / joining.

`vdjtools.preprocess.downsample`#

Read / unique / frequency downsampling (numpy multinomial).

Random down-sampling of clonotype frames (pure polars + numpy).

Reimplements the legacy vdjtools DownSampler / TopSampler family. Two resampling regimes, matching the legacy --unweighted switch:

reads (legacy default, weighted): draw size reads without replacement from the multiset of reads implied by duplicate_count. The legacy DownSampler shuffles a flattened per-read array and keeps the first size entries; the exact equivalent is the multivariate hypergeometric distribution (numpy.random.Generator.multivariate_hypergeometric). The task brief phrased this as “multinomial”, but sampling a sequencing library to a fixed depth is a without-replacement operation — a multinomial (with replacement) could return more reads of a clonotype than were observed — so the hypergeometric is used to stay faithful to the legacy behaviour and to the biology.
clones (legacy --unweighted): draw size unique clonotypes uniformly at random without replacement, keeping each one’s original count. Note the legacy clonotype-level mode is uniform, not count-weighted (weighting by count is what the read-level mode does).

vdjtools.preprocess.downsample.downsample(df, size, by='reads', seed=0)[source]#

Randomly down-sample a clonotype frame to a target size.

Parameters:

df (DataFrame) – A clonotype frame with a duplicate_count column.
size (int) – Target size — number of reads (by="reads") or number of unique clonotypes (by="clones").
by (str) – "reads" (default) draws size reads without replacement, weighted by duplicate_count (multivariate hypergeometric); "clones" draws size unique clonotypes uniformly without replacement, keeping their original counts.
seed (int) – Seed for the numpy random generator (reproducible output).

Returns:

A new clonotype frame with frequency recomputed. Clonotypes that drew zero reads (by="reads") are dropped. If size is greater than or equal to the available size the input is returned unchanged (legacy guard).

Raises:

ValueError – If by is not "reads" or "clones", or size < 0.

Return type:

DataFrame

vdjtools.preprocess.downsample.select_top(df, n, renormalize=True)[source]#

Select the top n clonotypes by duplicate_count.

Reimplements the legacy SelectTop / TopSampler (take the n largest clonotypes). Ties are broken by the frame’s existing order (a stable sort).

Parameters:

df (DataFrame) – A clonotype frame with a duplicate_count column.
n (int) – Number of top clonotypes to keep. If n is greater than or equal to the number of clonotypes, all are returned.
renormalize (bool) – If True (legacy default), recompute frequency within the selected subset so it sums to 1; if False, preserve the input frequencies (legacy --save-freqs).

Returns:

The top-n clonotype frame, sorted by descending duplicate_count.

Return type:

DataFrame

`vdjtools.preprocess.filter`#

Segment, frequency and functional (coding) filters.

Clonotype filtering (pure polars).

Reimplements the legacy vdjtools clonotype-filter family:

filter_productive() — AIRR-productive rearrangements (supersedes filter_functional, the legacy FunctionalClonotypeFilter / isCoding).
filter_frequency() — FilterByFrequency (FrequencyFilter + QuantileFilter).
filter_segment() — FilterBySegment (VFilter / DFilter / JFilter).
filter_by_sample() — ApplySampleAsFilter (IntersectionClonotypeFilter).

Every filter recomputes frequency within the surviving subset by default, and filter_productive() exposes that as recompute_frequencies=False for a caller who wants the file’s own frequencies left alone.

On the word “functional”. It is IMGT’s, and IMGT applies it to a germline gene (F / ORF / P), not to a rearrangement. What this module filters is AIRR productivity — a property of the rearranged sequence: no stop codon, junction in frame. The two are orthogonal; a perfectly productive rearrangement can use a pseudogene V. filter_productive() is therefore the name, and filter_functional() is kept as a deprecated alias.

vdjtools.preprocess.filter.productive_mask(df)[source]#

The productivity predicate for this frame, and the evidence it rests on.

Returns (expr, source) where source is the AIRR column(s) used, or "junction_aa" when none were present and productivity had to be derived from the amino-acid string.

The fallback reads a stop codon as * and an out-of-frame junction as one of the legacy markers. It is a proxy: it cannot see a defect in a splicing site or a regulatory element, which AIRR’s productive can.

Parameters:: df (DataFrame)
Return type:: tuple[Expr, str]

vdjtools.preprocess.filter.filter_productive(df, keep='productive', *, recompute_frequencies=True)[source]#

Keep only AIRR-productive rearrangements (or only the complement).

A rearrangement is productive when it can encode a receptor chain. Where the frame carries the AIRR annotation columns (productive, or stop_codon / vj_in_frame) those are authoritative; otherwise productivity is derived from junction_aa, where a stop codon is * and an out-of-frame junction carries one of the legacy markers [atgc#~_?]. A null junction_aa is treated as non-productive.

This is not IMGT functionality. That is a property of the germline gene (F / ORF / P) and is orthogonal to this one — see filter_functional_genes().

Parameters:

df (DataFrame) – A clonotype frame.
keep (str) – "productive" (default) or "nonproductive" for the complement.
recompute_frequencies (bool) – Renormalise frequency over the survivors. Default ``True``, which is the legacy behaviour and what almost every caller wants. Pass False to leave the file’s own frequencies untouched — useful when the frequencies are the quantity of interest and must stay comparable to the unfiltered file.

Returns:

The filtered frame.

Raises:

ValueError – If keep is not "productive" or "nonproductive".

Return type:

DataFrame

vdjtools.preprocess.filter.filter_functional(df, keep='coding')[source]#

Deprecated alias for filter_productive(). Use that instead.

“Functional” is IMGT’s word for a germline gene; this function filters rearrangements, which AIRR calls productive. Kept for one release so existing callers do not break silently.

Parameters:

df (DataFrame)
keep (str)

Return type:

DataFrame

vdjtools.preprocess.filter.filter_functional_genes(df, *, segments=('V', 'J'), keep=('F',), organism='human', locus=None, recompute_frequencies=True)[source]#

Keep rearrangements whose germline gene calls are IMGT-functional.

This is the other axis, and the one that actually deserves the word functional. IMGT classifies a germline gene as F (functional), ORF (an open reading frame, but a defect in splicing, regulatory elements or conserved-residue hydropathy — not functional), or P (pseudogene: a defect in the ORF itself). See https://www.imgt.org/IMGTindex/functionality.php

It is orthogonal to filter_productive(). A rearrangement can be perfectly in frame with no stop codon — AIRR-productive — while using a pseudogene V; and a functional V gene can rearrange out of frame. Filtering one says nothing about the other.

A call this function cannot resolve against the germline reference is kept, not dropped: an unrecognised gene name means our reference is incomplete or the caller uses a different nomenclature, and silently discarding those rows would be a vocabulary bug reported as biology.

Parameters:

df (DataFrame) – A clonotype frame.
segments (tuple[str, ...]) – Which calls to check — any of "V", "D", "J".
keep (tuple[str, ...]) – IMGT functionality codes to keep. ("F",) is strict; ("F", "ORF") is the common looser choice.
organism (str) – Passed to the germline reference.
locus (str | None) – Locus to load the reference for. Inferred from the V calls when omitted.
recompute_frequencies (bool) – Renormalise frequency over the survivors. Default True.

Returns:

The filtered frame.

Return type:

DataFrame

vdjtools.preprocess.filter.MIN_JUNCTION_AA = 5#

Default junction_aa length bounds, INCLUSIVE, in amino acids.

A CDR3 shorter than 5 aa cannot span the Cys104..Phe118 anchors with any diversity between them, and a junction longer than 60 aa is beyond anything the germline can produce – both are almost always a misparse or a chimeric read rather than a receptor. Deliberately wide: this is a sanity bound, not a biological filter, and real junctions sit far inside it.

vdjtools.preprocess.filter.MAX_JUNCTION_AA = 60#

Default junction_aa length bounds, INCLUSIVE, in amino acids.

vdjtools.preprocess.filter.filter_length(df, *, min_len=5, max_len=60, keep='within', recompute_frequencies=True)[source]#

Keep clonotypes whose junction_aa length is within bounds, inclusive.

Both bounds are inclusive: min_len=5 keeps a 5-mer, max_len=60 keeps a 60-mer.

This is a data-sanity filter, not a receptor-biology one. It catches misparses, truncated reads and chimeras; it does not encode a claim about what lengths are immunologically interesting. Note that neither this nor filter_productive() filters on length by default – nothing upstream in this package has ever imposed a length bound, so switching this on will change counts on any corpus that carries junk.

CDR3 vs junction. These bounds are on junction_aa, which includes the Cys104 and Phe118 anchors and is therefore two residues longer than the IMGT CDR3. Subtract 2 if you are reasoning in CDR3 lengths.

A null junction_aa is dropped by keep="within" – an absent junction has no length.

Parameters:

df (DataFrame) – A clonotype frame.
min_len (int) – Shortest junction_aa to keep, inclusive.
max_len (int) – Longest junction_aa to keep, inclusive.
keep (str) – "within" (default) or "outside" for the complement – useful for inspecting what a bound would discard before committing to it.
recompute_frequencies (bool) – Renormalise frequency over the survivors. Default True.

Returns:

The filtered frame.

Raises:

ValueError – If keep is unknown, or min_len exceeds max_len.

Return type:

DataFrame

vdjtools.preprocess.filter.filter_frequency(df, min_freq=None, top_quantile=None)[source]#

Keep abundant clonotypes by frequency threshold and/or top quantile.

Reimplements FilterByFrequency (a composite of FrequencyFilter and QuantileFilter). Both criteria, when given, are combined with AND:

min_freq: keep clonotypes with frequency >= min_freq.
top_quantile: keep the top clonotypes (by duplicate_count) whose cumulative original frequency, including the clonotype itself, is at most top_quantile of the full-sample total frequency. This matches the legacy QuantileFilter: it walks the count-sorted sample accumulating frequency and drops the first clonotype that would push the running fraction above the threshold (so top_quantile=0.25 keeps roughly the top 25% of the read mass). The denominator is the full-sample frequency total (~1.0), and only clonotypes that already passed min_freq contribute to the cumulative (legacy filters short-circuit in the order count/freq/quantile).

Parameters:

df (DataFrame) – A clonotype frame with duplicate_count and frequency columns.
min_freq (float | None) – Minimum per-clonotype frequency (e.g. legacy default 0.01). None disables it.
top_quantile (float | None) – Top read-mass quantile to retain (e.g. legacy default 0.25). None disables it.

Returns:

The filtered frame, sorted by descending duplicate_count, with frequency recomputed.

Return type:

DataFrame

vdjtools.preprocess.filter.filter_segment(df, v=None, d=None, j=None, keep=True)[source]#

Keep or remove clonotypes by V/D/J segment membership.

Reimplements FilterBySegment. A clonotype matches when its V segment is in v and its D segment in d and its J segment in j (only the lists that are supplied constrain; unsupplied loci always pass). Matching is a prefix match, so incomplete names act as wildcards and are allele-insensitive (TRBV12 matches TRBV12-3*01).

Parameters:

df (DataFrame) – A clonotype frame with v_call / d_call / j_call columns.
v (list[str] | None) – V-segment query names (prefixes). None leaves V unconstrained.
d (list[str] | None) – D-segment query names. None leaves D unconstrained.
j (list[str] | None) – J-segment query names. None leaves J unconstrained.
keep (bool) – If True (default) keep matching clonotypes; if False remove them (legacy --negative).

Returns:

The filtered frame with frequency recomputed.

Return type:

DataFrame

vdjtools.preprocess.filter.filter_by_sample(df, other, keep=True, key=('junction_aa', 'v_call', 'j_call'))[source]#

Keep or remove clonotypes by exact-key presence in another sample.

Reimplements ApplySampleAsFilter / IntersectionClonotypeFilter: build the key set of other and keep (or, with keep=False, remove) the clonotypes of df whose key is present in it. Matching is an exact match on the key columns.

Parameters:

df (DataFrame) – The clonotype frame to filter.
other (DataFrame) – The filter sample; only its key columns are used.
keep (bool) – If True (default) keep clonotypes present in other; if False remove them (legacy --negative).
key (tuple[str, ...]) – Columns forming the match key (default ("junction_aa", "v_call", "j_call") — legacy “strict”-style at the aa level).

Returns:

The filtered frame with frequency recomputed.

Return type:

DataFrame

`vdjtools.preprocess.correct`#

Frequency-based sequencing-error correction.

Frequency-based sequencing-error correction (polars; neighbour search via seqtree).

Reimplements the legacy vdjtools Corrector. Clonotypes whose junction_nt sit within a few substitutions of a much more abundant clonotype are treated as PCR / sequencing errors of that parent and merged into it.

The <= max_mismatches substitution neighbour search is delegated to seqtree (the “delegate search to seqtree” convention), which builds an immutable fuzzy index and returns, per query, every within-budget neighbour and its substitution count. Only the abundance-ratio decision is done here.

vdjtools.preprocess.correct.correct(df, max_mismatches=2, ratio=0.05, same_vj=False)[source]#

Merge low-frequency sequencing-error clonotypes into their parents.

Reimplements Corrector. Clonotype pairs whose junction_nt are within max_mismatches substitutions are compared by abundance: a smaller clonotype is merged into a larger one when its count is below ratio ** m times the larger’s (m = number of substitutions). All decisions use the original counts (a single, order-independent pass, as in the legacy parallel stream).

Parameters:

df (DataFrame) – A clonotype frame with junction_nt and duplicate_count columns.
max_mismatches (int) – Maximum substitutions for two clonotypes to be neighbours (legacy default 2); insertions/deletions are not considered.
ratio (float) – Per-mismatch parent/child count ratio (legacy default 0.05). A child is merged when child_count < ratio ** m * parent_count.
same_vj (bool) – If True (the opt-in “match-segment” mode, legacy -a / --match-segment) only clonotypes sharing the exact v_call and j_call are compared; if False (default, matching legacy fidelity — legacy Corrector is segment-agnostic by default) all clonotypes are compared regardless of segment.

Returns:

The corrected frame (errors dropped, parents’ counts increased), sorted by descending duplicate_count with frequency recomputed. Rows with a null junction_nt pass through uncorrected.

Raises:

ImportError – If seqtree is not importable (it is a base dependency).

Return type:

DataFrame

`vdjtools.preprocess.decontaminate`#

Cross-sample contamination removal.

Cross-sample contamination removal (pure polars).

Reimplements the legacy Decontaminate / RatioFilter: drop clonotypes that occur in another sample at a much higher abundance, on the assumption that such a clonotype leaked (index hopping / carry-over) into the current sample.

vdjtools.preprocess.decontaminate.decontaminate(df, others, ratio=20.0, by='freq', key=('junction_nt', 'v_call', 'j_call'))[source]#

Remove clonotypes dominated by an exact match in another sample.

Reimplements RatioFilter. A clonotype of df is removed when some sample in others carries the same clonotype (exact match on key) at an abundance >= ratio times its abundance here (legacy keeps a clonotype iff max_other < here * ratio, i.e. removes on >=). Abundance is the within-sample frequency (by="freq") or read duplicate_count (by="reads").

Note

The legacy Decontaminate --read-based branch was a no-op — both CLI branches constructed the same frequency-based RatioFilter (a documented TODO). by="reads" here implements the intended read-count comparison.

Parameters:

df (DataFrame) – The clonotype frame to decontaminate.
others (list[DataFrame]) – The other samples that may be contamination sources.
ratio (float) – Parent-to-child abundance ratio (legacy default 20).
by (str) – "freq" (default) compares within-sample frequencies; "reads" compares raw read counts.
key (tuple[str, ...]) – Exact match key (legacy default is the strict key ("junction_nt", "v_call", "j_call")).

Returns:

The decontaminated frame with frequency recomputed.

Raises:

ValueError – If by is not "freq" or "reads".

Return type:

DataFrame

`vdjtools.preprocess.pool`#

Pool samples into a joint clonotype table.

Pool clonotypes across samples (pure polars).

Reimplements the legacy PoolSamples / SampleAggregator / PooledSample: collapse clonotypes across a set of samples by a chosen match key, summing counts and recomputing frequency, and annotate each pooled clonotype with its incidence, occurrence count and convergence.

vdjtools.preprocess.pool.KEY_SETS: dict[str, tuple[str, ...]] = {'aa': ('junction_aa',), 'aaV': ('junction_aa', 'v_call'), 'aaVJ': ('junction_aa', 'v_call', 'j_call'), 'nt': ('junction_nt',), 'ntV': ('junction_nt', 'v_call'), 'ntVJ': ('junction_nt', 'v_call', 'j_call'), 'strict': ('junction_nt', 'v_call', 'j_call')}#: Legacy OverlapType match keys -> the clonotype columns they compare on.

vdjtools.preprocess.pool.resolve_key(key)[source]#

Resolve a match-key name (or explicit column tuple) to key columns.

Parameters:: key (str | tuple[str, ...]) – One of the legacy overlap-type names ("strict", "nt", "ntV", "ntVJ", "aa", "aaV", "aaVJ") or an explicit tuple/list of column names.
Returns:: The list of clonotype columns forming the key.
Raises:: ValueError – If key is an unknown name.
Return type:: list[str]

vdjtools.preprocess.pool.pool_samples(samples, key='aa', sample_col=None)[source]#

Pool clonotypes across samples, summing counts and recomputing frequency.

Reimplements PooledSample (legacy default match key "aa"). For each distinct match key across all samples the pooled clonotype carries:

duplicate_count — summed read count across all samples.
frequency — pooled count over the pool’s total reads.
incidence — number of distinct samples the clonotype occurs in.
occurrences — total number of clonotype rows aggregated (a clonotype can appear once per sample, and, for amino-acid keys, via several nucleotide variants within one sample).
convergence — number of distinct nucleotide variants (junction_nt + v_call + j_call, the legacy strict key) collapsed into the pooled clonotype. This is 1 for nucleotide-level keys and counts convergent recombination for amino-acid-level keys.

The representative non-key fields (e.g. the junction_nt of an amino-acid pool) are taken from the most abundant contributing row (legacy MaxClonotypeAggregator).

Parameters:

samples (list[DataFrame]) – A list of per-sample clonotype frames, or a single long frame (pass [df]) split by sample_col.
key (str | tuple[str, ...]) – Match key name or explicit column tuple (see resolve_key()).
sample_col (str | None) – If the input is a single long frame carrying a sample-id column, its name; used only to count incidence correctly.

Returns:

A pooled clonotype frame in the canonical schema plus incidence, occurrences and convergence, sorted by descending duplicate_count.

Return type:

DataFrame

`vdjtools.preprocess.join`#

Join samples on shared clonotypes.

Join clonotypes across samples into a joint table (pure polars).

Reimplements the legacy JoinSamples / JointSample / JointClonotype: build the table of clonotypes present in at least min_samples samples, keep each member’s per-sample frequency, and summarise each joint clonotype by the geometric mean of its member frequencies.

vdjtools.preprocess.join.JITTER = 1e-09#: Legacy MathUtil.JITTER — the epsilon added to each per-sample frequency so that a geometric mean over samples where the clonotype is absent stays finite.

vdjtools.preprocess.join.join_samples(samples, key='aa', min_samples=2, names=None)[source]#

Join clonotypes across samples, keyed by an overlap-type match key.

Reimplements JointSample (legacy default key "aa", times-detected 2). A joint clonotype is kept when it is present in at least min_samples samples (legacy OccurrenceJoinFilter). Its joint frequency is the geometric mean of the per-sample frequencies (absent samples contribute JITTER),

base = (∏_i (freq_i + JITTER))**(1/n_samples)

then normalised so passing joint frequencies sum to 1 (legacy calcFreq = base / Σ base). The joint count is floor(base / min base) so that the smallest joint clonotype has count 1 (legacy calcCount = base / min base).

Parameters:

samples (list[DataFrame]) – A list of per-sample clonotype frames.
key (str | tuple[str, ...]) – Match key name or explicit column tuple (see vdjtools.preprocess.pool.resolve_key()).
min_samples (int) – Minimum number of samples a clonotype must occur in to be kept.
names (list[str] | None) – Optional per-sample names for the freq_* / count_* columns; defaults to the sample indices 0..n-1.

Returns:

the key columns, one freq_<name> and count_<name> column per sample, incidence (number of samples present), frequency (normalised geometric-mean joint frequency) and duplicate_count (normalised joint count, smallest = 1), sorted by descending frequency.

Return type:

A joint clonotype frame

Raises:

ValueError – If names is given but its length differs from samples.

`vdjtools.preprocess.batch`#

VJ-usage batch-effect correction.

VJ-usage batch-effect correction + clonotype-table application (pure polars/numpy).

Different sequencing batches carry systematic V/J gene-usage biases (primer mixes, amplification, extraction). This module removes that batch-specific offset so that per-sample VJ usage becomes comparable across batches, and can then push the corrected usage back onto a sample’s clonotype table (rescale + resample).

Two stages#

correct_vj_usage() — batch-correct per-sample V-J gene usage. Two transforms:
- transform="location" (default) — the classic location adjustment (the location term of ComBat; Johnson, Li & Rabinovic, Biostatistics 2007) on gene-usage log-probabilities. Per (locus, gene, batch) the batch mean mu_batch of log p is replaced by the grand mean mu_grand: log_corrected = log_p - mu_batch + mu_grand, then exp and renormalise. Location only — no scale term.
- transform="sigmoid" — the σ-standardised, grand-mean-preserving correction of Vlasova, Nekrasova, Komkov, … Britanova, Shugay, Genome Medicine 2026;18:20 (DOI 10.1186/s13073-025-01589-4). Per (locus, gene, batch) compute a z-score Z = (log p - mu_batch) / sigma_batch (capped at ±z_cap), then map it back to a probability with a sigmoid that preserves the pooled grand-mean usage P_avg(gene): P_final = 2 * P_avg / (1 + exp(-Z)) (Z=0 → P_avg), renormalised per (sample, locus). The paper uses the plain mean/σ of the log-normal (Shapiro–Wilk-validated), i.e. winsor_q=None (default). Winsorization (winsor_q=0.025) is an optional robustness knob for the noisy usage-as-features regime (many shallow / RNA-seq repertoires), not for deep-repertoire correction; legacy mirpy v2 mir.basic.gene_usage.compute_batch_corrected_gene_usage winsorized by default (and used p*exp(Z) for its own map — the 2*P_avg*sigmoid(Z) map here is the paper’s Methods formula).
apply_vj_correction() — apply a sample’s corrected usage back to its clonotype table: reweight each clonotype by P_final(G) / P(G) for its gene, then (by default) roulette-wheel resample to a new integer-count table with the corrected usage. Port of legacy mirpy v2 mir.common.sampling.resample_to_gene_usage.

Alleles are stripped for the usage key (mirpy convention); clonotypes keep their original allele calls.

vdjtools.preprocess.batch.correct_vj_usage(samples_or_df, batch_col, sample_col='sample_id', weighted=True, pseudocount=1.0, transform='location', z_cap=6.0, winsor_q=None)[source]#

Batch-correct per-sample V-J gene usage.

Parameters:

samples_or_df (DataFrame | list[DataFrame]) – A single long clonotype frame carrying sample_col and batch_col columns, or a list of such frames (concatenated).
batch_col (str) – Column naming each sample’s batch.
sample_col (str) – Column naming the sample (default "sample_id").
weighted (bool) – If True (default) usage counts reads (duplicate_count); if False it counts clonotypes.
pseudocount (float) – Laplace smoothing constant added per gene (default 1.0).
transform (str) – "location" (default) for the ComBat location adjustment (log_p - mu_batch + mu_grand), or "sigmoid" for the σ-standardised, grand-mean-preserving z-score map P_final = 2*P_avg/(1+exp(-Z)) (Vlasova et al. 2026). See the module docstring.
z_cap (float) – For transform="sigmoid", clip the z-score to ±z_cap (default 6.0); ignored for "location".
winsor_q (float | None) – If set (e.g. 0.025), winsorize each (gene, batch) group’s log p to [winsor_q, 1-winsor_q] before taking the batch mean/σ. The default None uses the plain mean/σ of the log-normal — the paper’s method (Vlasova et al. 2026, Shapiro–Wilk-validated). Winsorization is a robustness knob for the noisy usage-as-features regime (many shallow / RNA-seq repertoires used for UMAP etc.), not for deep-repertoire correction.

Returns:

A long frame with one row per (sample, locus, v_call, j_call) over the union of genes per locus, with columns: sample_id, batch, locus, v_call, j_call, count, p (smoothed raw usage probability) and p_corrected (batch-corrected, renormalised usage probability), sorted by sample_id, locus, v_call, j_call. Feed this straight into apply_vj_correction().

Raises:

ValueError – If sample_col / batch_col is missing or transform is unknown.

Return type:

DataFrame

vdjtools.preprocess.batch.apply_vj_correction(sample_df, corrected_usage, *, scope='vj', weighted=True, resample=True, sample_id=None, seed=0)[source]#

Apply batch-corrected V/J usage back to a sample’s clonotype table.

Each clonotype is reweighted by its gene’s correction factor P_corrected(G) / P(G) (G = the clonotype’s V/J gene at scope), then either resampled or rescaled to a new clonotype table whose V/J usage matches the corrected usage. Port of legacy mirpy v2 resample_to_gene_usage.

Parameters:

sample_df (DataFrame) – One sample’s canonical clonotype frame (v_call, j_call, duplicate_count). Original allele calls are preserved on output.
corrected_usage (DataFrame) – Output of correct_vj_usage(). If it covers more than one sample, pass sample_id (or slice it) to select this sample’s rows.
scope (str) – Gene scope for the correction key: "vj" (default, per V-J pair), "v", or "j" (the VJ usage is marginalised for "v"/"j").
weighted (bool) – If True (default) the resampling weight is duplicate_count * factor (roulette-wheel over reads); if False it is factor (over clonotypes).
resample (bool) – If True (default), multinomial roulette-wheel resample to a new integer-count table at the sample’s original total read count (Vlasova et al. 2026). If False, deterministically rescale to the expected counts.
sample_id (str | None) – The sample to select from corrected_usage when it holds several.
seed (int) – Seed for the multinomial resample (resample=True).

Returns:

A canonical clonotype frame (same columns as sample_df) with corrected duplicate_count and recomputed frequency; zero-count clonotypes dropped.

Raises:

ValueError – If corrected_usage covers multiple samples and none is selected, if scope is unknown, or if no clonotype has a positive corrected weight.

Return type:

DataFrame

Repertoire dynamics (`vdjtools.dynamics`)#

Longitudinal within-donor clonotype testing across timepoints. Three complementary tools: a per-clonotype paired test (Ayestaran 2024) with per-pair effective sample size, deterministic downscale, two-tailed Fisher, and a five-way classification (emergent / expanded / persistent / contracted / vanishing); the same test on 1-Hamming / 1-Levenshtein CDR3 metaclonotype groups; and the VDJtrack size-bucket recapture model with an edgeR NB-exact caller (Pavlova, Zvyagin & Shugay 2024). The sibling of vdjtools.biomarker — that tests incidence across subjects, this tests frequency across timepoints.

`vdjtools.dynamics.paired`#

Per-pair N_eff estimation and the paired clonotype test.

Paired within-donor clonotype test: per-pair N_eff, downscale, two-tailed Fisher.

The method is Ayestaran (PhD, Cambridge 2024), Ch. 2. Sequencing is a two-step sampling process — a true frequency f is sampled into N_S1 molecules, which are sampled into N_seq reads — so the sample size driving the noise is the harmonic sum (Eq. 2.5–2.6):

1 / N_eff = 1 / N_S1 + 1 / N_seq

and is dominated by the smaller step. Where N_S1 << N_seq (the usual case: ~1e6 PBMCs into ~2e7 reads) the library is heavily oversampled, the observed counts are far noisier than N_seq implies, and a Fisher test that assumes N_seq is the sample size calls enormous numbers of clones changed when nothing changed (thesis Fig. 2.13).

N_eff is a property of the pair, not of a sample or a cohort: it is read off the mean–variance relationship between two samples, so a global N_eff is not merely a bad idea, it is undefined. That is why _downscale() is private — a public “rescale to N” would invite exactly the cohort-wide normalisation this method exists to avoid.

Estimation fits the thesis’s Eq. 2.3, log σ² = log f + log(1/N), with the slope fixed at 1 — which is just a weighted mean of log σ² − log f, so it needs numpy, not a regression library. It differs from the thesis in which variance it fits, and that choice is load-bearing enough to record here.

The thesis (p. 19) bins clones by their frequency in sample B and takes the mean and variance of the A frequencies in each bin, then multiplies by 2: binning on the noisy f_B instead of the true f mixes a range of true frequencies into each bin, roughly doubling the apparent variance, so the fitted value is half the real N_eff.

We instead bin on the mean f̄ = (f_A + f_B)/2 and fit var(f_A − f_B) = 2f/N. The difference cancels the true frequency, so there is no mis-binning to correct and no correction factor to mis-set — arguably closer to Eq. 2.3’s own derivation, which assumed f known.

Why, measured on simulated two-step data against a planted N_eff (see tests/python/test_dynamics_paired.py): the ×2 is exact only when the clone-size distribution is flat, because only then is the bin contamination equal to the sampling variance. Real repertoires are heavy-tailed, where the contamination is smaller and a fixed ×2 overshoots by ~25%. That is not cosmetic — it inflates the false-positive rate ~1.27× and fails the thesis’s own p-value-uniformity acceptance test (Fig. 2.13) in 3 of 4 regimes, including precisely the oversampled ones (N_S1 << N_seq) the method exists to handle:

regime (N_S1/N_seq)   planted    thesis x2   gate      difference   gate
200k / 2M             181,818    1.25x       FAIL      1.07x        pass
 50k / 2M              48,780    1.22x       FAIL      1.02x        pass
200k / 1M             166,667    1.23x       FAIL      0.99x        pass
  1M / 300k           230,769    0.97x       pass      0.95x        pass

The test itself is not at fault: pinned to the planted N_eff it is perfectly calibrated (conservative, as an exact test on discrete counts should be). The bias was entirely the estimator’s.

The test itself is Eq. 2.4: conditioning on the total r1 + r2 cancels the unknown true frequency exactly, leaving a hypergeometric — i.e. a two-tailed Fisher exact test on [[c_a, c_b], [R_a - c_a, R_b - c_b]]. It uses the minimum-likelihood two-sided convention (R / scipy), because rounding leaves the two library sizes near-equal rather than equal, and on near-equal margins the doubling convention differs by up to 2x.

vdjtools.dynamics.paired.DEFAULT_KEY = ('junction_aa', 'v_call', 'j_call')#: Default clonotype match key (CDR3 aa + V + J), as in vdjtools.overlap.track.

vdjtools.dynamics.paired.CLASSES = ('emergent', 'expanded', 'persistent', 'contracted', 'vanishing', 'untested')#: The five dynamics classes, plus untested for clones below the testability floor.

vdjtools.dynamics.paired.estimate_neff(a, b, *, key=('junction_aa', 'v_call', 'j_call'), min_count=2, bins=25, min_bin=10, ref_n=None)[source]#

Estimate the pair’s effective sample size from its own mean–variance scaling.

Parameters:

a (DataFrame) – First sample of the pair (canonical clonotype frame).
b (DataFrame) – Second sample of the pair.
key – Clonotype match key (default CDR3 aa + V + J).
min_count (int) – Discreteness floor for the fit — clonotypes whose mean frequency across the pair is below this many counts are excluded from the mean–variance fit, because at 1–2 counts the frequency is too discrete to estimate a variance from. This is not the testability floor used by test_pair() (see its min_total); it is the most sensitive knob here, so the benchmark sweeps it rather than trusting a default.
bins (int) – Number of log-spaced frequency bins.
min_bin (int) – Minimum clonotypes in a bin for it to contribute to the fit.
ref_n (float | None) – Reference sample size for the outlier pre-filter. None (default) derives it from the data with an unfiltered first pass — the thesis hard-codes 200,000, which is a property of its cohort (2x one of our datasets’ whole libraries and 33x another’s), not a constant. Pass a number to pin it.

Returns:

The estimated N_eff.

Raises:

ValueError – If fewer than two frequency bins are usable — a pair too shallow, or with too few shared clonotypes, to fit. Never returns a silent fallback: a wrong N_eff silently mis-scales every downstream p-value.

Return type:

vdjtools.dynamics.paired.test_pair(a, b, *, key=('junction_aa', 'v_call', 'j_call'), neff='auto', min_total=6, alpha=0.01, **neff_kw)[source]#

Test every clonotype for a within-donor frequency change between two samples.

Downscales both samples to the pair’s N_eff and applies a two-tailed Fisher exact test to [[c_a, c_b], [R_a - c_a, R_b - c_b]]. Conditioning on c_a + c_b cancels the unknown true frequency exactly (thesis Eq. 2.4), so the test needs no estimate of it.

Parameters:

a (DataFrame) – The earlier sample (e.g. pre-vaccination).
b (DataFrame) – The later sample (e.g. post-vaccination). Direction is reported b-relative-to-a.
key – Clonotype match key.
neff (float | str) – "auto" estimates it from the pair via estimate_neff(); a float pins it; None skips the downscale entirely — correct only when the counts are already molecule/UMI counts rather than reads (thesis p. 86), since then there is no oversampling to undo.
min_total (int) – Testability floor — clonotypes with a combined downscaled count below this are classed untested rather than tested and called non-significant. Below ~6 the discrete hypergeometric cannot reach a small p at all, so testing them only costs multiple-testing burden. Distinct from estimate_neff’s min_count.
alpha (float) – FDR threshold (Benjamini–Hochberg q) for calling a change significant.
**neff_kw – Forwarded to estimate_neff() when neff="auto".

Returns:

One row per clonotype with the key columns, count_a/count_b (downscaled), f_a/f_b (original within-sample frequencies), p_value, q_value and dynamics — one of emergent (absent from a), expanded, persistent (no evidence of change), contracted, vanishing (absent from b), or untested. The classes partition the frame.

Raises:

ValueError – If neff="auto" and the pair cannot be fit (see estimate_neff()).

Return type:

DataFrame

`vdjtools.dynamics.groups`#

Metaclonotype (1-Hamming / 1-Levenshtein CDR3 ball) grouping before the paired test.

Metaclonotype-grouped dynamics — cluster near-variant clonotypes, then test the group.

A vaccine- or antigen-driven response is often convergent: many T-cell clones with slightly different CDR3s recognise the same epitope, so the expansion signal is spread thin across a family of near-variants and a per-exact-clonotype test (vdjtools.dynamics.test_pair()) sees each one only weakly. Collapsing a CDR3-neighbourhood into one metaclonotype first, then running the same paired test on the group, concentrates that signal.

The grouping is delegated to vdjtools.biomarker.metaclonotypes() (native seqtree / vdjmatch fuzzy search, single-linkage components) — scope="1,0,0,1" is a 1-substitution (Hamming) ball, "1,1,1,1" a 1-edit (Levenshtein) ball. This is the VDJtrack 1-mismatch annotation-expansion (example.Rmd / vaccination.Rmd) generalised into the test itself.

vdjtools.dynamics.groups.test_metaclonotypes(a, b, *, scope='1,0,0,1', match_v=True, match_j=True, threads=0, **test_kw)[source]#

Group CDR3-neighbour clonotypes into metaclonotypes, then test_pair() the groups.

The union of both samples’ clonotypes is clustered once into metaclonotypes; each sample’s clonotypes are relabelled by meta_id; and the paired within-donor test runs on the per-meta_id summed counts — so a convergent expansion counts as a single feature.

Parameters:

a (DataFrame) – The earlier sample (canonical clonotype frame).
b (DataFrame) – The later sample.
scope (str) – vdjmatch edit scope "subs,ins,dels,total" — "1,0,0,1" (default) is a 1-substitution (Hamming) ball, "1,1,1,1" a 1-edit (Levenshtein) ball, "0,0,0,0" reduces to exact grouping.
match_v (bool) – Require the same v_call to group two clonotypes.
match_j (bool) – Require the same j_call to group.
threads (int) – Worker threads for the native search (0 = all cores).
**test_kw – Forwarded to test_pair() (e.g. neff, min_total, alpha).

Returns:

test_pair()’s output keyed by meta_id, plus a representative junction_aa (and v_call/j_call when matched) and n_variants per group, sorted by q_value.

Raises:

ImportError – If vdjmatch is not importable (it is a base dependency).

Return type:

DataFrame

`vdjtools.dynamics.capture`#

The VDJtrack recapture model: size buckets, Poisson capture, Beta credible intervals, log-linear group test.

Clonotype recapture model — the VDJtrack size-bucket capture test.

Where vdjtools.dynamics.test_pair() calls each clonotype expanded/contracted, the capture model asks a group question: does an annotated set of clonotypes (e.g. antigen- specific, or “emerging”) persist / recapture across a time course more than the rest, after controlling for what drives recapture on its own — clonotype size and repertoire diversity?

Under Poisson sampling a clonotype of population frequency f is recaptured in a sample of depth R with probability P = 1 − exp(−f·R) (poisson_capture()), so recapture rate rises monotonically with pre-sample size — singleton < doubleton < tripleton < large. The model bins clonotypes into those four size classes, measures the recapture fraction per (donor, group, size class), puts a Beta(captured, missing) posterior on it for uncertainty, and tests the group effect with a log-linear model log(recapture) ~ size + group + log(div_ratio) (and a per-bucket paired test across donors).

Port of the group’s VDJtrack R pipeline (github.com/antigenomics/vdjtrack); the method is Pavlova, Zvyagin & Shugay, Front Immunol 2024, DOI 10.3389/fimmu.2024.1321603. The two-step sampling correction and the per-clonotype test live in vdjtools.dynamics.paired.

vdjtools.dynamics.capture.SIZE_CLASSES = ('singleton', 'doubleton', 'tripleton', 'large')#: Pre-sample size classes, ordered (singleton=1, doubleton=2, tripleton=3, large=4+).

vdjtools.dynamics.capture.size_class(count='duplicate_count')[source]#

Polars expression classifying a clonotype count into a SIZE_CLASSES bucket.

Parameters:: count (str | Expr)
Return type:: Expr

vdjtools.dynamics.capture.poisson_capture(freq, depth)[source]#

Poisson recapture probability P = 1 − exp(−f·R).

The probability that a clonotype of population frequency freq yields at least one molecule in a sample of depth (total reads/UMIs) — the generative model the recapture curve is read against (Pavlova 2024; aging_capture_model.Rmd).

Parameters:

freq – Clonotype population frequency (scalar or array).
depth – Sampling depth R (scalar or array).

Returns:

The capture probability, same shape as the broadcast of freq and depth.

Return type:

ndarray

vdjtools.dynamics.capture.capture_rates(pre, post, *, key=('junction_aa', 'v_call', 'j_call'), group_col=None, donor=None)[source]#

Recapture fraction per (size class[, group]) between one pre and post sample.

Each pre clonotype is binned by its pre-sample size (SIZE_CLASSES) and marked captured if its key reappears in post. Rows are then counted per size class (and group_col level, if given): n_captured = α, n_total = α+β, and the recapture rate gets a Beta(α+1, β+1) posterior (Laplace) with a 95% credible interval.

Parameters:

pre (DataFrame) – The earlier sample (canonical clonotype frame).
post (DataFrame) – The later sample.
key – Clonotype match key (default CDR3 aa + V + J).
group_col (str | None) – A column of pre carrying a group label (e.g. antigen specificity). If None, all clonotypes are one group "all".
donor (str | None) – If given, added as a constant donor column (for capture_test() over a cohort assembled by concatenating per-donor calls).

Returns:

One row per (donor,) group, size_class with n_captured, n_total, capture_rate (posterior mean), ci_lo, ci_hi.

Return type:

DataFrame

vdjtools.dynamics.capture.capture_test(rates, *, group_col='group', div_ratio=None)[source]#

Test the group effect on recapture with a log-linear model.

Fits log(capture_rate) ~ size_class + group [+ log(div_ratio)] by OLS (the VDJtrack lm; example.Rmd), so the group coefficient’s p-value asks whether the annotated group recaptures differently than baseline after accounting for clonotype size and (if supplied) the per-donor diversity ratio. Size classes and group are dummy-coded (first level = reference); a numeric log(div_ratio) is added per donor when div_ratio is given.

Parameters:

rates (DataFrame) – Stacked capture_rates() output (one row per donor × group × size class).
group_col (str) – The group column (default "group"); coded against its first level.
div_ratio (dict | None) – Optional {donor: div_after/div_before} — adds log(div_ratio) as a covariate (needs a donor column in rates).

Returns:

A coefficient table (term, estimate, std_error, p_value); the group row is the effect of interest.

Return type:

DataFrame

vdjtools.dynamics.capture.capture_paired_test(rates, *, group_col='group', donor_col='donor')[source]#

Per-size-class paired t-test of log recapture, group vs baseline across donors.

The VDJtrack per-bucket test (vaccination.Rmd): within each size class, pair donors and compare the two group levels’ log(capture_rate) with a paired t-test — robust to the donor-to-donor baseline that the log-linear model absorbs into an intercept. Requires exactly two group levels.

Parameters:

rates (DataFrame) – Stacked capture_rates() output.
group_col (str) – The two-level group column.
donor_col (str) – The donor column to pair on.

Returns:

size_class, n_pairs, t, p_value (null where fewer than two complete donor pairs exist).

Return type:

One row per size class

`vdjtools.dynamics.expansion`#

The edgeR-style negative-binomial exact expansion caller (TMM + qCML + beta-binomial exact test).

edgeR-style negative-binomial exact test for clonotype expansion between two timepoints.

The complementary per-clone caller of Pavlova, Zvyagin & Shugay 2024 (§2.5): the classic edgeR pipeline — TMM library normalization, a common dispersion by quantile-adjusted conditional maximum likelihood (qCML), and the negative-binomial exact test — used to call vaccine-associated clonotypes (|log2 FC| >= 5 & p <= 0.01 in the paper).

For two libraries with one replicate each, the NB exact test collapses to a closed form: after equalizing library sizes, each count is NB(mean=λ, size=r=1/φ) and the conditional law of y_a given the total t = y_a + y_b is Beta-Binomial(t, r, r) — λ cancels — so both the qCML dispersion fit and the two-sided (“minimum-likelihood”) exact p-value are beta-binomial computations (scipy.stats.betabinom). This is the same reduction edgeR’s exactTest uses.

The per-clonotype vdjtools.dynamics.test_pair() (Ayestaran N_eff + Fisher) is the default caller; this one models over-dispersion explicitly instead of correcting the sample size.

vdjtools.dynamics.expansion.expansion_test(a, b, *, key=('junction_aa', 'v_call', 'j_call'), dispersion='auto', min_total=6, log2fc=1.0, alpha=0.05)[source]#

Call clonotypes expanded/contracted between two samples by the edgeR NB exact test.

Parameters:

a (DataFrame) – The earlier sample (canonical clonotype frame).
b (DataFrame) – The later sample (direction is b relative to a).
key – Clonotype match key (default CDR3 aa + V + J).
dispersion (float | str) – "auto" fits the common φ by qCML; a float pins it.
min_total (int) – Clonotypes with a combined equalized count below this are left untested (the discrete tail cannot reach a small p).
log2fc (float) – A clonotype is called only if |log2 FC| >= log2fc and q < alpha. The paper used 5 for vaccine-association; 1 (a doubling) is a gentler default.
alpha (float) – BH-FDR threshold.

Returns:

the key columns, count_a/count_b (raw), q_a/q_b (equalized), log2fc, p_value, q_value, dispersion, and call — one of expanded/contracted/unchanged/untested. Sorted by q_value.

Return type:

One row per clonotype

Biomarker association (`vdjtools.biomarker`)#

Incidence-based clonotype-association testing across a cohort of repertoires (Emerson 2017, Howie 2015, De Witt 2018, Vlasova 2026): feature-vs-condition association (Fisher / χ² / Bayesian / permutation; binary, category, or Cochran–Mantel–Haenszel stratified conditions) and feature-vs-feature co-occurrence (α-β pairing, same-chain co-specificity), with exact or 1-mismatch matching and metaclonotype grouping.

`vdjtools.biomarker.association`#

General incidence association and candidate selection.

General incidence-based biomarker association across a cohort of repertoires.

Generalises the Emerson-2017 Fisher test (vdjtools.biomarker.fisher_association()) along four axes, all sharing the same streamed subject-incidence table:

statistical test — Fisher, χ², a Bayesian log-odds posterior, a Beta-Binomial Bayes factor, or a label permutation (vdjtools.biomarker.stats);
condition type — binary, a category expanded one-vs-rest per level (HLA allele, zygosity), or a paired condition combined by Cochran–Mantel–Haenszel (built with vdjtools.biomarker.condition);
match scope — exact cdr3aa / +v / +v+j (the key), a fuzzy 1-mismatch search (each candidate keeps its identity and gains incidence), or 1mm metaclonotypes (candidates are merged into groups);
candidate set — all public features (min_incidence count or min_incidence_frac fraction of subjects), or an explicit candidates list.

The heavy step (the incidence table) is one streamed group_by over a vdjtools.io.scan_cohort() LazyFrame; every test is then vectorised numpy over the whole feature table. For match="fuzzy", that step is a cohort-wide search independent of the phenotype design — prepare_fuzzy_features() builds it once for reuse across many designs (association(..., features=...)) instead of paying its cost on every call.

class vdjtools.biomarker.association.FeatureFrame(feat_lf, idcols, rep)[source]#

Bases: object

A precomputed association() feature frame, reusable across many designs.

Built once by prepare_fuzzy_features(); pass to association() as features= to skip re-running the fuzzy search when testing many phenotype designs against the same cohort/key/candidates/scope.

Parameters:

feat_lf (LazyFrame)
idcols (list[str])
rep (DataFrame | None)

feat_lf: LazyFrame#

idcols: list[str]#

rep: DataFrame | None#

vdjtools.biomarker.association.prepare_fuzzy_features(cohort, key=('junction_aa', 'v_call', 'j_call'), *, candidates=None, scope='1,0,0,1', productive_only=True, strip_allele=True, threads=0)[source]#

Precompute a match="fuzzy" feature frame once, for reuse across many designs.

The fuzzy search (collecting the cohort’s full per-locus CDR3 universe, then vdjmatch.cluster.overlap()) depends only on (cohort, key, candidates, scope) — never on the phenotype design. A caller testing many designs against the same cohort (e.g. one HLA gene per call, same candidates every time) should build this once and pass it to every association() call via features=, instead of paying the full search cost on every call.

Example

>>> feats = prepare_fuzzy_features(cohort, key=("junction_aa", "v_call"),
...                                candidates=candidates)
>>> for design in designs:
...     association(cohort, design, features=feats)

Parameters:

cohort (LazyFrame | DataFrame)
key (tuple[str, ...])
candidates (DataFrame | None)
scope (str)
productive_only (bool)
strip_allele (bool)
threads (int)

Return type:

FeatureFrame

vdjtools.biomarker.association.select_candidates(cohort, *, key=('junction_aa', 'v_call', 'j_call'), match='exact', min_incidence=2, min_incidence_frac=None, productive_only=True, strip_allele=True, scope='1,0,0,1', threads=0)[source]#

Public features whose subject incidence clears a count and/or fraction threshold.

min_incidence_frac (e.g. 0.05 = 5% of subjects) is resolved against the cohort’s distinct sample_id count and combined with min_incidence by max. Returns the feature key columns (or meta_id + representative key for match="1mm") plus incidence, sorted by descending incidence.

Parameters:

cohort (LazyFrame | DataFrame)
key (tuple[str, ...])
match (str)
min_incidence (int)
min_incidence_frac (float | None)
productive_only (bool)
strip_allele (bool)
scope (str)
threads (int)

Return type:

DataFrame

vdjtools.biomarker.association.association(cohort, phenotype, *, pheno_col=None, level_col=None, stratum_col=None, test='fisher', key=('junction_aa', 'v_call', 'j_call'), match='exact', min_incidence=2, min_incidence_frac=None, candidates=None, alternative='greater', productive_only=True, strip_allele=True, scope='1,0,0,1', n_perm=1000, seed=0, threads=0, features=None)[source]#

Test each clonotype feature’s subject incidence against a condition.

Parameters:

cohort (LazyFrame | DataFrame) – Clonotype cohort — a streamed vdjtools.io.scan_cohort() LazyFrame or a polars.DataFrame with sample_id + the key columns.
phenotype (DataFrame | LazyFrame) – A design frame — one row per subject (or subject × level) with sample_id and either the reserved _pos/_level/_stratum columns (from vdjtools.biomarker.condition) or a plain binary pheno_col.
pheno_col (str | None) – Binary phenotype column (if phenotype has no _pos).
level_col (str | None) – Category-level column → one test per level (adds a level column). Cost is multiplicative, not additive, in the number of levels: the feature-join at the heart of association() duplicates every matched feature row once per level the design carries for that sample (a legitimate consequence of each level needing its own incidence table), so a design with many levels (e.g. one row per HLA allele) can dominate memory even though the feature frame itself costs the same regardless of how many levels the design has.
stratum_col (str | None) – Stratum column → the tests are combined by Cochran–Mantel–Haenszel.
test (str | list[str]) – One test or a list of {"fisher", "chi2", "bayes_logodds", "bayes_bf", "permutation"} (long output, one row per feature×level×test). Ignored when stratum_col is set (CMH is used).
key (tuple[str, ...]) – Feature key (subset of (junction_aa, v_call, j_call)) — the V/J match scope.
match (str) –
"exact", "fuzzy" or "1mm".
- "exact" — the key itself.
- "fuzzy" — a 1-mismatch search: a candidate’s incidence counts every subject carrying anything within scope of it, but the candidate keeps its own identity. Non-junction_aa key columns still have to match exactly, so key=("junction_aa", "v_call") pins the germline half while the CDR3 varies. This is the discovery mode.
- "1mm" — metaclonotypes: candidates within scope are merged and the group is tested (meta_id replaces the key). A downstream step, not discovery — merging dilutes a strong member into its neighbourhood.
min_incidence (int) – Minimum subjects a feature must appear in.
min_incidence_frac (float | None) – Alternative/added fraction-of-subjects threshold (e.g. 0.05).
candidates (DataFrame | None) – Restrict testing to these feature keys (a frame with the key columns).
alternative (str) – "greater" / "less" / "two-sided" (Fisher/permutation).
n_perm (int) – Permutation settings.
seed (int) – Permutation settings.
features (FeatureFrame | None) – A precomputed FeatureFrame from prepare_fuzzy_features(), reused across designs that share the same cohort/key/candidates/scope. When set, key, match, candidates, scope, productive_only and strip_allele are ignored — they are already baked into features.
productive_only (bool)
strip_allele (bool)
scope (str)
threads (int)

Returns:

feature key (or meta_id + representative key + n_members), optional level, then incidence, n_pos_present, n_neg_present, n_pos, n_neg, odds_ratio, log2_or, direction, test and the per-test statistics p_value, q_value (Fisher/χ²/ permutation), logor, logor_ci_lo, logor_ci_hi, p_or_gt1 (Bayesian log-odds), or log_bf (Beta-Binomial). CMH adds or_mh. Sorted by p_value where present.

Return type:

Long frame

`vdjtools.biomarker.cooccurrence`#

Feature-vs-feature co-occurrence (in-silico α-β pairing, same-chain co-specificity).

Feature-vs-feature co-occurrence across a cohort — in-silico α-β pairing & co-specificity.

Two clonotype features that co-occur across subjects more than expected under independence are candidate pairs: a TRA and a TRB from the same clone (Howie 2015 pairSEQ, with subjects playing the role of wells; Vlasova 2026 in-silico α-β pairing), or two same-chain clonotypes that recognise the same antigen (De Witt 2018 co-occurrence patterns). This is the same subject-incidence machinery as vdjtools.biomarker.association(), applied to a pair of features: per pair a 2×2 incidence table over the n subjects profiled for both chains,

has A

n_AB

…

no A

…

…


has A	`n_AB`	…
no A	…	…

is scored by the lift θ = n·n_AB/(n_A·n_B) (Vlasova; observed / expected co-occurrence), Fisher’s exact / χ², and Benjamini-Hochberg FDR. evalue=True adds the expected count and a Poisson upper-tail E-value (the classic control-calibrated co-occurrence significance). The candidate features per chain are bounded by an incidence threshold (and max_features); the returned pairs are those with at least min_cooccurrence co-occurrences.

Repertoire depth is a common cause and is corrected by default. A deep repertoire is more likely to contain any clonotype, so two entirely independent clonotypes co-occur across subjects purely because deep subjects tend to carry both. The induced lift is 1 + CV²(N) for rare clonotypes — independent of biology (0.899 → 1.81 on the FMBA covid19 cohort) — and a pooled Fisher test is badly miscalibrated by it: on simulated independent pairs at the ~11%-incidence regime max_features steers callers into, pooled Fisher declares 45–49% of them significant at p<0.05. depth_strata (default 10) therefore scores each pair by Cochran–Mantel–Haenszel across equal-count per-subject depth strata, which restores calibration (measured false-positive rate 0.023–0.057) at equal power. depth_strata=0 restores the pooled test. Note this corrects depth only — shared HLA, ancestry, batch and exposure also induce cross-subject co-occurrence, and none of them is fixed here; see the warning in docs/usage.rst.

vdjtools.biomarker.cooccurrence.cooccurrence(cohort, *, chain_a='TRA', chain_b='TRB', key=('junction_aa', 'v_call', 'j_call'), match='exact', test='fisher', min_incidence=2, min_incidence_frac=None, min_cooccurrence=2, candidates_a=None, candidates_b=None, evalue=False, alternative='greater', depth_strata=10, max_features=2000, productive_only=True, strip_allele=True, scope='1,0,0,1', threads=0)[source]#

Score co-occurrence of feature pairs across the subjects profiled for both chains.

Parameters:

cohort (LazyFrame | DataFrame) – Clonotype cohort with a sample_id column (both chains present per subject for α-β pairing).
chain_a (str) – Loci to pair (e.g. "TRA"/"TRB"). chain_b=None (or equal to chain_a) does same-chain pairs (upper triangle, self-pairs excluded).
chain_b (str | None) – Loci to pair (e.g. "TRA"/"TRB"). chain_b=None (or equal to chain_a) does same-chain pairs (upper triangle, self-pairs excluded).
key (tuple[str, ...]) – Feature definition and exact/1mm scope (as in association()).
match (str) – Feature definition and exact/1mm scope (as in association()).
test (str) – "fisher" or "chi2" for the per-pair 2×2 p-value. Ignored when depth conditioning is active (CMH is used), mirroring association()’s stratified branch. It is used when the guards below make conditioning fall back to pooled.
alternative (str) – "greater" (default; co-occurrence — the usual question), "less" (mutual exclusion) or "two-sided". Shapes the Fisher p (and, under CMH, the halving of the two-sided χ² toward or_mh). It does not shape a plain test="chi2" p on the pooled branch: Pearson χ² is inherently two-sided, so a depth-conditioning fallback with test="chi2" reports the same two-sided p regardless of alternative — read direction/odds_ratio for the sign there.
min_incidence (int) – Candidate-feature incidence threshold per chain.
min_incidence_frac (float | None) – Candidate-feature incidence threshold per chain.
min_cooccurrence (int) – Keep only pairs co-occurring in ≥ this many subjects.
candidates_a (DataFrame | None) – Restrict each chain’s features to these keys.
candidates_b (DataFrame | None) – Restrict each chain’s features to these keys.
evalue (bool) – Also report expected (= n_A·n_B/n) and a Poisson upper-tail e_value.
depth_strata (int) – Number of equal-count per-subject depth strata to condition on (Cochran–Mantel–Haenszel). Depth is a common cause of co-occurrence, so the pooled test is anticonservative; the default corrects it. 0 → pooled test (uncorrected; kept as the oracle and the De Witt-comparable path). Subject depth is the number of distinct key features that subject contributes to the analysed chains — derived from cohort, never supplied, so it cannot disagree with the data. Conditioning needs both a reason and enough data, so it falls back to the pooled test (with a warning) when depth barely varies (CV < 0.45, i.e. an induced lift under ~1.2) or the cohort cannot fill ≥20 subjects per stratum (<40 subjects). When that happens or_mh/chi2 are absent and test applies.
max_features (int) – Cap candidate features per chain (top by incidence); logs if it truncates.
productive_only (bool)
strip_allele (bool)
scope (str)
threads (int)

Returns:

a_<key> / b_<key> columns, n (subjects with both chains), n_a, n_b, n_ab, theta (lift), odds_ratio, log2_or, p_value, q_value (BH), and — with evalue — expected, e_value. With depth_strata > 0 also or_mh and chi2 (the CMH estimate the p_value comes from). Sorted by p_value.

theta/odds_ratio remain the pooled, depth-uncorrected lift — compare them against or_mh, which is conditioned on depth.

Return type:

One row per surviving pair

`vdjtools.biomarker.condition`#

Phenotype-design builders (binary, categorical, HLA alleles, zygosity, CMH strata).

Turn subject metadata into the design frame the association engine consumes.

A design frame has one row per subject (or per subject × level) with the reserved columns

_pos — boolean, the condition-positive indicator for that (subject, level);
_level — optional, the category level (HLA allele, homo/hetero, …); absent ⇒ a single binary test;
_stratum — optional, the Cochran–Mantel–Haenszel stratum; absent ⇒ unstratified.

vdjtools.biomarker.association() reads exactly these columns, so any of the builders below (or a hand-built frame with the same columns) can drive it. All keep the subject key column sample_id.

vdjtools.biomarker.condition.binary(meta, col, *, sample_col='sample_id')[source]#

Design frame for a binary condition (e.g. CMV +/−); unknown labels dropped.

Parameters:

meta (DataFrame)
col (str)
sample_col (str)

Return type:

DataFrame

vdjtools.biomarker.condition.categorical(meta, col, *, min_level_size=1, sample_col='sample_id')[source]#

One-vs-rest design for a single-label category (each subject in exactly one level).

Emits the full subject × level cross with _pos = (col == level). Levels carried by fewer than min_level_size subjects are dropped. Null category values are excluded.

Parameters:

meta (DataFrame)
col (str)
min_level_size (int)
sample_col (str)

Return type:

DataFrame

vdjtools.biomarker.condition.hla_alleles(meta, cols, *, resolution=None, min_level_size=1, sample_col='sample_id')[source]#

Multi-label design over HLA alleles: per allele, carriers vs non-carriers.

cols are the allele columns for a locus (e.g. ["HLA-A.1", "HLA-A.2"] or the 4-digit ["sample.HLA-A.1", "sample.HLA-A.2"]). A subject is _pos for every allele it carries. resolution trims the allele to that many colon-separated fields (resolution=1 → A*02); None keeps the field as written. Alleles carried by fewer than min_level_size subjects are dropped.

HLA-untyped subjects are dropped, not counted as non-carriers — consistent with binary() / zygosity(), which drop unknown phenotypes. Treating an untyped subject as a non-carrier of every allele silently inflates the negative arm and biases the odds ratio anticonservatively.

Parameters:

meta (DataFrame)
cols (list[str])
resolution (int | None)
min_level_size (int)
sample_col (str)

Return type:

DataFrame

vdjtools.biomarker.condition.zygosity(meta, locus_cols, *, sample_col='sample_id')[source]#

Binary homozygous(_pos=True)/heterozygous design for a locus’s two allele columns.

Parameters:

meta (DataFrame)
locus_cols (tuple[str, str])
sample_col (str)

Return type:

DataFrame

vdjtools.biomarker.condition.stratified(meta, pheno_col, stratum_col, *, sample_col='sample_id')[source]#

Design for a paired condition: binary pheno_col stratified by stratum_col (CMH).

E.g. CMV association conditioned on an HLA group. Subjects with an unknown phenotype or stratum are dropped.

Parameters:

meta (DataFrame)
pheno_col (str)
stratum_col (str)
sample_col (str)

Return type:

DataFrame

`vdjtools.biomarker.stats`#

Vectorised 2×2 test kernels (Fisher, χ², Bayesian, CMH, permutation, FDR).

Vectorised 2×2 contingency-test kernels for biomarker association / co-occurrence.

Every kernel takes the four cells of the subject-incidence table as numpy int arrays (one entry per feature / feature-pair) and returns numpy arrays — there is no per-feature Python loop, so millions of features score in a handful of scipy calls. The cell convention throughout is

feature present

a

b

feature absent

c

d


feature present	`a`	`b`
feature absent	`c`	`d`

so a = subjects that are condition-positive and carry the feature, n = a+b+c+d, n_pos = a+c, n_neg = b+d, present = a+b. Kernels: Fisher (hypergeometric tail), Pearson χ² (Yates), Haldane odds ratio, a normal-approximation Bayesian posterior over the log odds-ratio, a Beta-Binomial Bayes factor, Cochran–Mantel–Haenszel (stratified), a label-permutation null, and Benjamini-Hochberg FDR.

vdjtools.biomarker.stats.fisher_p(a, b, c, d, alternative='greater')[source]#

Fisher’s exact p-value per feature, as the hypergeometric tail (vectorised).

"greater" tests enrichment in condition+ (one-tailed, the CMV setting), "less" depletion, "two-sided" either (doubled smaller tail, capped at 1) — identical to the convention in vdjtools.biomarker.fisher_association().

"two-sided-minlike" is the minimum-likelihood two-sided convention used by R’s fisher.test and scipy.stats.fisher_exact(): the sum of the likelihood of every table at most as likely as the observed one. It differs from "two-sided" whenever the margins are not exactly symmetric — measured on paired-dynamics-shaped tables, exactly equal margins agree 300/300, but margins differing by rounding alone disagree 171/300 with "two-sided" up to exactly 2× larger. Use it when comparing against R or scipy, or when the margins are only near-equal; "two-sided" is kept as the default two-sided convention because the association/co-occurrence results were computed with it.

Parameters:: alternative (str)
Return type:: ndarray

vdjtools.biomarker.stats.chi2_p(a, b, c, d, yates=True)[source]#

Pearson χ² p-value per feature for the 2×2 table (Yates-corrected by default).

Matches scipy.stats.chi2_contingency([[a,b],[c,d]], correction=yates). Any table with an empty margin yields χ²=0, p=1.

Parameters:: yates (bool)
Return type:: ndarray

vdjtools.biomarker.stats.odds_ratio(a, b, c, d)[source]#

Haldane–Anscombe (+0.5) odds ratio per feature — never 0/∞.

Return type:: ndarray

vdjtools.biomarker.stats.direction(a, b, c, d)[source]#

"enriched" if the feature is over-represented in condition+, else "depleted".

Return type:: ndarray

vdjtools.biomarker.stats.bayes_logodds(a, b, c, d, alpha=0.5)[source]#

Normal-approximation Bayesian posterior over the log odds-ratio (Woolf).

With a flat prior on log OR the posterior is Normal(logor, se²) where logor is the Haldane-corrected log odds-ratio and se² = Σ 1/(cell+alpha) (Woolf standard error; alpha=0.5 = Haldane). Returns logor, se, a 95% credible interval (ci_lo/ci_hi), and p_or_gt1 = P(OR>1) = Φ(logor/se).

Parameters:: alpha (float)
Return type:: dict[str, ndarray]

vdjtools.biomarker.stats.bayes_bf(a, b, c, d, alpha=1.0, beta=1.0)[source]#

log Bayes factor BF₁₀ for association, Beta-Binomial (analytic, vectorised).

Compares H₁ (the feature’s presence rate differs between conditions: independent Beta(alpha,beta) priors on the two rates) against H₀ (a single shared rate). The binomial coefficients cancel, leaving

log BF₁₀ = betaln(a+α, c+β) + betaln(b+α, d+β) − betaln(α,β) − betaln(a+b+α, c+d+β)

(positive ⇒ evidence for association). alpha=beta=1 is a uniform prior.

Parameters:

alpha (float)
beta (float)

Return type:

ndarray

vdjtools.biomarker.stats.cmh(a, b, c, d)[source]#

Cochran–Mantel–Haenszel stratified test — combined OR + χ² over strata.

a,b,c,d are 2-D arrays of shape (n_features, n_strata) (the same 2×2 table per feature, one column per stratum). Returns the Mantel-Haenszel odds ratio or_mh, the continuity-corrected chi2 statistic, and its 1-df p_value. Strata with fewer than 2 subjects (n_k < 2) contribute nothing (their variance is undefined).

Return type:: dict[str, ndarray]

vdjtools.biomarker.stats.permutation_p(present, labels, *, n_perm=1000, seed=0, alternative='greater')[source]#

Label-permutation p-value per feature from the subject-incidence matrix.

present is a boolean (n_subjects, n_features) incidence matrix, labels a boolean (n_subjects,) condition vector. Permuting labels holds each feature’s incidence fixed, so the null is exactly the hypergeometric one — this converges to fisher_p() as n_perm → ∞ and is a robust check for sparse tables. Uses the add-one estimator (1 + #exceed)/(1 + n_perm). Seeded → reproducible.

A permutation’s counts are L @ present with L the 0/1 permuted-label matrix, so a chunk of permutations is one BLAS call instead of n_perm fancy-index copies (~15x). The chunk still draws rng.permutation(n_sub) once per permutation in the original order, so p-values are unchanged for a given seed; and summing 0/1 to an integer <= n_sub is exact in float64 whatever order BLAS accumulates in.

Parameters:

present (ndarray)
labels (ndarray)
n_perm (int)
seed (int)
alternative (str)

Return type:

ndarray

vdjtools.biomarker.stats.fdr_bh(p)[source]#

Benjamini-Hochberg FDR-adjusted q-values (empty-safe).

Parameters:: p (ndarray)
Return type:: ndarray

`vdjtools.biomarker.fisher`#

Fisher-exact incidence association (V/J-match, exact / 1mm) — the Emerson-2017 shortcut.

Incidence-based Fisher association — the Emerson-2017 method (back-compat entry point).

This is now a thin wrapper over vdjtools.biomarker.association() with test="fisher"; it exists so the original column schema and call signature keep working. New code that wants other tests (χ², Bayesian, permutation), category / stratified conditions, candidate selection, or co-occurrence should call vdjtools.biomarker.association() / vdjtools.biomarker.cooccurrence() directly.

For each clonotype feature, a 2×2 subject-incidence table (present/absent × phenotype±) is tested with Fisher’s exact test (Emerson et al., Nat Genet 2017, doi:10.1038/ng.3822).

vdjtools.biomarker.fisher.fisher_association(cohort, phenotype, *, pheno_col, key=('junction_aa', 'v_call', 'j_call'), match='exact', min_incidence=2, alternative='greater', productive_only=True, strip_allele=True, scope='1,0,0,1', threads=0)[source]#

Test each clonotype feature’s subject incidence against a binary phenotype (Fisher).

See vdjtools.biomarker.association() for the arguments (this fixes test="fisher").

Returns:

One row per tested feature — the key columns (or meta_id + a representative key and n_members for match="1mm") followed by incidence, n_pos_present, n_neg_present, n_pos, n_neg, odds_ratio, log2_or, p_value, q_value (Benjamini-Hochberg) and direction, sorted by p_value.

Parameters:

cohort (LazyFrame | DataFrame)
phenotype (DataFrame | LazyFrame)
pheno_col (str)
key (tuple[str, ...])
match (str)
min_incidence (int)
alternative (str)
productive_only (bool)
strip_allele (bool)
scope (str)
threads (int)

Return type:

DataFrame

`vdjtools.biomarker.metaclonotype`#

Metaclonotype grouping (fuzzy CDR3 + V/J).

Metaclonotype grouping — collapse edit-scope-neighbour clonotypes into one feature.

Emerson-style incidence association tests each exact public TCRβ. A metaclonotype instead groups clonotypes whose CDR3s are within an edit scope (default one substitution) — optionally requiring the same V and/or J — so a family of near-variants counts as one biomarker feature. The fuzzy search is not reimplemented: it is delegated to vdjmatch.cluster.overlap() (native seqtree engine), exactly as vdjtools.overlap.fuzzy does. Connected components of the within-scope neighbour graph (single-linkage, union-find) become the metaclonotypes.

Scale: the unique clonotype keys are clustered once (not per subject). Partitioning by V and/or J means the all-pairs search only ever runs within one gene group, and the native call releases the GIL — so ~1M unique CDR3s are grouped in a handful of multi-threaded passes.

vdjtools.biomarker.metaclonotype.metaclonotypes(clonotypes, *, scope='1,0,0,1', match_v=True, match_j=True, threads=0)[source]#

Group unique clonotype keys into metaclonotypes by CDR3 edit-scope neighbourhood.

Two keys share a meta_id iff their junction_aa are within scope and they share a V call (when match_v) and a J call (when match_j). Grouping is single-linkage (connected components of the neighbour graph).

Parameters:

clonotypes (DataFrame) – A clonotype frame; must carry junction_aa and, when the corresponding match_* flag is set, v_call / j_call.
scope (str) – vdjmatch edit-distance scope "subs,ins,dels,total" (default one substitution, length-preserving). "0,0,0,0" reduces to exact grouping.
match_v (bool) – Require the same v_call for two keys to be grouped.
match_j (bool) – Require the same j_call for two keys to be grouped.
threads (int) – Worker threads for the native search (0 = all cores).

Returns:

The distinct grouping keys (junction_aa plus v_call/j_call as applicable) with an added meta_id column (compact 0-based integer, singletons included).

Raises:

ImportError – If vdjmatch is not importable (it is a base dependency).

Return type:

DataFrame

Single-cell interop (`vdjtools.sc`)#

Single-cell AIRR Cell / 10x paired-chain interop: contig ingestion, chain resolution and pairing, doublet / mispairing QC, cluster evaluation, and bridges to scirpy, dandelion and scRepertoire. See Single-cell for the narrative guide.

`vdjtools.sc.read`#

10x / AIRR-Cell / arda ingestion and AIRR Data File export.

10x / AIRR-Cell ingestion → a long Rearrangement frame keyed on cell_id.

The load-bearing single-cell representation here is deliberately flat: one row per productive contig, canonical AIRR Rearrangement columns, plus a cell_id (the cell barcode) tying a cell’s chains together and umi_count / clone_id carried alongside. Everything downstream (pairing, QC, clustering) consumes that frame; the AIRR Data File export (write_airr_cell()) is a secondary interchange layer.

10x CellRanger all_contig_annotations.csv names the junction column cdr3 (with the conserved Cys/Phe-Trp anchors included) — content-identical to our junction_aa / AIRR junction_aa convention (see io.schema) — and cdr3_nt for the nucleotide junction; the reader maps those straight onto junction_aa / junction_nt.

vdjtools.sc.read.VALID_LOCI: tuple[str, ...] = ('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL')#: Receptor loci 10x can call; Multi / None / anything else is dropped.

vdjtools.sc.read.SC_COLUMNS: list[str] = ['cell_id', 'sequence_id', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'junction_aa', 'junction_nt', 'duplicate_count', 'umi_count', 'clone_id', 'productive']#

Canonical single-cell long-frame columns, in order.

productive is a mandatory AIRR Rearrangement field and is carried so vdjtools.sc.to_airr() can emit a schema-valid table; the readers here keep only productive contigs, so it is True wherever the source stated productivity at all.

vdjtools.sc.read.read_10x(all_contig, consensus=None, *, require_cell=True, require_high_conf=True)[source]#

Read a 10x contig-annotation CSV into the canonical sc long frame.

Accepts both all_contig_annotations.csv[.gz] and filtered_contig_annotations.csv[.gz] — the two share one writer in CellRanger and so one column layout; filtered_ is simply pre-restricted to is_cell && high_confidence, which this reader applies anyway.

Keeps only productive, cell-associated, high-confidence contigs on a real receptor locus (VALID_LOCI) with a resolved consensus id, one row per contig. When a consensus_annotations file is supplied, the per-cell contig’s V/D/J calls are replaced by the matched consensus calls (joined on (raw_clonotype_id, raw_consensus_id) == (clonotype_id, consensus_id)).

Column drift across CellRanger versions is tolerated: the fwr*/cdr1/cdr2 region columns exist only from CR6, exact_subclonotype_id from CR4+, and sample only in cellranger multi output — none is required here. raw_consensus_id is used when present and skipped when absent.

NOTE: For CellRanger 4.0+, prefer read_airr_cell() on airr_rearrangement.tsv — it is the same file the downstream tools read, so ingestion and interop agree.

Parameters:

all_contig (str | Path) – Path to all_contig_annotations.csv or filtered_contig_annotations.csv (.gz accepted). Columns follow CellRanger VDJ (barcode, is_cell, contig_id, high_confidence, chain, v_gene, d_gene, j_gene, c_gene, productive, cdr3, cdr3_nt, reads, umis, raw_clonotype_id, raw_consensus_id). *_call spellings are also accepted in place of *_gene.
consensus (str | Path | None) – Optional consensus_annotations.csv[.gz] to source consensus V/D/J calls from; if None the contig’s own calls are used.
require_cell (bool) – Drop contigs whose is_cell is not truthy (default True).
require_high_conf (bool) – Drop contigs whose high_confidence is not truthy (default True).

Returns:

A pl.DataFrame in the canonical sc long-frame layout (SC_COLUMNS) — one row per surviving productive contig.

Raises:

ValueError – If a required column (barcode, contig_id, chain, cdr3) is missing from all_contig.

Return type:

DataFrame

vdjtools.sc.read.read_airr_cell(path)[source]#

Read an AIRR Rearrangement TSV carrying a cell_id column into the sc frame.

This is also the reader for CellRanger’s ``airr_rearrangement.tsv`` (emitted since Cell Ranger 4.0, and under per_sample_outs/<sample>/vdj_t/ for cellranger multi) and for arda’s barcoded output — both are plain AIRR Rearrangement tables with a cell_id column. Prefer it over read_10x(): it is the same file scirpy, dandelion and scRepertoire read, so the ingestion and interop paths agree by construction.

Parameters:: path (str | Path) – Path to an AIRR Rearrangement TSV (.gz accepted) with at least a cell_id column plus the usual AIRR fields.
Returns:: A pl.DataFrame in the canonical sc long-frame layout (SC_COLUMNS); columns absent from the file are filled with nulls. junction_aa / junction are accepted as sources for junction_aa / junction_nt.
Raises:: ValueError – If the file has no cell_id column.
Return type:: DataFrame

vdjtools.sc.read.ARDA_STATUS = 'arda_status'#: Cell-level QC columns lifted from arda’s .chains.tsv, namespaced so they cannot be confused with vdjtools.sc.resolve_chains()’ own verdict.

vdjtools.sc.read.read_arda_cells(prefix, *, chains=True)[source]#

Read the output of arda’s single-cell pipeline (arda cells).

arda cells writes <prefix>.contigs.airr.tsv — an AIRR Rearrangement table with cell_id, molecules and reads — plus a cell-level <prefix>.chains.tsv carrying arda’s own per-chain verdict (status is one of primary, secondary, doublet_candidate, extra).

NOTE: arda’s status and resolve_chains()’ verdict are independent calls on the same question. This reader surfaces arda’s as arda_status rather than acting on it, so running resolve_chains afterwards does not silently discard the upstream judgement – compare them, don’t assume they agree.

Reading goes through arda’s own read_airr, not a plain CSV read: arda has written two TSV dialects and that reader is what normalises them. It matters because junction_quality is Phred+33 and character 34 is a double quote.

Parameters:

prefix (str | Path) – The arda cells output prefix (<prefix>.contigs.airr.tsv is read), or a direct path to a *.airr.tsv file.
chains (bool) – Join arda’s per-chain status / molecules from <prefix>.chains.tsv when that file exists (default True).

Returns:

A pl.DataFrame in the canonical layout (SC_COLUMNS), plus arda_status / arda_molecules when the chains table was joined.

Raises:

FileNotFoundError – If no contigs AIRR table is found for prefix.

Return type:

DataFrame

vdjtools.sc.read.write_airr_cell(rearr, cells_out, *, receptors=True, repertoire_id='')[source]#

Emit an AIRR Data File (YAML) with a Cell array (and optional Receptor).

Builds one Cell per cell_id (linking its sequence_id receptors) and, when receptors is set, one Receptor per paired heavy/light chain within a cell. The receptor_hash is sha256 of the two upper-cased domain sequences.

Note

The AIRR spec’s receptor_variable_domain_{1,2}_aa is the full mature V-domain amino-acid sequence. 10x contigs only expose the junction, so this field is populated with the junction (junction_aa) and that limitation is recorded in the file’s Info block. Downstream code should treat these as junction-level, not full-domain, sequences.

Parameters:

rearr (DataFrame) – A single-cell long frame (SC_COLUMNS), typically the paired output — one cell_id may carry several chains.
cells_out (str | Path) – Destination path for the AIRR Data File (.yaml / .json).
receptors (bool) – Emit the Receptor list pairing heavy (β/heavy) and light (α/light) chains per cell (default True).
repertoire_id (str) – Value stamped into each Cell.repertoire_id (default empty).

Returns:

The Path written.

Raises:

ImportError – If PyYAML is not installed (see the sc extra).

Return type:

`vdjtools.sc.airr`#

The AIRR Rearrangement interchange layer every downstream bridge is built on, plus the scRepertoire export.

The single-cell interchange layer: a flat AIRR Rearrangement table.

scirpy, dandelion and scRepertoire all read the same thing — one AIRR Rearrangement row per contig, carrying sequence_id and cell_id. None of them consumes AIRR Cell objects; cell-level state lives in adata.obs / Dandelion.metadata / Seurat meta.data. So the interop surface here is one emitter (to_airr()) and one inverse (from_airr()); every bridge in vdjtools.sc is a thin adapter on top, and the AIRR Cell Data File (vdjtools.sc.write_airr_cell()) stays a separate, spec-faithful export rather than an interop path.

Two spellings are reconciled here, in exactly one place:

vdjtools calls the nucleotide junction junction_nt; AIRR — and therefore every downstream tool — calls it junction.
consensus_count is what scRepertoire’s AIRR parser reads as the read count, while scirpy and dandelion prefer umi_count. Both are emitted, so one file feeds all three.

vdjtools.sc.airr.JUNCTION = 'junction'#: AIRR spelling of the nucleotide junction (vdjtools stores it as junction_nt).

vdjtools.sc.airr.CONSENSUS_COUNT = 'consensus_count'#: Read count backing the contig consensus; scRepertoire’s AIRR parser reads this one.

vdjtools.sc.airr.AIRR_COLUMNS: list[str] = ['sequence_id', 'cell_id', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'junction', 'junction_aa', 'productive', 'duplicate_count', 'umi_count', 'consensus_count', 'clone_id']#: Columns of the emitted AIRR Rearrangement table, in order.

vdjtools.sc.airr.to_airr(cells)[source]#

Convert the canonical sc long frame to a flat AIRR Rearrangement table.

This is the contract every downstream bridge is built on. sequence_id is synthesised as <cell_id>_contig_<n> when the frame does not carry one — dandelion derives cell_id back out of exactly that pattern when a file omits it, so the round-trip closes either way.

consensus_count is set from duplicate_count (the read count): for a 10x contig the reads backing the consensus are the duplicate observations, and emitting both spellings is what lets a single file satisfy scRepertoire (which reads consensus_count) and scirpy/dandelion (which prefer umi_count) at once.

Parameters:: cells (DataFrame) – Single-cell long frame — vdjtools.sc.read.SC_COLUMNS, as returned by read_10x() or read_airr_cell(). Missing optional columns are emitted as nulls.
Returns:: A pl.DataFrame with columns AIRR_COLUMNS, one row per contig.
Raises:: ValueError – If cells has no cell_id column.
Return type:: DataFrame

vdjtools.sc.airr.from_airr(airr)[source]#

Convert a flat AIRR Rearrangement table back to the canonical sc long frame.

The inverse of to_airr(), and the shared tail of every from_* bridge (from_scirpy(), from_dandelion()) — they each reduce their container to an AIRR frame and hand it here.

Parameters:: airr (DataFrame) – AIRR Rearrangement frame with at least cell_id. Both the AIRR junction and the vdjtools junction_nt spellings are accepted, as are consensus_count / duplicate_count for the read count.
Returns:: A pl.DataFrame in the canonical layout (vdjtools.sc.read.SC_COLUMNS).
Raises:: ValueError – If airr has no cell_id column.
Return type:: DataFrame

vdjtools.sc.airr.write_airr(cells, path)[source]#

Write the sc frame as an AIRR Rearrangement TSV.

The result is readable by read_airr_cell(), by scirpy (ir.io.read_airr), by dandelion (ddl.read_airr) and by scRepertoire (loadContigs(format="AIRR")).

Parameters:

cells (DataFrame) – Single-cell long frame (vdjtools.sc.read.SC_COLUMNS).
path (str | Path) – Destination .tsv.

Returns:

The path written.

Return type:

vdjtools.sc.airr.write_screpertoire(cells, path, *, format='airr')[source]#

Write a file R’s scRepertoire can load, for combineTCR / combineBCR.

Export only — no R code ships with vdjtools. Load the result with scRepertoire::loadContigs(path, format = "AIRR") (or "10X").

NOTE: two footguns on the R side, neither of which this function can prevent. .parseAIRR reads the read count from ``consensus_count``, not umi_count — which is why to_airr() emits both. And combineTCR(samples=, ID=) rewrites barcodes to sample_ID_barcode, the usual cause of a silent barcode-join failure against a Seurat meta.data; pass the same samples/ID there and to combineExpression, or omit both.

Parameters:

cells (DataFrame) – Single-cell long frame (vdjtools.sc.read.SC_COLUMNS).
path (str | Path) – Destination — .tsv for format="airr", .csv for format="10x".
format (str) – "airr" writes an airr_rearrangement.tsv-shaped file (the columns .parseAIRR requires); "10x" writes a filtered_contig_annotations.csv-shaped file.

Returns:

The path written.

Raises:

ValueError – If format is not "airr" or "10x".

Return type:

`vdjtools.sc.dandelion`#

dandelion bridge, including a dandelion-free .h5ddl reader.

Bridge the single-cell frame to and from dandelion.

A Dandelion object is two tables: .data, the contig-level AIRR frame indexed by sequence_id, and .metadata, one row per cell. The first is exactly what to_airr() emits, so the bridge is a format hand-off rather than a translation — dandelion’s own scirpy converter round-trips through the same flat AIRR table.

read_h5ddl() is the useful part: dandelion persists to .h5ddl, which is plain HDF5 (h5py structured arrays, with a sibling Zarr store for distances), so a dandelion result is readable without installing dandelion.

NOTE: dandelion also ships a polars backend (ddl.set_backend("polars"), DandelionPolars) whose .data/.metadata take polars frames directly. If you are on that backend, hand to_airr() output straight over; the pandas conversion here is only for the default backend.

vdjtools.sc.dandelion.to_dandelion(cells)[source]#

Wrap the sc long frame as a dandelion.Dandelion.

sequence_id and umi_count are both required by dandelion’s loader; to_airr() guarantees the first (synthesising <cell>_contig_<n> when absent) and carries the second.

Parameters:: cells (DataFrame) – Single-cell long frame (vdjtools.sc.read.SC_COLUMNS).
Returns:: A Dandelion whose .data is the contig-level AIRR table; .metadata is built by dandelion’s own update_metadata().
Raises:: ImportError – If dandelion is not installed (pip install sc-dandelion).

vdjtools.sc.dandelion.from_dandelion(vdj)[source]#

Convert a dandelion.Dandelion back to the canonical sc long frame.

Parameters:: vdj – A Dandelion object (or anything exposing a contig-level .data frame).
Returns:: A pl.DataFrame in the canonical layout (vdjtools.sc.read.SC_COLUMNS).
Raises:: AttributeError – If vdj has no .data.
Return type:: DataFrame

vdjtools.sc.dandelion.read_h5ddl(path)[source]#

Read the contig table out of a dandelion .h5ddl file, without dandelion.

.h5ddl is plain HDF5: data is a h5py structured array of the contig-level AIRR table (metadata holds the cell-level one, and distances may live in a sibling Zarr store). Only h5py is needed to read it.

Parameters:

path (str | Path) – Path to a .h5ddl file.

Returns:

A pl.DataFrame in the canonical layout (vdjtools.sc.read.SC_COLUMNS).

Raises:

ImportError – If h5py is not installed.
KeyError – If the file has no data dataset.

Return type:

DataFrame

`vdjtools.sc.pair`#

Chain resolution, alpha/beta pairing and mispairing flags.

Chain resolution, paired-receptor assembly, and doublet / mispairing QC.

A 10x cell’s contigs are noisy: barcodes collide (doublets), ambient mRNA leaks a spurious light chain, and a genuine α/β cell can legitimately carry two productive α. These functions clean that up on the flat single-cell frame from vdjtools.sc.read.

Locus roles (heavy / light) follow the standard receptor families:

heavy   = {TRB, TRD, IGH}        # one per cell
light   = {TRA, TRG}             # one, sometimes two (dual-α)
B-light = {IGK, IGL}             # one, sometimes two

The per-cell chain ranking key is (-duplicate_count, -umi_count, sequence_id) everywhere: most reads first, then most UMIs, then a stable id tie-break. The thresholds encode the mirpy-derived rule one heavy but possibly two light, and are reimplemented here on polars (no mirpy dependency).

vdjtools.sc.pair.LOCUS_PAIR_TO_LOCI: dict[str, tuple[str, str]] = {'IGH_IGK': ('IGK', 'IGH'), 'IGH_IGL': ('IGL', 'IGH'), 'TRA_TRB': ('TRA', 'TRB'), 'TRG_TRD': ('TRG', 'TRD')}#: locus-pair family -> (light/chain1 locus, heavy/chain2 locus).

vdjtools.sc.pair.resolve_chains(rearr, *, secondary_ratio=0.1, secondary_min_umi=2, secondary_min_dup=5)[source]#

Reduce over-expanded per-cell chains to one heavy and one (or two) light.

Per cell_id:

keep exactly the top-1 heavy chain (TRB / TRD / IGH);
keep the top-1 light chain, and a second light chain only when all of second_dup/first_dup > secondary_ratio, second_umi/first_umi > secondary_ratio, second_umi >= secondary_min_umi and second_dup >= secondary_min_dup hold (the dual-α allowance);
the same secondary rule applies jointly across IGK + IGL for B-cells.

Parameters:

rearr (DataFrame) – Single-cell long frame (vdjtools.sc.read.SC_COLUMNS).
secondary_ratio (float) – Minimum second/first ratio (on both reads and UMIs) to admit a second light chain.
secondary_min_umi (int) – Minimum absolute UMI count for a second light chain.
secondary_min_dup (int) – Minimum absolute read count for a second light chain.

Returns:

The cleaned per-cell contigs (same columns as the input), ordered by cell then rank. Contigs on loci outside the receptor roles are dropped.

Return type:

DataFrame

vdjtools.sc.pair.pair_chains(rearr, *, locus_pair='TRA_TRB', resolve=True)[source]#

Assemble paired receptors as the Cartesian product of a cell’s light × heavy.

After (optionally) resolve_chains(), each cell forms one paired receptor per (light, heavy) combination of its chains in the requested family — so a cell with two α and one β yields two pairs (<cell>_1, <cell>_2). Cells missing either side of the family are counted but not emitted (see chain_multiplicity()).

Parameters:

rearr (DataFrame) – Single-cell long frame.
locus_pair (str) – Family to pair — one of "TRA_TRB", "TRG_TRD", "IGH_IGK", "IGH_IGL". The first locus is the α/light side (alpha_* columns), the second the β/heavy side (beta_* columns).
resolve (bool) – Run resolve_chains() first (default True).

Returns:

One row per paired receptor with cell_id, pair_id, alpha_v_call, alpha_j_call, alpha_junction_aa, alpha_umi_count, alpha_duplicate_count and the matching beta_* columns.

Raises:

ValueError – If locus_pair is not a recognised family.

Return type:

DataFrame

vdjtools.sc.pair.chain_multiplicity(rearr, *, locus_pair='TRA_TRB')[source]#

Presence-quadrant histogram (n_light, n_heavy) -> cell_count for a family.

Counts, over cells, how many carry each (n_light, n_heavy) combination of chain multiplicities in locus_pair — the α/β quadrant table used to diagnose doublets and dropout. Cells with neither chain in the family contribute a (0, 0) row.

Parameters:

rearr (DataFrame) – Single-cell long frame.
locus_pair (str) – Family to tabulate (see pair_chains()).

Returns:

A pl.DataFrame with columns n_light, n_heavy, cell_count, sorted by n_light then n_heavy.

Raises:

ValueError – If locus_pair is not a recognised family.

Return type:

DataFrame

vdjtools.sc.pair.flag_mispairing(paired, *, max_slaves_per_master=None, drop=False)[source]#

Flag suspected mispaired / ambient α chains against a master(β) → slave(α) graph.

Builds, across all cells, how often each master (β) heavy chain co-occurs with each slave (α) light chain. For every master its canonical slave is the one with the most co-occurrences (ties broken by summed read+UMI support). Any pairing whose α is not the master’s canonical slave is flagged as suspected mispairing / contamination. If a master pairs with more than max_slaves_per_master distinct α across the dataset, the master itself is flagged as ambient (a β smeared across too many barcodes).

Chains are keyed on (v_call, j_call, junction_aa) per side, so identical clonotypes across cells are recognised as the same master / slave.

Parameters:

paired (DataFrame) – Output of pair_chains() (alpha_* / beta_* columns).
max_slaves_per_master (int | None) – Distinct-α ceiling above which a master is called ambient; None disables the ambient check.
drop (bool) – If True, remove flagged rows instead of annotating them.

Returns:

The paired frame plus mispairing_flag (bool) and mispairing_reason ("ok" / "noncanonical_alpha" / "ambient_master"). When drop is set, flagged rows are removed and the two columns omitted.

Raises:

ValueError – If max_slaves_per_master is not a positive integer.

Return type:

DataFrame

`vdjtools.sc.pgen`#

Paired-chain generation probability (Pgen(α)·Pgen(β)) via the native model.

Paired-chain generation probability for single-cell repertoires.

Under chain independence the paired generation probability of a cell is Pgen(α) · Pgen(β) — the product of each chain’s junction Pgen under the native vdjtools.model engine (bundled per-locus models). This is the single-cell paired-Pgen residual from Phase 7; it is computed entirely from the native model (no vdjmatch dependency).

The paired frame is the vdjtools.sc.pair.pair_chains() layout — alpha_v_call, alpha_j_call, alpha_junction_aa and the beta_* counterparts (α/light and β/heavy). Each chain’s locus is inferred from its V-call prefix (TRA/TRB, or IGK/IGL + IGH for BCR) unless given explicitly.

Conditioning on V/J requires the call to match a model allele (e.g. TRBV20-1*01); a gene-level or unmatched call marginalises over all V/J for that chain (still a valid, if less specific, Pgen). Pass condition_vj=False to marginalise unconditionally.

vdjtools.sc.pgen.paired_pgen(paired, *, source='olga', condition_vj=True, resolve_genes=True, alpha_locus=None, beta_locus=None)[source]#

Add pgen_alpha, pgen_beta and pgen_paired to a paired single-cell frame.

Parameters:

paired (DataFrame) – A paired-chain frame (vdjtools.sc.pair.pair_chains() layout).
source (str) – Bundled model set — "olga" (OLGA-derived) or "learned" (native EM).
condition_vj (bool) – Condition each chain’s Pgen on its V/J call. False marginalises over all V/J unconditionally.
resolve_genes (bool) – Resolve a gene-level call (TRBV10-3) to a representative model allele (TRBV10-3*01) before scoring – see _gene_to_allele(). Default True, because CellRanger reports genes and without this every 10x row scores None. Set False to score only exact allele matches.
alpha_locus (str | None) – Locus of the α/light chain (e.g. "TRA", "IGK"); inferred from the alpha_v_call prefix if None.
beta_locus (str | None) – Locus of the β/heavy chain (e.g. "TRB", "IGH"); inferred from the beta_v_call prefix if None.

Returns:

paired with three added Float64 columns. pgen_alpha / pgen_beta are null for a cell missing that chain’s junction, or carrying a V/J call the model does not know; pgen_paired is null unless both are set.

Warns:

UserWarning – If every chain of a locus scored null – the usual cause is a V/J naming the model does not recognise, which would otherwise be an entire column of silent nulls.

Return type:

DataFrame

`vdjtools.sc.cluster_eval`#

Clustering-quality evaluation against antigen labels.

Clustering-evaluation metrics (purity / homogeneity / parsimony / q-measure).

Given a ground-truth labelling (e.g. the antigen/epitope a clonotype binds) and a predicted clustering (cluster ids), these functions score how well the clustering recovers the truth. They are the clustereval family of information-theoretic and set-overlap metrics, reimplemented here on a contingency matrix so they can grade any clonotype clustering — TCRnet components, metaclonotype groups, GLIPH motifs, …

The contingency matrix n[i, j] counts items with true class i and predicted cluster j (rows = true classes, columns = clusters), built with sklearn.metrics.cluster.contingency_matrix(). All entropies use natural logs.

Singleton convention. An unclustered item must get its own unique cluster id rather than being dropped or lumped into a shared “noise” cluster — otherwise purity and homogeneity are silently inflated. assign_singleton_ids() maps a sentinel (None / -1) onto distinct negative ids for exactly this.

vdjtools.sc.cluster_eval.assign_singleton_ids(pred, *, sentinel=None)[source]#

Give every unclustered item (sentinel) its own unique negative cluster id.

Parameters:

pred – Iterable of predicted cluster ids; sentinel marks unclustered items.
sentinel – The value flagging “no cluster” (default None; -1 is common).

Returns:

A list of cluster ids with each sentinel replaced by a distinct negative integer (-1, -2, …), so no two unclustered items share a cluster.

Example

>>> assign_singleton_ids([5, None, 5, None])
[5, -1, 5, -2]

vdjtools.sc.cluster_eval.purity(labels_true, labels_pred)[source]#

Cluster purity: mean over clusters of the dominant true class fraction.

purity = (1/N) · Σ_j max_i n[i, j] — 1.0 when every cluster is class-pure.

Return type:: float

vdjtools.sc.cluster_eval.normalized_purity(labels_true, labels_pred)[source]#

Purity rescaled against its one-cluster floor pmin = max_i n_i. / N.

(purity - pmin) / (1 - pmin); returns 1.0 when pmin == 1 (a single true class). This maps the trivial “everything in one cluster” purity to 0 and the perfect clustering to 1.

Return type:: float

vdjtools.sc.cluster_eval.inverse_purity(labels_true, labels_pred)[source]#

Inverse purity: mean over true classes of the dominant cluster fraction.

(1/N) · Σ_i max_j n[i, j] — the completeness-flavoured dual of purity.

Return type:: float

vdjtools.sc.cluster_eval.normalized_inverse_purity(labels_true, labels_pred)[source]#

Inverse purity rescaled against its all-singletons floor.

(inv - imin) / (1 - imin) with imin = |unique(true)| / N (the inverse purity you get when every item is its own cluster); returns 1.0 when imin == 1.

Return type:: float

vdjtools.sc.cluster_eval.homogeneity(labels_true, labels_pred)[source]#

Homogeneity: 1 - H(C|K) / H(C) — do clusters contain a single true class?

Returns 1.0 when H(C) == 0 (a single true class, trivially homogeneous). H(C|K) = max(0, H(C,K) - H(K)) with all entropies in natural logs.

Return type:: float

vdjtools.sc.cluster_eval.parsimony(labels_true, labels_pred)[source]#

Parsimony: 1 - H(K|C) / (ln N - H(C)) — penalises fragmenting a class.

Returns 1.0 when the denominator ln N - H(C) is 0. H(K|C) = max(0, H(C,K) - H(C)). Falls to 0 when every item is its own cluster (maximal fragmentation).

Return type:: float

vdjtools.sc.cluster_eval.q_measure(labels_true, labels_pred, beta=1.0)[source]#

Weighted harmonic mean of homogeneity h and parsimony p.

(1 + beta) · h · p / (beta · h + p); 0 when either h or p is ≤ 0 (there is nothing to balance). beta weights homogeneity relative to parsimony.

Parameters:: beta (float)
Return type:: float

vdjtools.sc.cluster_eval.cluster_eval(labels_true, labels_pred, *, beta=1.0)[source]#

Compute the full clustering-evaluation metric suite in one pass.

Parameters:

labels_true – Ground-truth class label per item (e.g. bound epitope).
labels_pred – Predicted cluster id per item (same length as labels_true). Unclustered items should already carry unique ids — see assign_singleton_ids().
beta (float) – Homogeneity weight for q_measure().

Returns:

Dict with keys purity, normalized_purity, inverse_purity, normalized_inverse_purity, homogeneity, parsimony, q_measure.

Raises:

ImportError – If scikit-learn is not installed (see the sc extra).
ValueError – If the label arrays are empty or unequal length.

Return type: