OLGA: fast computation of generation probabilities of B- and T-cell receptor amino acid sequences and motifs

Motivation: High-throughput sequencing of large immune repertoires has enabled the development of methods to predict the probability of generation by V(D)J recombination of T- and B-cell receptors of any specific nucleotide sequence. These generation probabilities are very non-homogeneous, ranging over 20 orders of magnitude in real repertoires. Since the function of a receptor really depends on its protein sequence, it is important to be able to predict this probability of generation at the amino acid level. However, brute-force summation over all the nucleotide sequences with the correct amino acid translation is computationally intractable. The purpose of this paper is to present a solution to this problem. Results: We use dynamic programming to construct an efficient and flexible algorithm, called OLGA (Optimized Likelihood estimate of immunoGlobulin Amino-acid sequences), for calculating the probability of generating a given CDR3 amino acid sequence or motif, with or without V/J restriction, as a result of V(D)J recombination in B or T cells. We apply it to databases of epitope-specific T-cell receptors to evaluate the probability that a typical human subject will possess T cells responsive to specific disease-associated epitopes. The model prediction shows an excellent agreement with published data. We suggest that OLGA may be a useful tool to guide vaccine design. Availability: Source code is available at https://github.com/zsethna/OLGA


I. INTRODUCTION
The ability of the adaptive immune system to recognize foreign peptides, while avoiding self peptides, depends crucially on the specificity of receptor-antigen binding and the diversity of the receptor repertoire. Immune repertoire sequencing (Repseq) of B-and T-cell receptors (BCR and TCR) [16,21,39,46] offers an efficient experimental tool to probe the diversity of full repertoires in healthy individuals [11,18,26,28,32,33,45], in cohorts with specific conditions [4,9,10,17,19,20,29,43] and evaluate the response to specific fluorescent MHCmultimers [2,14]. Recent work has shown that responding clonotypes often form disjoint clusters of similar amino acid sequences, which has lead to the identification of responsive amino acid motifs [2,14]. In order for these techniques to have practical applications in therapy and vaccine design, one needs a fast and efficient algorithm to evaluate which specific amino acid sequences and sequence motifs are likely to be generated and found in repertoires. We present a solution to this problem in the form of an algorithm and computational tool, called OLGA, which implements an exact computation of the generation probability of any BCR or TCR sequence (nucleotide or amino acid), or motif. BCR and TCR are stochastically generated by choosing a germline genetic template in each of several cassettes of alternates (V, (D), or J) and then splicing them together with random nucleotide deletions and insertions at the junctions. Given a generative model, one can define the generation probability of any nucleotide sequence as the sum of the probabilities of all the generative events that can produce that sequence [6,7,25,27]. However, computing the generation probability of amino acid sequences by summing over all consistent nucleotide sequences is impractical: because of codon degeneracy, the number of nucleotide sequences to be summed grows exponentially with sequence length. OLGA is powered by an efficient dynamic programming method to exactly sum over generative events and obtain net probabilities of amino acid sequences and motifs.
We validate our algorithm by comparing its results and performance to Monte-Carlo sampling estimates. We present results using publicly available data for both TCR α (TRA, Pogorelyy et al. [28]) and β (TRB, Robins et al. [33]) chains and BCR heavy chains (IGH, DeWitt et al. [3] of humans), and TRB of mice [35]. We applied OLGA to a TCR database that catalogs the different CDR3 amino acid sequences responding to a variety of different epitopes associated with disease [37]. We computed the generation probability of particular CDR3 amino acid sequences, as well as the net generation probability of all the TCR that respond to a particular epitope. Finally, we discuss OLGA's applications in vaccine design and other therapeutic contexts.

