Skip to content
MuteJesterPublic

About

Probability model for biological sequences such as TCRβ CDR3s, learned from sequences alone with optional abundances: sequence probabilities, sequencing depth, publicness and repertoire comparison, summed over every combinatorially possible sequence without enumeration or simulation.

Topics

Resources

Contributing

Security policy

Stars

0 stars

Watchers

0 watching

Forks

Repository files navigation

BEAM

Documentation PyPI Python Tests DOI License: GPL v3

Every possible sequence, accounted for.

BEAM is a probability model of biological sequences, amino acid or nucleotide. It learns from the sequences alone, with their abundances if you have them, and then answers: how likely is this sequence? How many distinct sequences will a sample of a given depth hold, and how many will two samples share? How many people in a cohort will carry a sequence? How do two repertoires, cell types or conditions differ, and why? Each answer is summed over every combinatorially possible sequence, without listing or simulating a single one, so the same input always gives the same answer.

BEAM comes with a ready-made reference model of T cell receptor β (TCRβ) CDR3s, built from 38,538 human repertoires, and builds models from your own data in one line.

Paper: T. Konstantinovsky and G. Yaari. A computable probability landscape for T cell receptor repertoire sampling and comparison. (Submitted.)

What can I do with it?

Question What BEAM gives you
How likely is this sequence? Its probability, normalised over all possible sequences, including ones never seen, and which residues make it likely or unlikely
What do typical sequences look like? Realistic sequences drawn from the model
How deeply should I sequence? The expected number of distinct sequences at any depth, and how many two samples share
How many people share this sequence? The expected number of repertoires in a cohort that contain it, and its full distribution
How different are two repertoires, and why? A distance that separates differences in sequence length from differences in local sequence patterns, and names the patterns

Install

pip install beam-model

Ready-to-use packages exist for Linux (x86-64 and ARM), macOS (Apple Silicon, macOS 14 or newer) and Windows (x86-64), for Python 3.11 to 3.15.

Quick start

import BEAM

# downloads once (~90 MB), then cached
model = BEAM.load_pretrained("human_tcrb_cdr3_aa_foundation")

model.pseq("CASSLGQGYEQYF")        # probability of one sequence
# 3.3e-05

The reference model uses about 1 GB of memory once loaded.

The reference model arrives with BEAM 1.0.0, published together with the paper. Until then, everything else works on models you build from your own sequences: see Using your own data.

Examples

How likely is a sequence?

pseq returns the probability of a sequence under the model: for the reference model, the chance that one T cell drawn at random carries exactly this CDR3. It takes one sequence or a list.

model.pseq(["CASSLGQGYEQYF", "CASSPGTGELFF", "CWWWWWWWWF"])
# array([3.3e-05, 2.9e-05, 6.9e-23])

Sequences the model has never seen still get a sensible probability: a plausible new CDR3 gets a moderate one, an implausible one (like the last) a vanishingly small one. Use log_pseq_batch for millions of sequences at once.

BEAM reads every residue twice, once from the start of the sequence and once from the end, and explain shows how each residue contributes to the score:

e = model.explain("CASSLGQGYEQYF")
e["forward"]     # each residue's log probability, read from the start
e["backward"]    # ... and read from the end

Bar chart of how unexpected each residue of CASSLGQGYEQYF is, read from the start and from the end

How unexpected each residue of CASSLGQGYEQYF is, in bits, read from the start (teal) and from the end (orange), and the two combined. The conserved start and end are expected; the information sits in the junction between them.

What do typical sequences look like?

model.sample(5, seed=1)
# ['CASSPGTSSDTYEQYF', 'CASSLGTYSYNEQFF', 'CASRPGNSPLHF', ...]

Sequences are drawn in proportion to their probability. The model can also be drawn as sequence logos computed over every sequence of a length, without sampling (model.insights.logo):

Two sequence logos of length-15 CDR3s: common sequences above, with taller, more conserved stacks, and rare sequences below

Length-15 CDR3s under the reference model: the common sequences (probability above 10⁻⁷, a quarter of all the probability at this length) and the rare ones. Both keep the conserved CASS…F frame; the common ones are more constrained, with 37 bits of information against 29.

How deeply should I sequence?

If you sample n cells from a population, how many different sequences will you see? And if two labs each sample n, how many sequences will they have in common? BEAM computes both directly, without simulation.

from BEAM.analytics.exact_occupancy import ExactOccupancy

occupancy = ExactOccupancy(model)

