How likely is a sequence?

The question: how probable is this sequence in a population, and why?

A BEAM model gives every sequence a probability: the chance that one sequence drawn at random from the population the model describes is exactly that one. With the reference model, it is the chance that one T cell drawn at random from the 38,538 human repertoires it was built from has exactly that CDR3β; the examples on this page use it.

One sequence, or a few

pseq returns the probability; pseq(..., log=True) and log_pseq_batch return its natural logarithm, which is easier to compare because CDR3 probabilities are tiny.

one_and_many.py

import BEAM

model = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")

# One sequence: a probability, or its natural logarithm.
print(model.pseq("CASSLGQGYEQYF"))
print(model.pseq("CASSLGQGYEQYF", log=True))

# Many sequences at once.
cdr3s = ["CASSLGQGYEQYF", "CASSPGTGELFF", "CSARDRTGNGYTF"]
for cdr3, logp in zip(cdr3s, model.log_pseq_batch(cdr3s)):
    print(f"{cdr3:16} log P = {logp:7.3f}")

Output

3.29567201803869e-05
-10.320315366964564
CASSLGQGYEQYF    log P = -10.320
CASSPGTGELFF     log P = -10.453
CSARDRTGNGYTF    log P = -12.019

Full program, including the set-up hidden here

A log probability of −10.3 is a probability of about 3 × 10⁻⁵: roughly 3 cells in every 100,000 carry this sequence.

Millions of sequences

log_pseq_batch scores a whole list in one call and is the way to score a repertoire.

millions.py

import numpy as np
import BEAM

model = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")

# A large batch: 200,000 sequences, here drawn from the model itself.
cdr3s = model.sample(200_000, seed=7)
logp = model.log_pseq_batch(cdr3s)

print(f"scored {len(logp):,} sequences")
print(f"median probability:     {np.exp(np.median(logp)):.2e}")
print(f"5% are rarer than:      {np.exp(np.percentile(logp, 5)):.2e}")

Output

scored 200,000 sequences
median probability:     9.85e-09
5% are rarer than:      1.05e-13

Full program, including the set-up hidden here

Half of the sequences have a probability below about 10⁻⁸. That is the typical scale: most CDR3s are individually rare, and only a small set of sequences reach 10⁻⁵ or above.

Comparing sequences of different lengths

A longer CDR3 is usually less probable, partly because it is longer and partly because long CDR3s are rarer to begin with. The two parts separate exactly:

$$P(s) = P(L)\,P(s \mid L)$$

where $P(L)$ is the probability of the length and $P(s \mid L)$ the probability of the sequence among all sequences of that length. length_distribution gives $P(L)$.

lengths.py

import BEAM

model = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")
p_length = model.length_distribution          # {length: probability}

for cdr3 in ["CASSLGQGYEQYF", "CASSLGQGAYNEQFF", "CASSPRDRGLAGGYNEQFF"]:
    p = model.pseq(cdr3)
    L = len(cdr3)
    print(f"{cdr3:20} L={L}  P(s) = {p:.2e}   P(L) = {p_length[L]:.3f}"
          f"   P(s | L) = {p / p_length[L]:.2e}")

Output

CASSLGQGYEQYF        L=13  P(s) = 3.30e-05   P(L) = 0.171   P(s | L) = 1.93e-04
CASSLGQGAYNEQFF      L=15  P(s) = 6.44e-06   P(L) = 0.237   P(s | L) = 2.71e-05
CASSPRDRGLAGGYNEQFF  L=19  P(s) = 1.60e-09   P(L) = 0.012   P(s | L) = 1.28e-07

Full program, including the set-up hidden here

The 19-residue sequence is about 20,000 times less probable than the 13-residue one. Only about 14 times of that comes from its length; the rest is the sequence itself. Comparing $P(s \mid L)$ asks “how typical is this sequence for its length”.

Which residues drive a score

BEAM scores every residue twice, by two experts: one reads the sequence from the start, the other from the end. For residue $s_i$ of a sequence of length $L$, using the $k$ neighbouring residues as context ($k = 4$ for the reference model):

$$\psi_i \;=\; q_F\big(s_i \mid i,\; s_{i-k}^{\,i-1}\big)^{w_F}\; q_B\big(s_i \mid L-i+1,\; s_{i+1}^{\,i+k}\big)^{w_B}, \qquad P(s) \;=\; \frac{\pi(L)}{Z(L)} \prod_{i=1}^{L} \psi_i$$

• $q_F(s_i \mid i, \dots)$ is the probability the expert reading from the start gives residue $s_i$, from its position $i$ and the $k$ residues before it;
• $q_B(s_i \mid L-i+1, \dots)$ is the same for the expert reading from the end, from the position counted from the end and the $k$ residues after it;
• $\psi_i$ is the residue’s combined factor: the geometric mean of the two, with weights $w_F = w_B = \tfrac{1}{2}$;
• $\pi(L)$ is the probability of the length, and $Z(L)$ the normaliser: the product of the $\psi_i$ summed over every sequence of length $L$.