II. METHODS
A. Stochastic model of VDJ recombination V(D)J recombination is a stochastic process involving several events (gene template selection, terminal deletions from the templates, random insertions at the junctions), each of which has a set of possible outcomes chosen according to a discrete probability distribution. The probability P rec gen (E) of any generation event E, defined as a combination of the above-mentioned processes is, for the TRB locus: where (V, D, J) identify the choices of gene templates, (d V , d D , d D , d J ) are the numbers of deletions at each end of the segments, and (m 1 , . . . , m V D ) and (n 1 , . . . , n DJ ) are the untemplated inserted nucleotide sequences at the VD and DJ junctions. These variables specify the recombination event E, and are drawn according to the probability distributions (P V , P DJ , P delV , P delD , P delJ , P insVJ , P insDJ , p 0 , q 0 , S VD , S DJ ). The inserted segments are drawn according to a Markov process starting with the nucleotide distribution p 0 and with the transition matrix R, and running from the 5' side (left to right) for the VD segment, and from the 3' side (right to left) from the DJ segment. Similar models can be defined for the α chain or for BCR chains. Although here we describe the method for TRB only, it is also implemented for other chains in the software. Since the same nucleotide sequence can be created by more than one specific recombination event, the generation probability of a nucleotide sequence is the sum of the probabilities of all possible events that generate the sequence: P nt gen (σ) = E→σ P rec gen (E), where the sum is over all recombination events E that produce the sequence σ = (σ 1 , . . . , σ n ). The probability of generation of an amino acid sequence, a = (a 1 , . . . , a L ) is the sum of the probabilities of all nucleotide sequences that translate into the amino acid sequence: where the ∼ sign indicates that σ translates into a. We can generalize this approach to any scheme that groups nucleotide triplets, or codons, into arbitrary classes, which we still denote by σ ∼ a. In the formulation above, these classes simply group together codons with the same translation according to the standard genetic code. In an example of generalization, all codons that code for amino acids with a common chemical property, e.g. hydrophobicity or charge, could be grouped into a single class. In that formulation, (a 1 , . . . , a L ) would correspond to a sequence of symbols denoting that property. More generally, any grouping of amino acids can be chosen (including one where any amino acid is acceptable), and the partition can be position dependent. Thus, the generation probability of arbitrary "motifs" can be queried. In the following, for ease of exposition, we restrict our attention to the case where a is an amino acid sequence.
B. Dynamic programming computation of the generation probability of amino acid sequences We now give an overview of how OLGA computes Eq. 2 without performing the sum explicitly, using dynamic programming. Fig. S1-S2 give a graphical overview of the method, and details of the method implementation can be found in SI Secs. I and II and in the code manual. Given the genomic nucleotide sequences of the possible gene templates, together with a specific model of the type described in Eq. A1, the algorithm computes the net probability of generating a recombined gene with a given CDR3 amino acid sequence under a given set of V and J gene choices.
Each recombination event implies an annotation of the CDR3 sequence, assigning a different origin to each nucleotide (V, N1, D, N2, or J, where N1 and N2 are the VD and DJ insertion segments, respectively) that parses the sequence into 5 contiguous segments (see schematic in Fig. 1). The principle of the method is to sum over the probabilities of all choices of nucleotides consistent with the known amino acid sequence, over the possible locations of the 4 boundaries (x 1 , x 2 , x 3 , and x 4 ) between the 5 segments, and over the possible V, D, and J genomic templates (Fig. 1). We do this in a recursive way using matrix operations by defining weights that accumulate the probabilities of events from the left of a position x (i.e. up to x), and weights that accumulate events from the right of x (i.e. from x + 1 on). Specifically, we define the following index notation: X x with a subscript called left index, accumulates weights from the left of x; Y x , with a superscript called right index, accumulates weights from the right of x; a matrix X x y corresponds to accumulated weights from position x + 1 to y (as will be explained shortly, these objects may have suppressed nucleotide indices as well). P aa gen is calculated recursively by matrix-like multiplications as: (3) The vector V x corresponds to a cumulated probability of the V segment finishing at position x; M x y is the probability of the VD insertion extending from x + 1 to y; N x y is the same for DJ insertions; D x y (D) corresponds to weights of the D segment extending from x + 1 to y, conditioned on the D germline choice being D; J x (D) gives the weight of J segments starting at position x + 1 conditioned on the D germline being D. This D dependency is necessary to account for the dependence between the D and J germline segment choices [27]. All the defined vectors and matrices depend implicitly on the amino acid sequence (a 1 , . . . , a L ), but we leave this dependency implicit to avoid making the notation too cumbersome.
Because we are dealing with amino acid sequences encoded by triplet nucleotide codons, we need to keep track of the identity of the nucleotide at the beginning or the end of a codon. Depending on the position of the index x in the codon, the objects defined above may be vectors of size 4 (or 4 × 4 matrices) in the suppressed nucleotide index. We use conventions that depend on whether we are considering left or right indices, as follows.
If x is a multiple of 3, i.e. x = 0 (mod 3), then we do not keep nucleotide information and both X x and Y x are scalars (whether x is a left or a right index). If x = 1 (mod 3), then X x must be interpreted as a row vector of 4 numbers, X x (σ), σ = A, T, G, C, corresponding to the cumulated probability weight that the nucleotide at position x (first position of the codon) takes value σ. If x = 2 (mod 3), then X x is also a row vector of 4 numbers, X x (σ), but with a different interpretation: it corresponds to the cumulated probability up to position x, with the additional constraint that the nucleotide at position x + 1 (the last position in the codon) can take value σ (the value is 0 otherwise). For right indices, the interpretation is reversed and the entries are column vectors: when x = 1 (mod 3) the Y x is a column vector containing the cumulated weights from x + 1 onwards, with the constraint that the nucleotide at x can be σ, and when x = 2 (mod 3), it is the probability weight that the nucleotide at position x + 1 is σ. Generalizing to matrices, X x y is a 4x4, 4x1, 1x4, or 1x1 matrix depending on whether the x and y positions are multiples of 3 or not, with the same rules as for vectors for each type of index.
Entries with left indices are interpreted as row vectors, and entries with right indices as column vectors. Thus, in Eq. B2 contractions between left and right indices correspond to dot products over the 4 nucleotides when the index is not a multiple of 3, and simply a product of scalars when it is.
The entries of the matrices corresponding to the germline segments, V, D(D), and J (D), can be calculated by simply summing over the probabilities of different germline nucleotide segments compatible with the amino acid sequence (a 1 , . . . , a L ) with conditions on deletions to achieve the required segment length. For instance, the V matrix elements are given by: x is the u th nucleotide of the i th codon, s V the sequence of the V germline gene, and I the indicator function. The ∼ sign is generalized to incomplete codons so that it returns a true value if there exists a codon completion that agrees with the motif a.
Detailed formulas for the other segments are derived using the same principles and are given in the SI Appendix. The sums in Eq. 4 (and equivalent expressions for J) can be restricted to particular germline genes to compute the generation probability of particular VJ-CDR3 combinations. The entries of the insertion segment N1 are calculated using the following formula: with y = 3(j − 1) + v (and x = 3(i − 1) + u as in Eq. 4). The transfer matrix corresponds to the probability of inserting a codon coding for a and ending with nucleotide σ, knowing that the previous codon ended with nucleotide τ . L u a and R v a are vectors or matrices with different definitions depending on the values of x and y modulo 3, corresponding to the probabilities of inserting incomplete codons on the left and right ends of the insertion segment. Eq. 5 is only valid for j > i, but similar formulas describe the case i = j. The precise definitions of L and R, the i = j case, and the formulas for N and the N2 insertion segment, which is exactly equivalent, are all given in detail in the SI Appendix. The matrix product of Eq. 5 can be calculated recursively, requiring only 4 × 4 matrix multiplications. Thus, all M x y elements can be calculated in O(L 2 ) operations, instead of the exponential time that would be required using brute-force summation over nucleotides in degenerate codons. Finally, since the sums of Eq. B2 can also be done recursively through L × L matrix operations, the whole procedure has O(L 2 ) computational complexity.

