Recombination model workshop#
A model in vdjtools.model is three things: a manifest declaring the recombination Bayes
net, a set of tidy polars marginal tables (one per event), and the germline frames those
tables are keyed against. Because the probabilities are ordinary DataFrames, everything below is
table-in / table-out — nothing is hidden behind a binary format.
This page covers building a model from your own reference, fitting it to your own sequences, checking it, comparing two of them, and asking how much diversity one actually describes. For generation probability and sampling see User guide; for the complete surface see API reference.
Custom germline libraries#
from_germline() builds a model scaffold from any V(D)J library — your
own FASTA, a population-specific reference, or arda’s. from_arda is a thin wrapper over it.
from vdjtools.model import read_germline_fasta, validate_germline
from vdjtools.model.io import from_germline
germline = read_germline_fasta("V.fasta", "J.fasta", "D.fasta") # D optional -> VJ vs VDJ
validate_germline(germline) # tidy issue frame; empty means clean
template = from_germline(germline, locus="TRB", organism="human")
The frame needs allele, segment (V/D/J) and sequence; gene,
functional, cdr3_anchor and full_germline are filled in if absent. sequence is the
CDR3-region germline — for V from the conserved Cys104 codon to the 3’ end, for J from the 5’
end through the [FW]118 codon — or the full germline for D. If your FASTAs are full-length, pass
anchors= (OLGA’s *_gene_CDR3_anchors.csv format) and they are sliced for you.
Note
The single most damaging mistake in a custom library is a misplaced CDR3 anchor: it shifts
every deletion profile by a constant, and nothing downstream complains. That is why
validate_germline() checks that each V starts on a Cys codon and
each J ends on Phe/Trp, and warns when one does not.
The marginals of a fresh template are placeholders. Their support ranges bound what EM can later
learn, which is why ins_max is a parameter: a model whose insertion support stops at 40 can
never learn a 45-nt N-region.
Learning from your own sequences#
infer_frame() takes a clonotype frame, finds the junction column, and
builds the per-read V/J masks that make EM tractable on a D-bearing locus.
from vdjtools.model.infer import infer_frame, training_frame
model, report = infer_frame(template, clones, max_iter=15, gene_prior=1.0)
training_frame(model) # run, iter, loglik, n_scoreable, rel_change
Train on non-functional reads — out-of-frame or stop-codon, since both escaped selection. Keeping only the out-of-frame half conditions the training set on junction length modulo 3, which the insertion-length model would then happily learn.
gene_prior is a Dirichlet pseudocount over the germline’s functional alleles. 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 without the prior one unlucky iteration deletes a real gene for good.
Fine-tuning is the same call with a warm start:
tuned, report = infer_frame(model, new_clones, init="template", max_iter=5)
Watching, and surviving, a long fit#
EM on a D-bearing locus runs for a long time — IGH enumerates roughly 1,225 D pairs per read against TRB’s 9 — and the training log only becomes readable once the fit returns, so a slow run and a stuck one look identical. Two things fix that:
from vdjtools.model.infer import infer_native, print_progress, resume
model, rep = infer_native(template, seqs, progress=print_progress(), checkpoint="ckpt/IGH")
# ...interrupted...
model, rep = resume("ckpt/IGH", seqs, max_iter=10)
progress reports the log-likelihood and its relative change after every iteration — the exact
quantity compared against tol, so you can watch it approach. checkpoint saves the model
after each iteration (swapped into place, so a kill mid-write leaves the previous checkpoint
loadable), and resume() picks it up.
Resuming is exact: three iterations plus a resumed four reach the same log-likelihood, and the
same tables, as an uninterrupted seven. The checkpoint carries its training log and the resumed run
appends to it, so training_frame spans every attempt.
From the CLI: vdjtools model learn ... -v --checkpoint DIR, then
vdjtools model learn ... --resume DIR.
Each fit appends a run to model.training["runs"], which
save_model() writes beside the model as training.json. Models that were
never fitted here — every bundled one — simply have training is None.
Checking a model#
check_model() returns a tidy issue frame rather than raising, so every
problem in a model is visible at once and the result sorts, filters and writes like any other table.
from vdjtools.model import check_model
issues = check_model(model)
issues.filter(pl.col("severity") == "error") # empty == sound
severity is the contract: error means the model will crash or score wrongly, warn means
suspicious but usable, info is a note. The checks exist because each one has produced a silently
wrong answer at some point — a functional gene left at P = 0 makes Pgen exactly zero for every
clonotype using it; deletion mass past a germline’s length is unreachable and quietly lost; an
allele present in a marginal but not in the germline crashes the native packer far from the cause.
Unreachable deletion mass is reported as a fraction per allele, ranked, rather than one row per cell: a shared deletion-bin grid across alleles of different lengths always strands a little mass on the short ones, so what matters is how much.
Note
On a model imported from OLGA this particular finding is inherited, and is reported at
warn rather than error. OLGA’s own marginals carry the same mass — IGHV4-30-4*01
has Pgen identically zero in OLGA too — and vdjtools reproduces those arrays bit-faithfully
on purpose, so correcting it would break the exact-OLGA-Pgen invariant. The gene-collapsed
models (the default) are clean.
From the command line, vdjtools model check exits 1 on any error-severity issue, which makes it
usable as a build gate.
Information content and diversity#
Two different questions, reported side by side because they differ by orders of magnitude:
from vdjtools.model.analyze import total_entropy
from vdjtools.model.score import diversity
total_entropy(model) # per-event contribution, in bits
diversity(model, n=5000) # both entropies and both Hill numbers
Scenario entropy is the information in one recombination event, summed over the Bayes net. It
is an upper bound on the sequence entropy, because different scenarios can produce the same
junction. Sequence entropy is the entropy of the junction distribution itself, estimated by
Monte Carlo: sequences drawn from the model are distributed as Pgen, so E[-log2 Pgen] over
generated sequences is an unbiased estimator and the standard error comes free.
From those follow two Hill numbers — 2^H (the usual “~10^x distinct sequences” figure) and
1 / E[Pgen] (how many draws before two coincide, exact because Σ Pgen² = E[Pgen]). Human
TRB comes out at roughly 52 bits per scenario, 45 bits per sequence, and ~3·10¹³ effective
sequences.
Comparing two models#
from vdjtools.model.analyze import compare_models, compare_net_dot, compare_usage, render_dot
compare_models(olga, learned, by="gene") # per-event tv, tv_max, jsd_bits, support diffs
compare_usage(olga, learned, "v") # gene usage side by side
render_dot(compare_net_dot(olga, learned), "diff.pdf", fmt="pdf")
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 vanishing. Conditioned events
are averaged over parent groups weighted by the parent’s marginal, so a rarely-used V’s deletion
profile cannot dominate; tv_max reports the worst single group, which is what finds the one
broken gene an average hides.
Jensen-Shannon is the headline metric: symmetric, bounded by one bit, and finite when the supports differ — which is exactly the case here. KL is deliberately not reported, being infinite whenever one model assigns zero to something the other does not.
Use by="gene" to compare models built on different germline vintages or sources; an
OLGA-namespace model and an arda-namespace one only line up at gene level.
Likelihood and BIC#
from vdjtools.model.score import compare_pgen, model_fit, pgen_summary
model_fit(model, held_out_junctions) # loglik, k, AIC, BIC
pgen_summary(compare_pgen(a, b, sequences)) # how differently two models score one set
Two conventions matter. Likelihoods use nucleotide Pgen, because Σ Pgen_nt over all nt CDR3s
is 1, so log Pgen_nt is a proper log-likelihood. Pgen_aa sums only the in-frame, stop-free
nucleotide fiber of a translation, so an amino-acid log-likelihood is unnormalized and its missing
constant differs between models — it is a relative score on one fixed sequence set, never an
absolute one. And a sequence the model cannot generate gets pgen = 0 with a null
log-probability, never -inf; aggregates run over the scoreable subset and always report
n_scoreable beside n, so a flattering log-likelihood earned on 10% of the data is visible.
The free-parameter count behind BIC is support-based: per normalization group,
max(occupied cells - 1, 0), dropping groups that are undefined (all-zero) or unreachable (their
parent has zero marginal). Counting rows instead would put human TRB’s v_3_del at ~3,600
parameters when ~700 are real. It is a lower bound — a support count cannot distinguish a structural
zero from a parameter EM drove to zero — so BIC is comparable only between models counted the same
way, which is the case for any two compared through this function.
Warning
A log-likelihood computed on the sequences a model was trained on is that model’s own EM objective, which EM increases by construction. It validates nothing. Score a held-out set.
Known quirks of the OLGA models#
The bundled olga model set is a bit-faithful import of OLGA’s published models: vdjtools
reproduces their marginals exactly, and native Pgen matches olga’s own to machine precision
across all seven loci. Faithful means faithful to the defects too. These are properties of the
source models, verified against OLGA’s raw model_marginals.txt with OLGA’s own parser — not
import bugs, and not things to “correct” here, because doing so would break the exact-Pgen
invariant that makes the import checkable in the first place.
Deletion mass on trims that cannot be reached#
OLGA stores P(delV | V) on one deletion-bin grid per locus, sized for the longest allele. An
allele whose CDR3-region germline is shorter than that grid can carry probability on trims longer
than it has nucleotides. The Pgen DP never visits those, so the probability is not redistributed —
it is simply absent from every Pgen through that allele.
Measured on the shipped models (fraction of each allele’s own deletion distribution that is unreachable; identical in OLGA’s arrays and in ours):
allele |
unreachable |
consequence |
|---|---|---|
|
100 % |
|
|
80.9 % |
Pgen through this J is roughly 5x too low |
|
54.7 % |
Pgen through this V is roughly 2x too low |
|
52.0 % |
|
|
14.0 % |
|
|
10.6 % |
Across the whole set this affects a minority of alleles per locus (IGH V is the worst: 60 of 62
alleles carry some, most of it small). check_model() reports it as
deletion_unreachable at warn severity for an OLGA-sourced model, with the fraction, so you
can see whether a gene you care about is affected.
What to do about it. Use the gene-collapsed models — the default — where the effect is absorbed
by the collapse, or use the learned set, which is refit from data on arda germline and does not
inherit the grid. If you need allele resolution and an affected gene, be aware Pgen through it is
an underestimate.
Genomically impossible D–J pairs#
V(D)J recombination is deletional — it excises the DNA between the segments it joins — so a D can only be joined to a J lying 3’ of it. Almost every locus puts all D genes 5’ of all J genes, so this constrains nothing. TRB interleaves:
5' — TRBV… — TRBD1 — TRBJ1 cluster — TRBD2 — TRBJ2 cluster — TRBC2 — 3'
so a TRBD2–TRBJ1 join is physically impossible. OLGA’s model does not encode that: it gives
P(TRBD2*01 | TRBJ1-6*01) = 0.333, a third of that J’s D mass on a pair that cannot occur. The
bundled olga set reproduces it faithfully and check_model() reports
it as impossible_dj_pair at warn.
For models vdjtools fits, this is a hard genomic fact rather than a parameter, and the mask is
applied inside the M-step, so EM cannot relearn the pair on the next iteration. It needed to be:
a D match is only 10–18 nt — weak evidence on short reads — and an earlier build of the human TRB
model had reached P(TRBD2 | TRBJ1-1) = 0.091 before the constraint existed. Enforcing it also
improved the fit, converging in 9 iterations instead of 11 to a slightly better log-likelihood,
because mass that had been leaking onto impossible pairs went to the D that actually generated the
read.
enforce_dj_order() repairs a model fitted before the constraint existed.
Genes with no CDR3-region germline#
OLGA leaves the CDR3-region germline empty for a handful of ORF alleles while still giving them
usage. They can be drawn by the generator but score Pgen = 0, and human TRBV23-1 — 8.6 % of
a real TRB repertoire — is one of them. check_model reports these as unscoreable_gene_mass.
from_olga() accepts derive_orf=True to reconstruct the missing germline
from the full germline and the anchor; the learned models are built that way, while the olga
set keeps it off so it stays an exact Pgen oracle.
V/J usage is protocol-specific#
Not a defect, but the most common way to get a wrong answer with these models. OLGA’s were fit to
DNA-multiplex data; the learned set to 5’RACE reads. The two amplify different V genes at very
different rates — TRBV19 is 3.1 % of OLGA’s usage and 37 % of these 5’RACE reads — so neither
marginal is right for your library. The junction model (trims, insertions, dinucleotides) is the
shared, transferable part. Use rescale_usage() before scoring.
Out-of-frame input#
olga’s own compute_nt_CDR3_pgen rejects an out-of-frame junction outright (“Invalid
nucleotide CDR3 sequence”). vdjtools scores it, because the generation model is defined before
selection and out-of-frame rearrangements are exactly what it is trained on. This is a deliberate
difference in input validation, not in the probability: on any sequence OLGA accepts, the two agree
exactly.
Extending the allele library#
from vdjtools.model.infer import extend_alleles
bigger = extend_alleles(model, load_germline("TRB", "human"))
bigger, report = infer_frame(bigger, clones, init="template", max_iter=5)
A new allele of a gene the model already has is seeded from a gene-mate; a brand-new gene is seeded from the germline-nearest existing allele at a floor mass, because 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 unreachable mass.
Each pre-existing gene keeps its total usage: alleles of one gene are alternative versions of the same gene — a diploid carries at most two — so a richer library must split a gene’s mass more finely, never multiply it. Existing alleles are never modified, including their germline.
This seeds; it does not estimate. Follow it with a warm-start fit.
Rescaling V/J usage#
from vdjtools.model import rescale_usage
scored_with = rescale_usage(load_bundled("TRB", "learned"), my_sample)
V and J usage is protocol-dependent and the junction model is not: 5’RACE and DNA-multiplex
amplify different V genes at very different rates, so neither usage marginal is right in general.
Learn the junction model once, then set P(V) from the repertoire you are actually about to
score. Pass an out-of-frame sample only if that is what you will score — a pseudogene’s
rearrangements are never productive and are therefore enriched out of frame.
Exporting and importing tables#
from vdjtools.model.io import load_model, marginals_frame, save_model, set_marginals
marginals_frame(model) # every marginal as one long frame
set_marginals(model, edited_frame) # ...and back again
save_model(model, "mymodel", fmt="tsv") # a hand-editable model directory
load_model("mymodel") # format auto-detected, dtypes restored
Building from the read corpus#
The models shipped in the wheel are fitted to real 5’RACE reads. That whole pipeline —
fetch FASTQ, map with arda, collapse to unique clonotypes, run EM — is
build_all(), parallel across chains:
vdjtools model build --chains TRB,TRA,IGH --workers 4 -o models/
It needs HuggingFace access to the source dataset and arda’s mmseqs2.
For examples and tests, two arda-mapped clonotype sets (human TRA and TRB, out-of-frame) ship in the source tree as gzipped FASTA and load with no network, no arda and no mmseqs2:
from vdjtools.model.data import load_prepared
clones = load_prepared("human", "TRB") # junction, v_call, j_call, d_call, d2_call, count
model, report = infer_frame("TRB", clones, max_iter=10)
The V/J/D calls ride in the FASTA header (>{i}|{v_call}|{j_call}|{d_call}|{d2_call}|{count})
because EM needs them for its per-read masks. These files live under tests/ and are not
packaged into the wheel; from an installed vdjtools use
prepare() or pass path=.
Note
Real annotated junctions occasionally carry an ambiguous base, and the recombination model is
defined over A/C/G/T only. Both training entry points substitute A by default and warn with
the count (ambiguous=None drops the clonotype instead). Substituting keeps the clonotype —
one uncertain position in a junction that is otherwise good evidence — and on these reads it
affects ~0.01% of rows. It is a substitution, not a marginalization, so data with many ambiguous
positions should use ambiguous=None.
Command line#
Every operation above has a vdjtools model counterpart. A model is named either as a directory
or as LOCUS[:source[:organism]]:
vdjtools model check TRB:learned # exits 1 on any error
vdjtools model template --locus TRB -o template/
vdjtools model learn clones.tsv -t template/ -o fitted/
vdjtools model log fitted/ # loglik per iteration
vdjtools model entropy TRB:olga --table total
vdjtools model diversity TRB:olga -n 5000
vdjtools model compare TRB:olga TRB:learned --by gene --dot diff.pdf
vdjtools model compare-pgen TRB:olga TRB:learned seqs.tsv --summary
vdjtools model loglik seqs.tsv TRB:learned # loglik, k, AIC, BIC
vdjtools model extend fitted/ --locus TRB -o extended/
vdjtools model rescale TRB:learned my_sample.tsv -o rescaled/
vdjtools model export TRB:olga --format tsv -o trb_tsv/
vdjtools model net TRB:olga --format pdf -o bn.pdf