explain returns all of these as natural logarithms: for each residue forward ($\ln q_F$), backward ($\ln q_B$) and combined ($\ln \psi_i$), and for the sequence log_length_prior ($\ln \pi(L)$) and log_partition ($\ln Z(L)$). Below, each residue’s terms are shown as surprisal in bits, $-\log_2$: 0 means the residue was entirely expected, and each extra bit means half as likely.

explain.py

import math
import numpy as np
import BEAM

model = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")
e = model.explain("CASSLGQGYEQYF")

bits = lambda log_p: -np.asarray(log_p) / math.log(2)          # surprisal, in bits
print("  i  residue   from start   from end   combined (psi_i)")
for i, (res, f, b, c) in enumerate(zip(e["residues"], bits(e["forward"]),
                                       bits(e["backward"]), bits(e["combined"])), start=1):
    print(f"{i:3}  {res:^7}  {f:11.2f}  {b:9.2f}  {c:17.2f}")

# ln P(s) = ln pi(L) - ln Z(L) + sum_i ln psi_i
total = e["log_length_prior"] - e["log_partition"] + e["combined"].sum()
print(f"\nln pi(L) = {e['log_length_prior']:.3f}   ln Z(L) = {e['log_partition']:.3f}   "
      f"sum of ln psi_i = {e['combined'].sum():.3f}")
print(f"ln pi(L) - ln Z(L) + sum of ln psi_i = {total:.6f}")
print(f"log_pseq                             = {e['log_pseq']:.6f}")

Output

  i  residue   from start   from end   combined (psi_i)
  1     C            0.00       0.00               0.00
  2     A            0.17       0.05               0.11
  3     S            0.15       0.08               0.11
  4     S            0.25       0.77               0.51
  5     L            1.99       3.27               2.63
  6     G            2.14       0.92               1.53
  7     Q            3.08       3.41               3.24
  8     G            0.91       3.24               2.08
  9     Y            4.25       0.49               2.37
 10     E            1.73       1.07               1.40
 11     Q            0.14       0.15               0.15
 12     Y            0.07       1.42               0.74
 13     F            0.03       0.01               0.02

ln pi(L) = -1.765   ln Z(L) = -1.767   sum of ln psi_i = -10.322
ln pi(L) - ln Z(L) + sum of ln psi_i = -10.320315
log_pseq                             = -10.320315

Full program, including the set-up hidden here

Surprisal of each residue of CASSLGQGYEQYF under the two readings of the reference model (bars) and combined (line).

The reading from the start is sure of the opening CASS, the reading from the end is sure of the closing Q and F, and both are uncertain in the middle (LGQGY), where the junction between the gene segments makes CDR3s diverse. Where one expert is sure and the other is not, as at the Y in position 9, the combined factor sits between them. The terms reconstruct the score exactly, as the last lines show, so this is an explanation of the probability, not an illustration of it.

A $\psi_i$ is a factor of the score, not a probability of the residue on its own: the product of the two experts is normalised as a whole, by $Z(L)$. Because $Z(L)$ is the same for every sequence of length $L$, the $\psi_i$ of different residues, or of one position in different sequences of the same length, can be compared directly.

What the model cannot score

A sequence longer than the model’s maximum length (27 residues for the reference model), or containing a character outside the model’s alphabet (the 20 standard amino acids for the reference model), has probability 0. explain raises an error that says why.

invalid.py

import BEAM

model = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")

too_long = "CASS" + "G" * 30 + "YEQYF"     # 39 residues; the model covers up to 27
unknown = "CASSLGQGXEQYF"                 # X is not one of the 20 amino acids
print(model.pseq(too_long), model.pseq(unknown))

try:
    model.explain(unknown)
except ValueError as err:
    print("explain:", err)

Output

0.0 0.0
explain: sequence contains ['X'], which the 'aa' alphabet does not accept; its residues are 'ACDEFGHIKLMNPQRSTVWY'

Full program, including the set-up hidden here

Pitfalls

• Score sequences written as the model’s were. A model knows only the form of the sequences it was built from, so trimmed or extended ends, or a different reading frame, score as different, rarer strings. For the reference model, that form is junction_aa, with the conserved start (C) and end (F or W); see Writing CDR3 sequences.
• Probabilities describe the population the model was built from. Under the reference model, “likely” means common across many published human repertoires, not likely in any one person.
• Compare like with like. Rank sequences of one length by $P(s)$, or compare across lengths with $P(s \mid L)$, depending on whether length should count.