A. Method validation
To verify the correctness of the OLGA code, we compared its predictions for generation probabilities to those estimated by Monte Carlo (MC) sequence generation [29]. MC estimation is done by drawing events from a given generative model, binning according to the resulting CDR3 amino acid sequence, and normalizing by the total number of recombination events. The scatter plot of the estimated generation probabilities for these sequences against the values predicted by OLGA gives a direct test of the algorithm. As MC estimation is susceptible to Poisson sampling noise, it is important to ensure that enough events are drawn to accurately assess the generative probabilities of individual CDR3 sequences. For this reason, we made the comparison using a generative model inferred from a mouse, rather than human, T cell repertoire, because of the significantly lower entropy of mouse repertoires [35]. The specific model was inferred by IGoR [25] using ∼ 70000 out-of-frame TRB sequences from a mature mouse thymus. MC estimation was done by generating 5 × 10 11 recombination events, from which the first 10 6 unique CDR3 amino acid sequences are counted to serve as a sample for the comparison. This procedure provided good sequence coverage, with > 98% of sequences generated at least twice and > 95% of sequences generated at least 10 times. As Fig. 2 shows for mouse TRB (see Fig. S3 for human TRA), MC estimation and OLGA calculation are in agreement (up to Poisson noise in the MC estimate). The Kullback-Leibler divergence between the two distributions, a formal measure of their agreement, is a mere 4.82 × 10 −7 bits.

B. Comparison of performance with existing methods
We compared the performance of OLGA to other methods. Direct calculation of amino acid sequence generation probability using OLGA is orders of magnitude faster than the two possible alternative methods: MC estimation (as described above), or exhaustive enumeration of the generative events giving rise to a given amino acid sequence. OLGA took 6 CPU hrs to compute the generation probabilities of the 10 6 amino acid sequences, i.e. 47 seqs/CPU/sec for mouse TRB (see SI Sec. III and Table S1 for runtimes of other loci). By comparison, MC estimation required 4313 CPU hrs. The scaling for the MC estimation does not depend on the number of queried sequences, but instead is determined by the number of recombinations needed to control the Poisson noise, which scales inversely with generation probability. In practice, to determine the generation probability of a typical sequence (which can be as low 10 −20 , see Fig. 3 and below), one needs to generate very large datasets, and thus the generation probability of many sequences cannot be calculated by the MC method.
Alternatively, one could list all possible nucleotide sequences that translate to a particular amino acid CDR3 and sum the generation probabilities of each nucleotide sequence, using the IGoR algorithm [25]. Each amino acid sequence in the mouse validation sample is, on average, coded for by 1.84 billion nucleotide sequences (and much more for human TRB). Since IGoR computes generation probabilities of nucleotide sequences at the rate of ∼ 60 seqs/CPU/sec, it would take ∼ 8500 CPU hrs to compute the generation probability of a single amino acid sequence. A systematic comparison of OLGA with IGoR ( Fig. S4) and MC estimation (Figs. S4 and S5) as a function of the number of analysed sequences and their CDR3 lengths shows that OLGA is faster than both other methods for all practical purposes (see Sec. IV for details). , nucleotide CDR3 sequences (P nt gen ), and CDR3 amino acid sequences (P aa gen ) in different contexts. Each curve is determined by Monte Carlo sampling of 10 6 productive sequences for the indicated locus, and computing its generation probabilities at the three different levels. Entropies in bits (S) are, up to a ln(2)/ ln(10) factor, the negative of the mean of each distributions, indicated by dotted lines.