# distinct sequences in 100,000 cells, and in a million
occupancy.richness(100_000).value
# 96771.4...
occupancy.richness(1_000_000).value
# 879776.1...

# sequences shared by two samples of 100,000 cells
occupancy.sharing(100_000, 100_000).value
# 3678.4...

Curve of the share of one sample's sequences also found in a second sample, rising with the number of cells sampled

The share of one sample's distinct sequences that a second sample of the same size also contains, from a thousand to a billion cells per sample, under the reference model.

Each result also reports a guaranteed bound on its relative error (.bound), here below one part in ten billion. The first query takes a minute or two while BEAM prepares the integrals it needs; every later query, at any depth, is instant.

How many people will share a sequence?

A sequence is public when many people carry it. Describe a cohort by how deeply each person's repertoire was sampled, and BEAM tells you how many of them are expected to carry a sequence.

from BEAM import PublicnessModel

# 100 people, 50,000 cells sampled from each
cohort = PublicnessModel([50_000] * 100)

mean, variance = cohort.moments(model.pseq("CASSLGQGYEQYF"))
mean                     # expected number of people who carry it
# 80.8

cohort.pmf(p) gives the full distribution: the probability that exactly 0, 1, 2, ... people carry the sequence.

Distributions of the number of people, out of 100, who carry each of three CDR3s

How many of the 100 people carry each of three CDR3s: a public one, found in about 80 of them, one found in about 27, and a rare one that almost nobody carries.

How different are two repertoires, and why?

Build a model from each person's repertoire and compare them. The distance splits into two parts that add up to it, one due to different sequence lengths and one due to different local sequence patterns, and it names the patterns that differ most.

import pandas as pd
from BEAM import BeamModel
from BEAM.analytics.distance import compare

def repertoire_model(path):
    table = pd.read_csv(path, sep="\t")             # an AIRR-format table
    return BeamModel.from_dataframe(table, sequence_column="junction_aa",
                                    abundance_column="duplicate_count")

a = repertoire_model("donor_A.tsv")
b = repertoire_model("donor_B.tsv")

result = compare(a, b)
result.distance          # 0 = identical
result.patterns[:5]      # the patterns that differ most

# the share of the distance due to lengths
result.length_squared_distance / result.squared_distance

Stacked bars per CDR3 length of the length part and the local-pattern part of the squared distance between two repertoires

Where the distance between two repertoires comes from, length by length: a typical repertoire against one with more long CDR3s and more …NEQFF endings. Most of it is the difference in lengths (orange); the rest is local patterns (teal), led by NEQFF.

Comparing two repertoires takes under a minute.

Using your own data

from BEAM import BeamModel

# from Python lists
model = BeamModel(sequences, abundances=counts)

# from a file with one "SEQUENCE<TAB>COUNT" per line
model = BeamModel.from_sequence_file("cdr3.tsv")

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

Abundances are optional; without them every sequence counts once. Nucleotide sequences work too, with alphabet="nt".

To score sequences against the reference model, or to compare your models with it, write CDR3s as it does: including the conserved cysteine at the start and the phenylalanine or tryptophan at the end, for example CASSLGQGYEQYF. In AIRR-format files this is the junction_aa column; AIRR's cdr3_aa column leaves those two residues out.

Why no simulation?

Most questions about a repertoire are sums over every possible sequence, and there are far too many sequences to list. The usual workaround is simulation: draw many sequences and count. BEAM is built so that these sums can be computed directly, by dynamic programming over the whole sequence space. Its probabilities are normalised over all sequences, and quantities such as expected richness, overlap and publicness are evaluated numerically, with stated error bounds where a sum is approximated. The same input always gives the same answer, and small differences between groups are not lost in simulation noise.

Documentation

The documentation website, mutejester.github.io/BEAM, has guides for each question above, explanations of how BEAM works with interactive figures, and the API reference.

Questions and ideas are welcome in Discussions; bugs in Issues.

Citing BEAM

If BEAM is useful in your work, please cite the paper above, and the software itself by its DOI, 10.5281/zenodo.23470065 (every version; each release also has its own). GitHub's "Cite this repository" button (from CITATION.cff) gives the reference in BibTeX and APA.

License

BEAM is free software under the GNU General Public License v3.0 or later. Contributions are welcome: see CONTRIBUTING.md.

About

Probability model for biological sequences such as TCRβ CDR3s, learned from sequences alone with optional abundances: sequence probabilities, sequencing depth, publicness and repertoire comparison, summed over every combinatorially possible sequence without enumeration or simulation.

Topics

Resources

Contributing

Security policy

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Used by

Contributors

Languages