Which model do I need?

A BEAM model describes whatever set of cells it is built from. Choosing that set is the first decision of any analysis, and it decides which questions the model can answer.

One repertoire

R1

↓ build

model of R1

Describes one sample or one person. Compare samples, follow a person over time, plan the depth for this sample.

One group

R1R2R3

↓ pool the cells, then build

model of the group

Describes a cell type, a condition or a cohort. Compare groups, ask which group a sequence resembles, predict sharing within the group.

Everything: a reference

…………

↓ pool all, then build

reference model

Generalises over all the data given, like the reference model built from 38,538 repertoires. The baseline for "how likely in people in general".

One model per repertoire

Built from one sample, a model describes that sample: one person, one tissue, one time point. This is the model to use when the unit of comparison is the sample:

• comparing two samples, or many (distances and clustering);
• following one person over time, or matching samples of the same person;
• planning the depth for this kind of sample.

A model of one sample describes its common sequences from their counts and its rare ones through the prior; the more cells it is built from, the more detail it holds.

One model per group

A group is any set of repertoires that belong together: a cell type (CD4, CD8, naive, memory), a tissue, a condition, a patient cohort. Pool the cells of all its repertoires and build one model. It describes the group as a whole, averaged over its members, and it is more detailed than any one member’s model because it has seen more cells.

levels.py

import numpy as np
import BEAM
from BEAM import BeamModel

reference = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")          # 38,538 repertoires: the population
# donors: three repertoires from one group (stand-in population B: longer CDR3s)
per_repertoire = [BeamModel(cells) for cells in donors]   # one model per repertoire
group = BeamModel([s for cells in donors for s in cells])  # the group's cells, pooled


def long_share(model):
    """P(length >= 17) under a model."""
    return sum(p for L, p in model.length_distribution.items() if L >= 17)


for name, model in [("donor 1", per_repertoire[0]), ("donor 2", per_repertoire[1]),
                    ("donor 3", per_repertoire[2]), ("group", group), ("reference", reference)]:
    print(f"{name:10} P(length >= 17) = {long_share(model):.3f}   cells counted: {model.total_weight:>15,.0f}")

Output

donor 1    P(length >= 17) = 0.430   cells counted:          30,000
donor 2    P(length >= 17) = 0.426   cells counted:          30,000
donor 3    P(length >= 17) = 0.426   cells counted:          30,000
group      P(length >= 17) = 0.427   cells counted:          90,000
reference  P(length >= 17) = 0.115   cells counted:   5,629,439,400

Full program, including the set-up hidden here

The stand-in group here has longer CDR3s than the population (stand-in population B). Each donor’s model shows it, the pooled group model shows it with more data behind it, and the reference model, built from everyone, does not.

Use a group model to:

• compare groups directly, for example CD8 against CD4 (Comparing groups);
• ask which group a sequence or repertoire resembles, by scoring it under each group’s model;
• predict sharing within the group, with the group model as the population.

Pooling the cells and building, or building per repertoire and adding the models up, give the same model:

pooled_equals_merged.py

import numpy as np
from BEAM import BeamModel
# cells_1, cells_2: two repertoires of one group
pooled = BeamModel(cells_1 + cells_2)                         # build from all the cells
merged = BeamModel(cells_1).merged(BeamModel(cells_2))        # or add up two built models

difference = np.abs(pooled.log_pseq_batch(probe) - merged.log_pseq_batch(probe)).max()
print(f"largest difference in log probability over 1,000 sequences: {difference:.1e}")

Output

largest difference in log probability over 1,000 sequences: 0.0e+00

Full program, including the set-up hidden here

So group models can be assembled from per-repertoire models later, without going back to the raw data.

One reference over everything

Pool every repertoire you have, of every group, and the model generalises over all of them. The reference model shipped with BEAM is this level for human TCRβ: 38,538 repertoires from 447 studies. It is the baseline for “how likely is this sequence in people in general”, and the population for sampling and publicness calculations when nothing more specific is known.

You can build your own reference the same way, for example over all repertoires of one species, one chain or one sequencing platform.

Leaving a cohort out