C. Distribution of generation probabilities and diversity
V(D)J recombination produces very diverse repertoires of nucleotide sequences, with a very broad distribution of generation probabilities spanning up to 20 orders of magnitude [6,27]. This distribution gives a comprehensive picture of the diversity of the process, and can be used to recapitulate many classical diversity measures [26], and to predict the overlap between the repertoires of different individuals [8]. In particular, the opposite of the mean logarithm of the generation probability, − log 2 P gen , is equal to the entropy of the process. While previous work focused on nucleotide sequence generation, OLGA allows us to compute this distribution for amino acid sequences. Fig. 3 shows the distribution of P aa gen for 4 loci: human and mouse TRB, human TRA, and human IGH, and compares it to the distributions of nucleotide sequence generation probabilities, P nt gen , and recombination event probabilities, P rec gen . While all these datasets are based on DNA RepSeq, we checked that the generation probability distribution was robust to the choice of protocol by computing the TRB distribution for independent datasets generated by RNA RepSeq [38,44,47] (Figs. S6 and S7, and SI Sec. V). The generation models used here and elsewhere in this paper were taken from Marcou et al. [25], except for the human TRB model which was relearned using IGoR from one individual in Emerson et al. [9] as a check. Going from recombination events to nucleotide sequences to amino acid sequences leads to substantial shifts in the distribution, and corresponding drops in entropies, as the distribution is progressively coarse-grained. Higher generation probability of a given receptor sequence leads to higher chance of finding it in any given individual. Generation probabilities may be constrasted to the scale set by the inverse of the number of independent recombination events (estimated between 10 8 [31] and 10 10 [22] for human TCR). Generation probabilities above this limit (10 −10 to 10 −8 for human TCR) can be considered "large" as the corresponding receptor will almost surely exist in each individual [8]. Another relevant scale to distinguish small from large generation probabilities is given by their geometric mean (dashed lines in Fig. 3).

D. Cross-species generation probabilities
While distinct species differ in their generation mechanisms, they may yet be able to generate the same CDR3s. Using OLGA, we computed the probabilities of producing human TRB CDR3s by the mouse recombination model, and vice versa (details in SI Sec. VI). An impressive 72.6% of human CDR3s can theoretically be produced by mice, and 100% of mouse CDR3s can be produced by humans. While cross-species generation probabilities are lower than intra-species ones (Fig. S8), they are correlated (Fig. S9). These results suggest that CDR3s observed in the repertoires of humanized mouse models of human diseases could be relevant for predicting their presence in human repertoires as well. OLGA allows for evaluating this potential, and could be used to inform clinical trials.

E. Generation probability of specific TCR
We can use OLGA to assess the total fraction of the generated repertoire that is specific to any given epitope, simply by summing the generation probabilities of all TRB sequences known to bind specifically to that epitope: where "a | epitope" means that the amino acid sequence a recognizes the epitope. Many experiments, based e.g. on multimer sorting assays [2,14] or T-cell culture assays, have established lists of epitope-specific TCR sequences for a number of disease-related epitopes. We used the VDJdb database [37], which aggregates such experiments, to compute P func gen of all TRB known to be reactive against several epitopes. In Fig. 4 we show results for 4 epitopes associated with Hepatitis C, and 5 epitopes associated with Influenza A. The net fraction of the repertoire specific to these epitopes (10 −7 to 10 −4 ) is large in the sense defined above, meaning that any individual is likely to have many copies of reactive T cells in their naive repertoire. : Generation probabilities of human CDR3s that respond to hepatitis C and influenza A epitopes. P aa gen of sequences that respond to an epitope are plotted as circles (color encodes density of the points). The fraction of the repertoire specific to each epitope (P func gen as defined in Eq. 7 ) is obtained as the sum of the P aa gen for each of the corresponding sequences (values plotted as triangles).
The presence of any specific TCR in the repertoire will be affected by the recombination probability of both its α and β chains, and also by function-dependent selective pressures. Assessing accurately the fraction of reactive TCRs in the blood is beyond the scope of this method. However, it is still interesting to ask whether epitopespecific TRB sequences had higher generation probabilities than regular sequences, either because of observational biases, or because the immune system might have evolved to make them more likely to be produced. To answer that question, we display in Fig. 5 the P aa gen distribution of the sequences listed in VDJdb that are specific to any epitope of each of 6 commonly studied viruses. For comparison we plot the P aa gen distribution of the full TRB sequence repertoire of a healthy donor (data taken from Emerson et al. [9]).
The viral distributions are very similar to each other, and also to the healthy repertoire background, meaning that the ability of a CDR3 to respond to a particular disease epitope is not strongly correlated with its generation probability. To see whether this result was confirmed in the case of a real infection, we repeated the same analysis on TRB RepSeq data from T-cells responding to three different types of pathogens (fungus, bacteria, and toxin) [1]. Consistently, we found that their distribution of generation probability was identical to that of naive sequences (SI Fig. S10 and SI Sec. VII). To compare OLGA's predictions with sequence occurrence frequencies in real data, we used the aggregated TRB repertoire of 658 human subjects described in Emerson  cally, we measured the frequencies in this large dataset of the specific CDR3 sequences contained in the VDJdb database [37], and compared them to the values assigned by OLGA. When measuring frequencies we discarded read count information, recording only the presence or absence of nucleotide sequences in each individual in order to eliminate effects of clonal expansion and PCR amplification bias, averaging over the 648 individuals in the Emerson et al. [9] dataset to get reliable estimates of frequencies. Each sequence in the VDJdb database is displayed as a dot in Fig. 6, and the resulting distribution shows a strong correspondence between mean frequency in the large data set and the predicted P aa gen of that sequence.
We then measured the fraction of CDR3s in the aggregated repertoire that is specific to epitopes associated with 6 viruses (using lists of specific sequences in VD-Jdb), and compared it to OLGA's prediction, P func gen . The agreement was again excellent (triangles in Fig. 6). Again we observe that most epitope-specific sequence groups have large enough frequencies to be found in any individual. Thus, the model can be used to predict the size of repertoire subsets specific to any epitope, as long as specificity data are available for this epitope.

