Selection evidence as information#
Scope and status.
Reference profiles describe which receptor features distinguish a target from a declared background. The gene–length decomposition and V-conditioned motif scores are implemented diagnostics. A tested class-based selector regressed against the accepted original baseline and was not promoted. The equations below explain evidence and calibration; they do not establish an optimal scorer or change the annotation default.
Reading a reference profile#
Keep the following axes separate:
V support measures gene-usage enrichment. Conditional motif support asks whether a junction pattern adds information within V; it does not isolate a J-gene effect.
Length KL/JS describes the entire spectratype. JS remains finite when an unsmoothed KL is infinite because the background lacks a reference length.
Same-size neighbour comparisons distinguish concentrated local families from a large reference. Mean degree and the fraction of connected vertices answer different questions. Self-inclusive Hill2 re-expresses mean degree.
Background publicity measures sequence availability in the declared control, not sharing across donors or experimental specificity.
Occurrence rarefaction estimates represented mass and richness at a stated clonotype/community resolution. Graph acquisition interpolates observed connectivity; it does not estimate unseen receptors.
Alpha/beta support retains each marginal and observed pair links. A joint-class row grouping single-chain metrics is not a paired-score performance estimate.
full20 denotes the full amino-acid alphabet: core exact triplets and all
single-wildcard four-mer masks share one selection-aware shuffle maximum.
It does not denote untrimmed junctions. The current profile default is 999
shuffles; historical 199-shuffle index outputs retain their own identities.
At 999 shuffles, the inclusive 0.005 tail permits up to four shuffled exceedances.
For class/scorer heatmaps, prefer MCC and balanced accuracy to raw accuracy under class imbalance. State the number of defined epitopes for each macro mean, keep undefined MCCs unavailable, and compare scorers on matching denominators. Reference-label retrieval, experimental binding discrimination and final-event null calibration are different endpoints. A descriptive class does not prescribe a scorer, a molecular mechanism or a binding probability.
Likelihood ratios, selection factors and bits#
Fix species, locus and full peptide–MHC restriction, denoted by context \(c\); write \(e\) for the target within that context and suppress \(c\) in the equations. Let \(x=(L,V,J,Z)\) describe an anchor-inclusive junction, its length \(L\), gene calls \(V,J\), and remaining sequence features \(Z\). Define \(P_e\) as the distribution in the observed target repertoire and \(P_0\) as the declared background distribution, using the same sampling unit on both sides. For example, distinct junction/gene keys on both sides measure repertoire composition; reference keys divided by control event frequencies instead compare two explicitly different sampling schemes. The empirical selection factor and pointwise information are
A positive \(i_e(x)\) denotes enrichment, a negative value depletion, and one bit doubles the likelihood ratio. For known distributions, \(Q_e\) is the likelihood ratio for the simple hypotheses \(P_e\) versus \(P_0\). If the target prior in an explicitly defined two-population mixture is \(\pi_e\), Bayes’ rule gives
Without that sampling model and prior, neither an empirical ratio nor a calibrated tail is a binding posterior. Observed repertoire enrichment includes biological selection, sampling and ascertainment; it does not by itself isolate an antigen-specific causal selection factor. In particular, a selected AIRR control is not the recombination distribution \(P_{\rm gen}\).
The mean information per sampled receptor is the Kullback–Leibler divergence,
where \(H(P_e)=-\sum_xP_e(x)\log_2P_e(x)\) and \(H(P_e,P_0)=-\sum_xP_e(x)\log_2P_0(x)\). All these quantities are in bits per receptor under the declared sampling unit. A low target entropy alone does not show antigen-specific selection: a concentrated feature can be equally concentrated in the control. These identities assume \(P_e\) has no positive mass where \(P_0\) is zero. A zero empirical control count is unresolved support, not evidence that the population probability is exactly zero. Report the affected counts and leave the comparison unresolved; silently dropping those categories, renormalising them away or treating them as infinite evidence changes the question.
Sampling likelihood and mutual information.
For \(n\) independent observations categorised by a fixed feature map, let \(n_y\) be the category counts and \(\widehat p_y=n_y/n\). With a known background probability \(p_{0y}\), the multinomial log-likelihood ratio of the fitted categorical distribution to that background is
The multinomial coefficient cancels; zero-count terms contribute zero. This connects descriptive information to statistical likelihood, but is not a P-value: the alternative was fitted to those counts, the background is usually estimated, and donor/expansion dependence can violate independent sampling. For a binary population label \(U\) with target probability \(\pi\), set \(M=\pi P_e+(1-\pi)P_0\). Then mutual information is
It measures label information under that explicit mixture and differs from \(D_{\rm KL}(P_e\Vert P_0)\). Neither quantity supplies an automatic performance-optimal weighting of overlapping features.
Which feature carries the information?#
The probability chain rule gives an exact decomposition, without feature independence:
Taking the target expectation decomposes \(I_e\) into length information and target-weighted conditional divergences for V, J and the sequence remainder. Each expected term is nonnegative, although individual query contributions can be negative. The allocation depends on the conditioning order; it must be reported, not interpreted as a unique causal partition. A V-gene contribution is constant across queries with the same V and conditioning variables: it cannot resolve their junction-specific recognition.
More generally, for a deterministic feature map \(Y=\phi(x)\),
This follows by factoring each probability into its feature marginal and conditional remainder. Thus a projection retains some of the discriminating information and loses the nonnegative conditional remainder. It motivates comparing V/J, length, motif and neighbourhood projections against the same control. It does not justify using their marginal divergences as additive score weights: correlated features can encode the same information. Finite-sample plug-in divergences are descriptive estimates, not unbiased information estimates or proofs that a rule will generalise.
For a genuine pair \((a,b)\), define \(C_j(a,b)=P_j(a,b)/[P_j(a)P_j(b)]\) for \(j\in\{e,0\}\). Then
This separates marginal chain evidence from excess pairing dependence. Independent random TRA/TRB pairing sets \(C_0=1\) in the specified null; it does not assert independence in natural paired repertoires. When one chain is absent, use its observed partner’s marginal evidence; an unobserved alpha-chain germline signal cannot be assigned to a beta-only query.
Masked four-mers as an information-retrieval score#
The implemented motif-evidence diagnostic uses the MF–IKF recipe; its equations below specify the calculation without requiring that precursor code. Map amino acids to the eleven classes AILV, RHK, CM, ST, DE, NQ, G, F, P, W and Y. Slide a four-residue window along the full junction and replace each of its four positions in turn by a wildcard. Let \(\mathcal T(x)\) be the resulting set of distinct tokens; repeated occurrences of one token in a junction count once. V alleles are pooled to gene identity; the reduction acts on amino acids, not on V-gene names.
For token \(t\), let \(x_{lve}(t)\) count distinct reference junctions with length \(l\), V gene \(v\) and full-pMHC label \(e\) containing \(t\). Write \(x_{ve}=\sum_lx_{lve}\), \(x_e=\sum_{l,v}x_{lve}\) and \(x(t)=\sum_e x_e(t)\) over the declared labelled reference universe of \(X\) keys. A key includes its target label, so a receptor with several labels contributes to several labelled keys. The fixed support set is
For ordinary annotation the score is
The implementation takes the maximum supported length-stratum contribution per token before summing; the displayed unfiltered weight is identical across those strata. MF is the target/V share of token-bearing labelled keys, not the usual within-target term frequency \(x_{ve}/N_{ve}\). IKF is the surprisal of token presence under the labelled reference-key distribution, in nats; dividing by \(\ln2\) expresses it in bits. This is a TF–IDF analogue, not a fitted likelihood ratio against AIRR control. Rare tokens receive more weight only when supported and concentrated in the target/V combination. Overlapping windows are dependent; their sum is a retrieval score, not a product of independent probabilities. Its support cutoffs are inherited implementation choices, not biophysical constants or a newly fitted parameter grid.
In benchmark mode only exact-hit contributions are filtered from the numerator: replace \(x_{V_qe}(t)\) there by the remaining hit count. The support set, \(x(t)\), \(X\), priors, graph and all fitted weights remain fixed. No-hit token contributions vanish; no reference model is refitted for an individual query. The motif score can expose a useful sequence axis, but selecting it for an epitope still requires evidence that this axis differs from the declared control.
Control-relative information in a motif projection#
For a fixed masked token \(t\), define the binary feature \(Z_t(x)=\mathbf1\{t\in\mathcal T(x)\}\). Within a V-gene/length stratum \((v,l)\), let \(n,N\) be reference keys with the motif and all reference keys, and \(m,M\) the corresponding control counts. Each population uses distinct full-junction/V-gene keys; repeated windows, V alleles and alternative J annotations do not multiply a key. With \(p=n/N\) and \(q=m/M\), the presence contribution and binary information are
The first can be negative for a depleted motif; the second is the Kullback–Leibler divergence of one Bernoulli projection, in bits per conditional junction/V key. Zero target-probability terms contribute zero by continuity. An empty control stratum, unobserved control motif with positive target mass, or unobserved control absence with positive target absence leaves the corresponding empirical comparison unavailable; the counts remain explicit. The all-length scope conditions on V only, while the exact-length scope also controls the opportunity to contain more windows in a longer junction. J usage and motif position are not additional conditioning variables in this diagnostic.
The shipping reference-structure --control command writes these reference-observed token projections to .motifs.tsv. It reuses the existing tokenizer in bounded blocks and reduces once-per-key presence counts, with one control model per species/locus and one fixed model per target. Tables are written one target at a time without retaining every target’s motif rows. Overlapping tokens are dependent: their binary divergences are not an additive estimate of joint sequence information. These descriptive projections help identify informative families and motifs; the retained MF–IKF retrieval score, final rule selection and complete-score calibration remain distinct calculations.
--motif-mode selects exact3 (full-AA 3-mers), wild4 (full-AA
4-mers with one wildcard in each possible position), imgt3 or
imgt_wild4 (default). The latter two reduce residues to the eleven IMGT
classes before counting. motif_mode and motif_alphabet identify the
projection in every output row. The same binary-information formulas apply
to all four views; these values are distinct from TF-IDF retrieval weights.
A depleted germline category can still contain a strongly selected junction solution. Let \(N_{\rm all},M_{\rm all}\) count all reference/control keys in the same declared populations. The joint presence evidence factors exactly as
Only matched population units and the same stratum permit this addition. A
pooled marginal V factor cannot replace the V/length term in an exact-length
comparison. A positive conditional term can outweigh stratum depletion;
neither term alone is a receptor veto. This is one joint feature projection,
not a sum over overlapping motifs or a new score-selection rule. Unobserved
control support remains unavailable instead of receiving an invented floor.
The same .motifs.tsv export appends total_reference_keys,
total_control_keys, stratum_bits and joint_presence_bits. The first
two are whole-population totals, not within-stratum counts; the last two are
the stratum term and direct joint-count term in the identity above. Existing
conditional presence bits, binary KL and availability statuses remain separate.
Germline-loop information with allele ambiguity#
reference-structure --control also writes .germline-map.tsv, .germline.tsv and .germline-information.tsv. These project the existing V-gene counts onto individual aligned CDR1, CDR2 and CDR2.5 columns, using full amino acids and the eleven IMGT classes separately. Column numbers are one-based columns in the bundled alignment, not structural residue numbers. A gene has a confirmed residue only where every bundled allele agrees. Missing genes, missing loop templates, stop/unknown symbols and allele-dependent columns remain explicit; an omitted allele never selects *01. A dot denotes aligned absence within a populated loop, while an entirely empty template is unavailable.
The sampling unit remains a distinct full-junction/V-gene key. Keys from two genes retain their separate mass after those genes map to the same residue, even if their junctions agree. Let \(w_e(g)\) and \(w_0(g)\) count such keys for gene \(g\), and let \(N=\sum_gw_e(g)\) and \(M=\sum_gw_0(g)\). For one loop column and alphabet, let \(\mathcal A_g\) be the set of residues across the modeled alleles; missing models or template symbols admit the entire alphabet. For residue \(a\), define
For any allele assignments within those allowed sets, the actual residue fractions obey \(p_e(a)\ge n_e(a)/N\) and \(p_0(a)\le u_0(a)/M\). Thus, when the displayed counts are positive,
residue_bits_lower is this lower bound: every compatible or unresolved control gene can contribute to its denominator. It is an empirical bound over allele assignments under the declared template model, not a sampling confidence interval, population guarantee or P-value. Zero compatible control counts leave it unavailable instead of assigning infinite evidence. The exact-length scope uses the corresponding length-stratum counts.
consensus_bits and consensus_information_bits instead describe the deterministic gene-consensus annotation, with unresolved reference mass left unavailable. Where that expectation is defined, the projection cannot carry more information than the original V-gene distribution on the same key population. This annotation information differs from the unknown allele-resolved residue information; the lower-bound field handles compatible background ambiguity explicitly. Correlated residues, loops and reduced alphabets are separate views and must not be summed as independent evidence. None of these diagnostics selects a score automatically.
Geometry as a projection of selection#
For a fixed dissimilarity \(d\) and radius \(r\), define \(B_d(q,r)=\{x:d(q,x)\le r\}\). Local enrichment is another likelihood ratio, now for the binary neighbourhood event:
For a nonnegative kernel \(K\), the corresponding ratio is
when the denominator is positive. These are population definitions of candidate geometry evidence, not identities with every shipped reducer. Distances, radii, kernels and sampling units must match on both sides. A smooth substitution score can retain within-neighbourhood variation; a discrete motif or edit graph can expose separated recognition families. Reduced alphabets can collapse distinct junctions to zero distance, and an alignment dissimilarity need not define a Euclidean manifold.
For \(n\) reference keys, the empirical pair-mass curve and its scale-dependent slope are
The slope describes local growth only where the curve has an estimable nonzero range; discrete plateaux, small samples and disconnected families do not yield a single intrinsic dimension. Compare target and control curves under the same geometry and units before attributing density to antigen selection. Nearest-neighbour distances \(u_i=\min_{j\ne i}d(x_i,x_j)\) give a complementary sorted curve. A Kneedle diagnostic locates a change in its rescaled shape; its knee is a geometric scale proposal, not a significance threshold or a confidence probability. A constant curve has no identified knee. The existing nearest-distance diagnostic and the pair-mass slope are distinct summaries; the latter definition does not claim that it has already been computed.
At a declared radius, graph components with masses \(p_k=n_k/n\) have entropy \(H_C=-\sum_kp_k\log_2p_k\). This describes fragmentation, not antigen information without a matched control comparison. The cross-length nearest-neighbour fraction at query length \(l\) is
Strict ties are not cross-length wins; an empty denominator is unavailable. Unlike \(1-\sum_lP_e(L=l)^2\), this statistic measures an actual geometric relation, not length diversity alone. Neither statistic is, by itself, a validated gap-score weight.
Human-readable rules and comparable tails#
A germline-anchored class has informative V/loop usage relative to control, a motif-constrained class has informative conditional junction motifs, and a distributed-neighbourhood class has informative smooth local enrichment. Multiple recognition families can coexist within a target; these labels describe measured axes rather than force each peptide into one exclusive class. MHC class, allele, peptide chemistry and available structural contacts motivate an interpretation of those axes; peptide composition alone does not determine receptor recognition or its score. Species, chain coverage, length and family support remain explicit conditions. No informative supported projection means use the retained general score and report the available evidence, rather than invent specificity.
Let \(A\) denote the complete procedure, including data-dependent rule choice, and \(T_e(q,R)=T_{A,e}(q,R)\) its higher-is-better score. The two directional tails are
The reference direction uses composition-matched mock references; the query direction uses control junctions matched to query Pgen strata. A stratum conditions on generation probability, not on exact sequence identity or all biological selection effects. Each mock reference gets its own fixed model and the same rule-selection algorithm. The real reference model remains fixed across queries; the query direction applies the same procedure to every random query. Ties enter the upper tail, all null no-hits remain in the denominator, and observed no-hits receive one. Under the respective conditional exchangeability assumptions, these are conservative Monte Carlo rank P-values. Matching Pgen alone does not prove those assumptions.
If both directional P-values are valid for nulls \(H_R,H_Q\), then \(p_{{\rm both},e}=\max(p_{R,e},p_{Q,e})\) is valid for \(H_R\cup H_Q\). Indeed, under \(H_R\), \(\Pr(p_{{\rm both},e}\le\alpha)\le\Pr(p_{R,e}\le\alpha)\le\alpha\), and the same holds under \(H_Q\). No independence assumption is needed. This gives a common tail scale for different epitope scoring rules, conditional on their own correctly reproduced nulls; it does not make the scores themselves universal likelihoods. Selecting the smallest tail across labels or the strongest paired/marginal axis requires calibration of that complete selection event, not reuse of a component tail. The same principle applies when reporting any confident specificity in a repertoire: the reported query/label search scope must match the calibrated event.
With \(B=199\), a plus-one tail cannot be smaller than \(1/200=0.005\). Many distinct strong raw scores therefore share the same empirical tail. Preserve raw within-epitope ranking alongside cross-epitope tails; finite calibration resolution must not silently erase ranking information. The present formulation supplies exact identities and conditional validity statements, not a guarantee of future ranking accuracy, biological binding or state-of-the-art performance.
Implementation entry points#
Use Command-line interface for motif-evidence and reference-structure --control. The latter produces pointwise selection bits in .genes.tsv and conditional/marginal information in .information.tsv.
See Score definitions for the current default reducers and Two directions of control calibration
for the two control directions. The information-directed selector described here
is a mathematical development objective, not an additional CLI mode.
Database occurrence coverage#
reference-coverage reuses vdjtools’ iNEXT abundance estimators and
augmented-assemblage bootstrap. Let \(n_i\) be the VDJdb occurrence count
of taxon \(i\), \(n=\sum_i n_i\), \(S\) the number of observed taxa,
and \(f_1,f_2\) the numbers with abundance one and two. Taxa are either exact
full-junction/V/J-call keys or Leiden communities; a community sums its members’
counts. These are database occurrences, not cells or reads.
The finite-sample Chao1 estimate is
For \(f_1>0\), set \(A=nf_0/(nf_0+f_1)\). Estimated coverage at observed effort is \(\widehat C=1-(f_1/n)A\); with no singletons it is one under this estimator. The zero-doubleton branch is explicitly labelled. In particular, a single observed taxon with one occurrence has no information about unseen taxa: the estimator’s boundary value is not a sufficiency finding.
Unseen probability mass is \(1-\widehat C\). Estimated unseen-type fraction is \((\widehat S-S)/\widehat S\); these quantities answer different questions. vdjtools provides interpolation/extrapolation through twice the observed effort, with normal 95% bootstrap intervals for richness. Coverage receives a 2.5/97.5 percentile envelope from the same owning augmented-assemblage probability law, evaluated through vdjtools’ bounded native rarefaction batches. Bootstrap replicates and seeds are explicit. vdjtools samples a large equal-probability unseen pool through an aggregated unseen category and uniform virtual taxon IDs. This preserves the augmented-assemblage law with memory proportional to observations rather than estimated unseen richness. The adapter bounds actual bootstrap working memory and requires total occurrences below \(2^{31}\).
The formulas follow the iNEXT abundance implementation. The graph partition is held fixed within occurrence rarefaction; its conditional plateau does not establish exhaustion of unseen communities. The community adaptation is distinct from distance-weighted functional diversity. This command does not implement donor/study incidence Chao2 or a numeric adequacy threshold. The requested source-incidence route remains a separate sampling model.