When you evaluate a cohort against a reference that contains it, the cohort’s own sequences look more probable than they would to an independent reference. without removes a cohort’s counts from a model; the result is the model built without the cohort. The example uses a stand-in corpus of 13 repertoires, three of which form the cohort; how much the cohort’s sequences gain depends on how large a share of the corpus the cohort is, and is far smaller in a reference of 38,538 repertoires:

leave_out.py

import numpy as np
from BEAM import BeamModel
# corpus_cells: the repertoires a reference is built from; cohort_cells: your cohort,
# which is part of that corpus
reference = BeamModel(corpus_cells + cohort_cells)
cohort = BeamModel(cohort_cells)

held_out = reference.without(cohort)              # the reference with the cohort's counts removed
rebuilt = BeamModel(corpus_cells)                 # the same thing, built from scratch

cohort_seqs = sorted(set(cohort_cells))
difference = np.abs(held_out.log_pseq_batch(cohort_seqs) - rebuilt.log_pseq_batch(cohort_seqs)).max()
print(f"largest difference from rebuilding without the cohort: {difference:.1e}")
gain = np.mean(reference.log_pseq_batch(cohort_seqs) - held_out.log_pseq_batch(cohort_seqs))
print(f"mean log P of the cohort's sequences is {gain:.2f} nats higher when the cohort is included")

Output

largest difference from rebuilding without the cohort: 0.0e+00
mean log P of the cohort's sequences is 0.75 nats higher when the cohort is included

Full program, including the set-up hidden here

With the reference model this is BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation").without(cohort), when your cohort is one of the studies it was built from; the paper evaluated its predictions for the Emerson cohort this way. The subtraction is only right if the cohort’s counts are exactly the counts the reference holds for it: the same filtered sequences (3 to 27 standard amino acids) and the same weights (template counts where the platform reported them, otherwise one per distinct sequence per repertoire). Build the cohort’s counts with these rules, or pass them directly with without(sequences, abundances=...).

Merging counts or mixing models

Models of several repertoires or groups can be combined in two ways, and the two describe different things.

Merging adds the members’ counts and builds one model from the total: the same model as building from all their cells at once (a.merged(b), above). Every row of each expert then holds the members’ counts together,

$$n(\text{row}, x) \;=\; \sum_m n_m(\text{row}, x)$$

Each row sees only a position and a few residues of context, so the merged model learns the members’ local patterns and recombines them: it can start a sequence the way one member does and end it the way another does.

A mixture keeps the members’ models whole. With weights $w_m$ that add up to one (equal by default),

$$P_{\text{mix}}(s) \;=\; \sum_m w_m\, P_m(s)$$

pick member $m$ with probability $w_m$, then draw a whole sequence from that member’s model. A sequence that no member writes stays unlikely.

Merged: add the counts

group ACASS······YEQYF

group BCSA·······TEAFF

↓ one table: every row holds n_A + n_B

writesCASS······YEQYFCSA······TEAFFCASS······TEAFFCSA······YEQYFdashed: recombined, in neither group

Each row sees only a position and a few residues, so the model joins one group's start to the other's end.

Mixture: keep the models

group ACASS······YEQYF

group BCSA·······TEAFF

↓ pick a group (½ each), then a whole sequence from it

writesCASS······YEQYFCSA······TEAFF

P(s) = ½ P_A(s) + ½ P_B(s): a sequence neither group writes stays unlikely.

Two groups that differ at both ends show the difference:

pool_or_mix.py

import numpy as np
from BEAM import BeamModel, Mixture
# group_a: 20,000 CDR3s that start CASS and end YEQYF
# group_b: 20,000 CDR3s that start CSA and end TEAFF
merged = group_a.merged(group_b)         # the counts added: one model
mixture = Mixture([group_a, group_b])    # pick a group (1/2 each), then a sequence from it

print("sequence           group A   group B   mixture    merged")
for s in ["CASSLGQGYEQYF", "CSARGNTEAFF", "CASSLGQGTEAFF"]:
    print(f"{s:16}  " + "  ".join(f"{m.pseq(s):8.1e}" for m in (group_a, group_b, mixture, merged)))

