Building a model from your data

The question: how do I get a model of my repertoire, so I can score, sample and compare it?

A model is built by counting short stretches of sequence across a repertoire, weighted by how many cells carry each sequence. Building takes seconds: on a 16-core machine, 4 seconds for a repertoire of 50,000 cells and 8 seconds for 2 million cells (1.7 million distinct sequences).

The examples on this page use a stand-in repertoire: 50,000 cells drawn from the reference model. Replace it with your own sequences.

From Python lists

Give the sequences, here CDR3s, and, optionally, how many cells (or templates, or reads) carry each. Without abundances every sequence counts once.

from_lists.py

from collections import Counter
import BEAM
from BEAM import BeamModel

# sequences: list of CDR3 strings; abundances: how many cells carry each
model = BeamModel(sequences, abundances=abundances)

print(model)
print(f"distinct sequences: {len(sequences):,}   cells: {sum(abundances):,}")
print(f"P(CASSLGQGYEQYF) = {model.pseq('CASSLGQGYEQYF'):.2e}")

Output

BeamModel(order=4, sequences=49,063, states=194,481, memory=53MB)
distinct sequences: 49,063   cells: 50,000
P(CASSLGQGYEQYF) = 1.87e-05

Full program, including the set-up hidden here

The model of this repertoire gives CASSLGQGYEQYF a probability of 1.9 × 10⁻⁵, against 3.3 × 10⁻⁵ under the reference model: 50,000 cells carry less information than 38,538 repertoires, and the model reflects what this sample contains.

How long a sequence can be

A model represents sequences of up to max_length residues, and skips longer ones when it is built. The default is 27 for amino acids and 81 for nucleotides. max_length="data" takes the length of the longest sequence you give, and a number sets it yourself.

max_length.py

from collections import Counter
import BEAM
from BEAM import BeamModel

# sequences include two of 31 residues, longer than the default maximum of 27
default = BeamModel(sequences, abundances=abundances)                    # up to 27 residues
fitted = BeamModel(sequences, abundances=abundances, max_length="data")  # the longest given
chosen = BeamModel(sequences, abundances=abundances, max_length=40)      # a number you choose

for name, m in [("default", default), ('"data"', fitted), ("40", chosen)]:
    print(f"max_length={name:9}  represents up to {m.config['max_length']} residues,"
          f" skipped {m.n_skipped}")

Output

max_length=default    represents up to 27 residues, skipped 2
max_length="data"     represents up to 31 residues, skipped 0
max_length=40         represents up to 40 residues, skipped 0

Full program, including the set-up hidden here

Time and memory grow with max_length, and every sequence up to it gets a probability, so set it no longer than your sequences need. from_sequence_file and from_parquet read a file once without holding it, so they cannot take the length from the data: give them a number.

From an AIRR table

AIRR-format rearrangement tables (.tsv) can be read with pandas and passed directly. Name the columns that hold the sequence and the count.

from_airr.py

import pandas as pd
from BEAM import BeamModel

table = pd.read_csv("repertoire.tsv", sep="\t")      # an AIRR rearrangement table
model = BeamModel.from_dataframe(table, sequence_column="junction_aa",
                                 abundance_column="duplicate_count")
print(model)
print(f"sequences skipped: {model.n_skipped}")

Output

BeamModel(order=4, sequences=49,063, states=194,481, memory=53MB)
sequences skipped: 0

Full program, including the set-up hidden here

n_skipped counts the rows that were left out: an empty sequence, one longer than the maximum length, or one with a character outside the alphabet. Checking it catches a wrong column straight away.

Two other readers avoid holding the table in memory:

model = BeamModel.from_sequence_file("cdr3.tsv")   # one "SEQUENCE<TAB>COUNT" per line
model = BeamModel.from_parquet("repertoire.parquet",
                               sequence_column="junction_aa",
                               abundance_column="duplicate_count")   # needs pyarrow

Saving and loading

A saved model keeps its counts and settings, so it reloads to exactly the same probabilities.

save_load.py

from BEAM import BeamModel

model.save("donor_A")                          # writes donor_A.beam
again = BeamModel.from_file("donor_A.beam")

print(again)
print(again.pseq("CASSLGQGYEQYF") == model.pseq("CASSLGQGYEQYF"))

Output

BeamModel(order=4, sequences=20,000, states=194,481, memory=32MB)
True

Full program, including the set-up hidden here

Combining models

There are two different ways to combine repertoires, and they answer different questions.

combine.py

from BEAM import BeamModel
# a and b: models built from two donors' repertoires
pooled = a.merged(b)                        # pool the counts: one bigger repertoire

s = "CASSLGQGYEQYF"
mixed = 0.5 * a.pseq(s) + 0.5 * b.pseq(s)   # mix the probabilities: pick a donor, then a cell
print(f"donor A {a.pseq(s):.3e}   donor B {b.pseq(s):.3e}")
print(f"pooled  {pooled.pseq(s):.3e}   mixed   {mixed:.3e}")

Output

donor A 2.079e-05   donor B 1.907e-05
pooled  2.026e-05   mixed   1.993e-05

Full program, including the set-up hidden here

• Pooling the counts (a.merged(b)) gives the model of one combined repertoire, as if all the cells had been sequenced together.
• Mixing the probabilities (here written out as the average) describes picking a donor at random and then a cell from that donor. It is the right object for questions about a group of people.

The two differ because pooling counts can combine stretches of sequence seen in different donors into sequences no single donor had. For many donors, BEAM.Mixture works with probability mixtures of models directly.

Nucleotide sequences

The same model works on nucleotide sequences: name the alphabet when building.

model = BeamModel(nucleotide_cdr3s, alphabet="nt")   # A, C, G, T and N

Pitfalls

• Choose the column deliberately. Without sequence_column, from_dataframe and from_parquet read cdr3_aa, which in AIRR tables leaves out the conserved C and F/W. To compare with the reference model, pass sequence_column="junction_aa".
• Abundance units differ between platforms. Templates, UMIs and reads are not the same thing. Within one study use one unit; across studies keep in mind that counts are not directly comparable.
• Long sequences are dropped, not truncated. The default maximum length is 27 for amino acids and 81 (27 codons) for nucleotides; set max_length, or max_length="data", as in How long a sequence can be. Check n_skipped after building.