G. Generation probability of sequence motifs
OLGA can also compute the generation probability of any sequence motif, encoded by a string of multiple choices of amino acids. We apply this feature to calculate the net frequency of epitope-specific motifs, and of  motifs that define the TRA sequence of invariant T-cells.
T-cell sequences that can bind a given epitope are often closely related to each other, and this similarity can sometimes be partially captured by sequence motifs. We evaluated the probabilities of motifs derived from a recent study of CDR3 sequence specificity to a variety of epitopes [2]. We took two motifs corresponding to TRA and TRB VJ-CDR3 combinations of TCRs that are known to bind the Epstein-Barr virus HLA-A*0201-BMLF 1280 (BMLF) and the influenza virus HLA-A*0201-M 158 (M1) epitopes. The motifs and generation probabilities are reported in Table I. As a second application, we estimated the probabilities of generating a TRA chain corresponding to one of the motifs associated with Mucosal associated invariant T cells (MAIT) and invariant natural killer T cells (iNKT).
The motifs, which were collected from Gherardin et al. [13], and their probabilities are shown in Table II. The relatively high values for these motifs imply that these invariant chains are generated with high frequency in the primary repertoire and shared by all individuals, confirming the conclusions of Venturi et al. [42].

IV. DISCUSSION
Because the composition of the immune repertoire results from a stochastic process, the frequency with which distinct T-and B-cell receptors are generated is a quantity of primary interest. This frequency is computationally difficult to evaluate because each amino acid sequence can be created by a very large number of recombination events. Our tool overcomes that challenge with dynamic programming, allowing it to process ∼ 50 sequences per second on a single CPU. In its current state OLGA can compute the probabilities of CDR3 sequences and motifs, with or without V/J restriction, of 4 chain loci (human and mouse TRB, human TRA, and human IGH), but the list can readily be expanded by learning recombination models for other loci and species using IGoR [25] which shares the same model format. Obvious additions include the light chains of BCR [40], and more mouse models. While the algorithm evaluates the probability of single chains, recent analyses show that chain pairing in TCR is close to independent [5,15]. The probability of generating a whole TCR receptor can thus be computed by taking the product over the two chains.
OLGA can be used to compute baseline receptor frequencies and to identify outlying sequences in repertoire sequencing datasets. In Elhanati et al. [8] we used it to shed light on the question of public repertoirescomposed of sequences shared by many individualsand predict quantitatively its origin by convergent recombination [23,24,41]. Deviations from the baseline expectancy have been used to identify disease-associated TCR from cohorts of patients [9,10,12,34,48], and to identify clusters of reactive TCRs from tetramer experiments [14] and vaccination studies [30]. Such estimates could be made faster and more reliable by OLGA, especially for rare sequences, and without the need for a negative control cohort [29]. In the future, OLGA could be useful in vaccine and therapy design by focusing atten-tion on clonotypes that are likely to be present in every individual.
We applied OLGA to an experimental database of TCR responding to a variety of disease-associated epitopes. These selected TCR do not differ in their generation probabilities from those of random TCR found in the blood of healthy donors. However, some viral epitopes bind a much larger fraction of the repertoire than others. This observation has potentially important consequences for vaccine design. Since vaccine epitopes stimulate TCR in a pre-existing repertoire, epitopes targeting receptor sequences that are more likely to be generated will have a higher success rate in a wider range of individuals. OLGA can be used to identify such epitopes by computing their specific repertoire fractions, P func gen . While our examples are restricted to TCR, OLGA can also handle BCR and could be used to compute the generation probabilities of BCR precursors of highly reactive or broadly neutralizing antibodies, and thus guide vaccine design in that case as well. The algorithm does not yet handle hypermutations, and extending it to include them would be a useful development.
Acknowledgements. The work of TM and AMW was supported in part by grant ERCCOG n. 724208. The work of ZS and CC was supported in part by NSF grant PHY-1607612. The work of CC was also supported in part by NSF grant PHY-1734030. The work of YE was supported by a fellowship from the V Foundation. The authors declare no conflicts of interest. Recall that the generative VDJ model is defined as: with P aa gen (a 1 , . . . , a L ) = σ∼a P nt gen (σ 1 , . . . , σ 3L ) = E→σ∼a P rec gen (E). (A2) As described in the main text, the dynamic programming algorithm can be summarized by the summation over the positions x 1 , x 2 , x 3 , and x 4 of the following matrix multiplication: The interpretation of the left (subscript) and right (superscript) indices are detailed in the main text, and schematized in Fig. S1. The sums are performed iteratively using matrix multiplications, as detailed in Fig. ??. As in the main text, the nucleotide indices will often be suppressed along with the implicit dependence on the amino acid sequence (a 1 , . . . , a L ). For a given nucleotide position x j , it will be convenient to refer to the amino acid index, and the position in the codon (from both the left and the right), so we introduce the following (graphically shown in the cartoon below): x j = 3(i j − 1) + u j , and u, so that i j encodes the codon that index x j belongs to, and u j its position (from 1 to 3) within that codon, while u * j denotes the position taken from the right of index x j + 1 within its codon, so that u * j = 2 if u j = 1, u * j = 1 if u j = 2, and u * j = 3 if u j = 3. We now define the explicit forms for each of the matrices (note that we retain the indexing x j from Eq A3): a. Vx 1 Contribution from the templated V genes. V x1 can be a 1x1 or 1x4 matrix depending on u 1 . s V is the sequence of the V germline gene (read 5 to 3 ) from the conserved residue (generally the cysteine C) to the end of the gene. l V is the length of s V . These equations are given in the main text.
Contribution from the non-templated N1 insertions (VD junction). M x1 x2 is defined as the product of transfer matrices, and can be a 1x1, 1x4, 4x1, or 4x4 matrix depending on u 1 and u 2 . The transfer matrices are defined by the summed contributions of the Markov insertion model of all codons consistent with the amino acid a (thus summations are over nucleotides y, y 1 , and y 2 to consider all allowed codons): T a (τ, σ) = (y1,y2,σ)∼a S VD (σ|y 2 )S VD (y 2 |y 1 )S VD (y 1 |τ ) (A5) If i 2 > i 1 : where: Contribution from the templated D genes. D(D) x2 x3 can be a 1x1, 1x4, 4x1, or 4x4 matrix depending on u * 2 and u * 3 . s D is the sequence of the D germline gene (read 5 to 3 ) with length l D .
Contribution from the non-templated N2 insertions (DJ junction). N x3 x4 is defined as the product of transfer matrices, and can be a 1x1, 1x4, 4x1, or 4x4 matrix depending on u * 3 and u * 4 . The transfer matrices are defined by the summed contributions of the Markov insertion model of all codons consistent with the amino acid a (thus summations are over nucleotides y, y 1 , and y 2 to consider all allowed codons): If i 4 > i 3 : where: Contribution from the templated J genes. J (D) x4 can be a 1x1 or 4x1 matrix depending on u * 4 . s J is the sequence of the J germline gene (read 5 to 3 ) and l J gives the length of the sequence up to the conserved residue (generally u=1, u*=2 u=2, u*=1 u=3, u*=3 10 11 12 The model used for VJ recombination is quite similar to the model for VDJ recombination with the main differences being the lack of a D segment and an N2 insertion segment. However, a strong correlation between V and J templates is observed in the TRA chain, so we include a joint V, J distribution to allow for this correlation. Due to this similarity, the algorithm used to compute P gen is very similar. The VJ generative model is: with nucleotide and amino acid P gen s being defined the same as for the VDJ recombination model (Eq A2). The dynamic programing algorithm also has a similar form to Eq A3, and can be summarized as (retaining all notation conventions from before): N2 insertions: D alignment: J alignment: V alignment: Contribution from the templated V genes. where This algorithm is validated in the same manner to the VDJ algorithm, i.e. comparing to Monte Carlo (MC) estimation (Fig S3).

Appendix C: Dependence on model parameters and structure
In order to efficiently compute the summation in Eq. A3 the summations of the model contributions from each of the 5 segments of a CDR3 (V genomic, N1 insertions, D genomic, N2 insertions, and J genomic) are performed in a specific order (summarized in Fig S2). Specifically, we start at the left and right ends of the CDR3 read and move inwards, summing over positional indices at each step. As the D and J segments are correlated, it is useful to consider the V and N1 contributions separately from the D, N2, and J and to do the final summation over the index x 2 after the D, N2, and J components are summed over all D alleles (notice the D dependencies in Fig S2). This breakdown is useful to highlight the most computationally intensive steps: N2 insertions and the D alignment. These steps (along with the N1 insertions) require considering that the associated segment could begin and end at each allowed position. This is mathematically seen as the summation over positions and computing a matrix indexed by two indices, leading to an O(L 2 ) complexity. The N2 insertions and D alignments are further aggravated due to model correlations between the D and J genes requiring repeating the steps for N2 insertions and D alignment for each D allele. The runtime of OLGA is thus most sensitive to the maximum number of N2 insertions and the length and number of the D alleles. The effects of varying these parameters is best illustrated by comparing runtimes for mouse TRB, human TRB, and human IGH models (Table S1). In a similar fashion, the most computationally intensive step of computing P gen of a VJ model (e.g. human TRA) is the insertion step, and due to correlations between the V and J genes this is repeated for each J allele in a similar fashion as the D alleles However, as the J region of a human TRA is fairly large, many of these J genes can be excluded from alignment (if they contribute 0 probability), yielding the much faster computation rate of 184 seqs/CPU second.  In order to analyze OLGA's computational performance as a function of CDR3 length, and to compare to other hypothetical methods, we use the human TRB model as an example.
As discussed in the previous section, the most computationally intensive steps of OLGA (N1, N2, and D) require at most O(L 2 ) operations. In practice, OLGA's scaling of the computation speed as a function of CDR3 length, even for the worst case sequences, i.e. fully ambiguous amino acids of a given length, is closer to linear in the relevant regime due to the finite parameterization of the model (maximum number of insertions, maximum size of D sequences, etc). This is shown in Fig S4A. We also compare OLGA to runtimes of IGoR (i.e. direct enumeration of recombination events) and a hypothetical Monte Carlo computation ( Fig S4). As we will explain, neither the IGoR nor the MC are precise comparisons to OLGA, yet OLGA is faster than either.
The IGoR runtimes are for nucleotide sequences not amino acid sequences. In order for IGoR to compute the P gen of an amino acid sequence, it would need to compute and sum the P gen of each nucleotide sequence that codes for the amino acid sequence. These sequences can be enumerated for extremely short CDR3 lengths, however the number explodes exponentially in CDR3 length. Even for a CDR3 length of 4, by enumerating all nucleotide sequences for an amino acid sequence IGoR computes 0.33 seqs/CPU second compared to the 122 seqs/CPU second for OLGA. For longer CDR3 lengths we approximate how long IGoR would take by computing the average number of CDR3 nucleotide sequences per CDR3 amino acid sequence for a given length. OLGA not only heavily outperforms this exponential blowup, but actually outperforms IGoR when IGoR is computing a single nucleotide sequence of a given amino acid sequence.
The Monte Carlo runtime estimate comes from the setup of estimating the P gen of 100,000 sequences. These P gen would be estimated by simulating enough recombination events such that 66% of CDR3 sequences of a given length would be expected to have at least one count. There is a CDR3 length scaling due to the trend that shorter sequences tend to have higher P gen (Fig S5B). The P gen estimated using this methodology will be extremely noisy (Poisson noise on the expected number of counts) and not even give reliable estimates for many sequences.  It is true that the computation time for MC estimates scale as 1/P gen and not with the number of sequences. Thus, there is a hypothetical number of sequences when MC is faster than OLGA if we are willing to accept noisy estimates and to entirely miss some fraction of the CDR3s. This 'crossing point' number of sequences is plotted in Fig S5A and corresponds to completely unrealistic numbers of sequences, highlighting the fact that OLGA will not only give a more reliable P gen , even for very unlikely sequences, but is also much faster than MC even for short, high P gen , sequences. So, even overlooking the drawbacks and imprecision of MC estimation, for plausible sized datasets OLGA is still dramatically faster than MC.
Appendix E: Generation probability distributions from RNA-derived repertoires The analyses described in the main text were mostly concerned with datasets derived by sequencing the genomic DNA contained in a sample of immune cells to directly obtain sequences of the rearranged TCR genes. Immune repertoires can alternatively be obtained by sequencing the mRNA expressed from the same genes, and many such RNA-based data sets exist. Given a TCR sequence, OLGA evaluates the probability of the primitive recombination event (or events) that must have occurred to create the initial T cell carrying that sequence, and the applicability of OLGA is independent of how the sequence was obtained (i.e. from DNA or RNA sequencing). OLGA relies on the availability of a suitable recombination model but that model is thought to vary very little with time (and disease status) for each individual subject and only moderately from individual to individual in a given species. The probability that a given sequence, once generated in a primitive event, will be captured in a sequencing experiment is at best roughly constant across sequences, and may vary substantially between different capture protocols.
For these reasons, it is interesting to investigate how these generation probability distributions vary across CDR3 repertoires obtained using different sequencing protocols in different biological contexts. In Fig. S6 we plot the results of running OLGA on a few recently published human TRB repertoires that were obtained using RNA sequencing. These samples comprise a study of patients with glioblastoma disease (Sims et al. [38]), a study of patients with Crohn's disease and ulcerative colitis (Wu et al. [47]), and a comprehensive study of the dynamics of TCRs in healthy individuals (Wang et al. [44]). Fig. S6 shows the generation probability distribution of data sets from these three sources, for comparison plotted together with the distribution obtained from DNA sequencing of the large human sample of Emerson et al. [9]. As can be seen, two of the three RNA data sets give results quite consistent with the DNA-based results. The glioblastoma data (Sims et al. [38]) gives a distribution broadly similar to the other three, but with a systematic shift to higher frequency of occurrence of lower generation probability sequences. We do not know whether or not this difference is biologically significant, or an artifact of the used protocol. The difference does not seem to be due to sampling depth, as can be seen in Fig. S7, where multiple samples from Sims et al. [38] are plotted: the distributions derived from smaller samples are noisier than, but statistically consistent with, the distributions based on the largest samples.
Appendix F: Cross-species Pgen TCR sequence repertoires are different in detail between species, both because the genomic templates differ and because of differences in the parameters of the recombination process itself. As a result, there are clear interspecies differences in CDR3 length distribution and amino acid composition. Nevertheless, the TRB CDR3 regions of different vertebrate species have the same overall structure and the same conserved residues at the two ends of the CDR3. As a result, a CDR3 from one species usually has a non-zero probability to be produced within a different species, a fact of some interest in the context of studies of cross-species sharing of T cell types. We explored this concept with OLGA by feeding TRB CDR3s produced by the human generation model to a mouse generation model and vice versa. The resulting generation probability distributions are plotted in Fig. S8.
The sequences that are produced in one species have substantially lower probability of being generated in the other species (Fig. S8). The effect is strongest for finding human sequences in a mouse repertoire (compare the black dashed curve with red solid curve): the bulk of the human sequences have extremely low generation probabilities in the mouse model. The effect is less strong for finding mouse sequences in a human repertoire (compare the red dashed curve with the black solid curve): a small fraction of the mouse sequences have generation probabilities that are as high as the highest generation probabilities of human sequences. The results are not symmetric -while the mouse TRBs processed using the human model have a distinct bi-model distribution, the human TRBs have a very flat and low mouse generation probabilities. Furthermore, mouse TRB sequences always have a non-zero probability of being generated in a human TRB context, however 27.4% of human TRB sequences have P gen = 0 as defined by a mouse TRB model. This asymmetry is primarily due to differences in the insertion profiles (humans may have many more , compared to DNA RepSeq data from Emerson et al. [9] (black curve with standard deviation), using a model inferred from [9]. Color indicate sample size: larger datasets are blue, while small ones are red. The only effect of decreasing the sample size is increasing the noise, but the shape stays the same. All TRB datasets from the study are plotted.
inserted N1 and N2 nucleotides) and by extension CDR3 length. Nonetheless, these results suggest that there will be a non-negligible amount of sharing, entirely due to chance statistics, of CDR3 sequences between mouse and human repertoires. A more detailed view of this structure can be seen in a scatter plot of the generation probabilities between the two (Fig. S9). While there are many sequences with high generation probabilities in both the actual generative model and the cross species model, the cross species generation probabilities are much more variable and span many orders of magnitude, without much correlation to the correct species model.

Appendix G: Generation probability distributions from additional pathogen response datasets
In the main text, we displayed the distribution of generative probabilities for T cells known to respond to various pathogens, and even specific epitopes of particular pathogens. The T cell sequences are taken from databases that compile results from multiple experiments. We found that these distributions were, within statistical noise, indistinguishable from the background P gen distribution of PBMCs drawn from the blood. In other words, it would seem that there is no correlation between ease of generation of a T cell and its likelihood to respond to a particular pathogen or epitope. A defect of this analysis is that the database agglomerates sequences from different experimental protocols, so that there is no way of knowing what biases might have affected the inclusion of any given sequence in the database. Obviously, it would be better to do a single well-controlled experiment in which T cells from a single donor are stimulated to expand by selected pathogens, and the expanded T cells sequenced. Such an experiment was reported by Becattini et al. [1] several years ago. In their experiment, CD4+ helper T cells were separated from peripheral blood samples, autologous monocytes from the same samples were incubated with three different pathogens (a fungus, a bacterium, and a toxin) in order to load pathogen epitopes, and helper T cell subsamples (typically containing several million T cells, and hundreds of thousands of clonotypes) were incubated with the prepared monocytes (this was done independently for samples from several donors). The T cells in the various samples that had proliferated under this treatment were separated out (typically yielding millions of cells) and their TRB sequences obtained using the Adaptive Biotechnology genomic DNA protocol. The result is a collection of lists of clones (defined by CDR3 amino acid sequence) from the blood of individual donors that can be said to have expanded under stimulation by the three different pathogens. The responses obtained in this way are quite polyclonal, with a few thousand clonotypes in each list of responding clones (the polyclonality perhaps being due to the fact that stimulation is with preparations of whole pathogens, as opposed to particular pathogen peptides). The P gen distributions of the pathogen-responsive clones for different individuals and pathogens are plotted in Fig. S10 P gen distribution derived from blood samples of healthy individuals, which is also plotted (along with its two-sigma variance across a population of individuals) for reference. These data further strengthen the conclusion that pathogen response activity is uncorrelated with P gen . . For comparison, the background distribution from human peripheral blood TRB sequences from Emerson et al. [9] (with its two sigma variation across multiple individuals) is also plotted. The plotted curves are averages over data from individual donors. The sizes of the responsive T cell repertoires are quite variable: the CA dataset has 39934 clonotypes from 5 donors, the TT dataset has 26573 clonotypes from 4 donors, and the MT dataset has 5082 clonotypes from 2 donors. The generation model was infered from Emerson et al. [9].