# How much of the merged model's probability falls on sequences that start like one group
# and end like the other, summed over every length
def share(model, start, end):
    total = 0.0
    for L, pi in model.length_distribution.items():
        if L >= len(start) + len(end):
            clamp = {**dict(enumerate(start)), **{L - len(end) + i: x for i, x in enumerate(end)}}
            total += pi * np.exp(model.log_clamp_probability(L, clamp))
    return total
mixed = share(merged, "CASS", "TEAFF") + share(merged, "CSA", "YEQYF")
print(f"\nmerged model: {mixed:.1%} of its probability starts like one group and ends like the other")

# Diversity: Hill number of order 2, the inverse chance that two draws are the same sequence
print(f"D2: group A {group_a.hill_number(2):,.0f}, group B {group_b.hill_number(2):,.0f}, "
      f"mixture {mixture.hill_number(2):,.0f}, merged {merged.hill_number(2):,.0f}")
beta = mixture.beta_diversity(2.0)
print(f"beta diversity of the mixture: {beta.beta:.2f} (gamma {beta.gamma:,.0f} / alpha {beta.alpha:,.0f})")

Output

sequence           group A   group B   mixture    merged
CASSLGQGYEQYF      4.9e-04   4.1e-31   2.5e-04   1.7e-04
CSARGNTEAFF        1.1e-32   9.3e-04   4.7e-04   3.0e-04
CASSLGQGTEAFF      2.8e-19   2.6e-15   1.3e-15   6.2e-05

merged model: 47.5% of its probability starts like one group and ends like the other
D2: group A 44,304, group B 38,210, mixture 82,063, merged 127,208
beta diversity of the mixture: 2.00 (gamma 82,063 / alpha 41,032)

Full program, including the set-up hidden here

CASSLGQGTEAFF starts like group A and ends like group B. The mixture gives it almost nothing; the merged model gives it almost 50 billion times more, and puts almost half of all its probability on such recombined sequences. That is also why the merged model is more diverse than the two groups together.

Which one to use:

• Merge to describe a group as one population, with all of its evidence behind every row. The group models above and the reference model are merged models, and recombining the members’ patterns is what lets them generalise to sequences no member happened to contain.
• Mix when the members must stay distinct, above all to measure the diversity between them.

Diversity between members

Diversity is measured with Hill numbers. At order $q$,

$$D_q(P) \;=\; \Big(\sum_s P(s)^q\Big)^{1/(1-q)}$$

is the effective number of sequences; at $q = 2$ it is $1 / \sum_s P(s)^2$, the inverse of the chance that two draws are the same sequence. For a mixture of $M$ members:

• gamma, $D_q(P_{\text{mix}})$, is the diversity of the group as a whole;
• alpha is the diversity of a typical member: the members’ $D_q$ averaged by a power mean of order $1 - q$, which at $q = 2$ is their harmonic mean;
• beta $= \text{gamma} / \text{alpha}$ is the effective number of distinct members, from 1 when all members are identical to $M$ when they share no sequence.

The two groups above share no sequence, so beta is 2.00: gamma 82,063 against alpha 41,032. The merged model cannot take gamma’s place: its 127,208 effective sequences include the recombined ones, which belong to no member, so the ratio would no longer count members.

With stand-in groups that overlap, beta falls between the two extremes:

mixture.py

from BEAM import BeamModel, Mixture
# group_a, group_a2: two groups from the same population; group_b, group_c: two others
same = Mixture([group_a, group_a2])
different = Mixture([group_a, group_b, group_c])

# beta diversity: the effective number of distinct groups, from 1 (identical) to the count
print(f"two groups, one population:      beta = {same.beta_diversity(2.0).beta:.3f}")
print(f"three groups, three populations: beta = {different.beta_diversity(2.0).beta:.3f}")

Output

two groups, one population:      beta = 1.035
three groups, three populations: beta = 1.563

Full program, including the set-up hidden here

Two groups drawn from one population count as one (1.04); three groups drawn from three different populations count as 1.56. At order 2 the measure is driven by the most common sequences, which the three stand-in populations largely share; their differences in rarer sequences show at lower orders and in the distance.

Rule of thumb

Compare samples with per-repertoire models, groups with group models (or with the distances between their members’ models), and use a reference for questions about the population as a whole. Merge models to describe a group as one population; mix them to keep its members apart.