text:
Zachary Sethna

Giulio Isacchini

Thomas Dupic

Thierry Mora

Aleksandra M Walczak

Yuval Elhanati

Population variability in the generation and thymic selection of T-cell repertoires

The diversity of T-cell receptor (TCR) repertoires is achieved by a combination of two intrinsically stochastic steps: random receptor generation by VDJ recombination, and selection based on the recognition of random self-peptides presented on the major histocompatibility complex. These processes lead to a large receptor variability within and between individuals. However, the characterization of the variability is hampered by the limited size of the sampled repertoires. We introduce a new software tool SONIA to facilitate inference of individual-specific computational models for the generation and selection of the TCR beta chain (TRB) from sequenced repertoires of 651 individuals, separating and quantifying the variability of the two processes of generation and selection in the population. We find not only that most of the variability is driven by the VDJ generation process, but there is a large degree of consistency between individuals with the inter-individual variance of repertoires being about ∼2% of the intra-individual variance. Known viral-specific TCRs follow the same generation and selection statistics as all TCRs.

I. INTRODUCTION

Most organisms live in a similar environment, facing common pathogenic threats. However, the adaptive immune system, based on the stochastic VDJ recombination process, is a naturally diverse system, supporting both repertoire variability within the individual, and variability across the population [START_REF] Jr | Immunobiology: The Immune System in Health and Disease[END_REF]. Quantifying both types of variability, and understanding how they support a robust immune response, are still open questions. Determining the variability under normal healthy conditions is a crucial step for understanding the immune system in compromised situations such as infections, autoimmune diseases, and cancer. The adaptive immune system reacts specifically against a variety of different threats to the organism. This is achieved by maintaining a large ensemble of T cells, each having a different receptor that binds distinct subsets of antigens. The adaptive immune system maintains this diversity by generating a large repertoire of cells with different receptors [START_REF] Hozumi | Evidence for somatic rearrangement of immunoglobulin genes coding for variable and constant regions[END_REF][START_REF] Venturi | Sharing of T cell receptors in antigen-specific responses is driven by convergent recombination[END_REF][START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF] and then selecting them according to their binding properties. The first step of selection occurs in the thymus. Cells carrying receptors that bind too strongly or too weakly to the host's own proteins do not pass this selection [START_REF] Kyewski | A central role for central tolerance[END_REF][START_REF] Starr | Positive and negative selection of t cells[END_REF]. The remaining cells are let out into the periphery and undergo selection for binding of foreign antigen which results in cell prolif-eration. In all cases, T cell receptors (TCR) bind to antigen fragments presented as short peptides on the major histocomptability complex (MHC) of presenting cells [START_REF] Clambey | Molecules in medicine mini review: the αβ t cell receptor[END_REF]. Each human individual has 6 types of MHC molecules encoded by the very polymorphic human leukocyte antigen (HLA) locus. All of these processes -receptor generation, selection, and peptide presentation -are stochastic in nature and depend on the host's genetic background.

High-throughput T cell repertoire sequencing (RepSeq) provides a census of the T cell repertoire found in a blood or tissue sample [START_REF] Heather | Highthroughput sequencing of the T-cell receptor repertoire: pitfalls and opportunities[END_REF][START_REF] Lindau | Advances and applications of immune receptor sequencing in systems immunology[END_REF][START_REF] Six | The past, present and future of immune repertoire biology -the rise of next-generation repertoire analysis[END_REF][START_REF] Woodsworth | Sequence analysis of T-cell repertoires in health and disease[END_REF]. These samples are generally indicative of the true repertoire, and comparing them over a population yields similarities predicted based on MHCs, pathogenic history and general properties of the generation process [START_REF] Dewitt | Human t cell receptor occurrence patterns encode immune history, genetic background, and receptor specificity[END_REF][START_REF] Elhanati | Predicting the spectrum of tcr repertoire sharing with a data-driven model of recombination[END_REF]. Due to the large diversity of possible TCRs, different samples, even ones taken from the same individual under the same conditions, will often differ substantially due to statistical noise. As a result, characterization of a repertoire sample is often more reliably done by statistically modeling the underlying generation and selection processes instead of working with raw TCR sequences and read counts. In this paper we take such an approach to characterize the diversity of the human T-cell receptor beta chain (TRB) repertoire. This approach allows us to disentangle the two processes of generation and selection, and to quantify their relative contribution to the overall variability across individuals. Our results provide a quantification of natural TCR diversity which is essential for studying adaptive immunity in clinical contexts. [START_REF] Emerson | Immunosequencing identifies signatures of cytomegalovirus exposure history and HLAmediated effects on the T cell repertoire[END_REF]. For each person i = A, B, C, ... we define a personalized TRB generation model P i gen , and a personalized thymic selection model Q i (σ), as both processes are expected to vary across individuals as a function of their genetic background, in particular their HLA type. The generation model allows us to evaluate the probability of generating each receptor sequence σ in each individual i. Q i (σ) tells us how likely a given receptor amino-acid sequence σ is to pass thymic selection in a given individual. Combined together, the two models give the probability of a given TRB amino acid sequence in the repertoire of a given person

II. RESULTS

A. Data source and modeling strategy

P i post (σ) = Q i (σ)P i gen (σ). (B)
To learn these models, for each individual we separated sequences into productive and nonproductive sequences. Nonproductive sequences are free of selection effects and were used to learn the generation model, Pgen, using the IGoR software [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]. Most productive sequences are subject to selection and were used to learn the selection model, Q, by matching the statistics of the data with those of sequences generated synthetically with Pgen (using the OLGA software [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF]) and weighted by Q. Once the model is learned, the probabilities of amino-acid TRB sequences preand post-selection can be calculated using OLGA and SONIA.

We analyzed previously published RepSeq data from a large cohort study [START_REF] Emerson | Immunosequencing identifies signatures of cytomegalovirus exposure history and HLAmediated effects on the T cell repertoire[END_REF] consisting of TRB nucleotide sequences from blood samples of 651 healthy individuals. Sample sizes ranged from 50,000 to 400,000 unique CDR3 amino acid beta chains. For each individual i, we learned an individual-specific generation model, which describes the probability of generating a given aminoacid sequence σ by VDJ recombination, P i gen (σ), and an individual-specific selection model, Q i (σ), defined as the fitness of each sequence upon thymic selection. The resulting probability distribution of receptor sequences is P i post (σ) = Q i (σ)P i gen (σ) (Fig. 1A). To learn these models, TRB nucleotide sequences were divided into productive and non-productive sequences, where productive sequences are defined as being in frame with no stop codon. The pipeline is summarized in Fig. 1B. We applied the IGoR algorithm [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF] to non-productive TRBs of each individual to learn P i gen (σ). Productive sequences were used to learn individual-specific selection models Q i (σ) as in Ref. [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF], by comparing them with simulated productive sequences generated from the individual specific generation models. A new software package, SONIA, was developed to perform the Q i (σ) inference. For each sequence in each individual the algorithm computes two probabilities: its generation probability, P i gen , and its post-selection probability in the periphery, P i post . We then use them to estimate the intra-and inter-person variability.

B. Individual variability of VDJ recombination statistics

The model of VDJ recombination, P gen , assigns a probability to each VDJ recombination scenario [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF], where a scenario is a particular choice of the various recombination events: germline gene choice (V, D, and J), the number of deletions to those germline genes at the V-D and D-J junctions, and the number and identities of the untemplated, inserted nucleotides at each of the junctions (called N1 for the V-D junction, and N2 for the D-J junction). A detailed description of the model is given in the Methods section. Each recombination scenario determines a particular nucleotide sequence. The generation probability of a sequence is then the sum of all recombination scenarios that result in that sequence. Since the scenario is a hidden variable of the observed nucleotide sequence, we can use the Expectation Maximization algorithm to infer the maximum-likelihood estimator of the model parameters [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF] using IGoR [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF].

Productive sequences then translate into an aminoacid sequences σ, and we denote by P gen (σ) the probability of generation of σ conditioned on it being productive, equal to the sum of the generation probabilities of all possible nucleotide variants divided by the probability to generate a productive sequence (an abuse of notation relative to the strict definition of P gen with no conditioning on being productive).

We find that the generation models learned from dif- ferent individuals in our cohort, P i gen , are consistently similar to each other, with more variation in the gene usage than in the junctional diversity statistics (Fig. 2). The distributions of the number of inserted N1 and N2 nucleotides vary little (Fig. 2A). The biases of the untemplated inserted nucleotides, governed by a Markov model where the choice of each inserted base pair depends stochastically on the previous insertion [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF], is also conserved across individuals (Fig. 2B). Note that these probabilities are also similar for the N1 and N2 insertions provided that N2 is read in the anti-sense. Likewise, gene specific deletion profiles have very low variability (Fig. 2E). By contrast, gene usage shows greater yet moderate inter-individual variability (Fig. 2C-D). Overall, these results confirm the large level of reproducibility of the generation process over a large cohort.

Probability

N1 (V-D junction) N2 (D-J junction)

A A A C A G A T C A C C C G C T G A G C G G G T T A T C T G T T 0.0 0.1 0.2 0.3 0.4 0.5 Transition Probability N1 N2 T T T G T C T A G T G G G C G A C T C G C C C A A T A G A C A A T R B J1 -1 T R B J1 -2 T R B J1 -3 T R B J1 -4 T R B J1 -5 T R B J1 -6 T R B J2 -1 T R B J2 -2 T R B J2 -3 T R B J2 -4 T R B J2 -5 T R B J2 -6 T R B J2
TRBJ2-7 T R B V 1 T R B V 1 0 T R B V 1 1 T R B V 1 2 T R B V 1 3 T R B V 1 4 T R B V 1 5 T R B V 1 6 T R B V 1 7 T R B V 1 8 T R B V 1 9 T R B V 2 T R B V 2 0 T R B V 2 3 T R B V 2 4 T R B V 2 5 T R B V 2 6 T R B V 2 7 T R B V 2 8 T R B V 2 9 T R B V 3 T R B V 3 0 T R B V 4 T R B V 5 T R B V 6 T R B V 7 T R B V
We then asked whether these small individual variations in the recombination statistics were correlated as a result of shared biological mechanisms or genetic factors. We found that the numbers of insertions at the two junctions were highly correlated with each other (Pearson's r = 0.79), meaning that individuals that tend to have longer N1 insertions also tend to have longer N2 insertions on average (Fig. 3A). N1 insertions were also slightly longer by ∼ 0.17 insertions on average. The variance of the number of insertions calculated over the repertoire of one individual is extremely correlated to its mean (Pearson's coefficients of 0.88 and 0.87, Fig. 3B), suggesting a single individual-specific parameter controlling both N1 and N2 length distributions. This parameter is likely linked to the activity of the Terminal Deoxynucleotidyl Transferase (TdT) enzyme responsible for Mean N length N insertions [START_REF] Bogue | Regulation of N-region diversity in antigen receptors through thymocyte differentiation and thymus ontogeny[END_REF].

To quantify other correlations we calculated Pearson's correlation coefficient over the population between combinations of various parameters. In order to determine significance and account for the finite cohort size we also compute a 'shuffled' Pearson's coefficient for each parameter combination by scrambling the individuals to destroy correlations. Fig. 3C shows the normalized distribution of Pearson's correlation for the combinations of the marginal distributions of V, D, and J usages. Correlations between V -V , J -J, and V -J marginals all show substantial excess of positive and negative values relative to the shuffled control. Full parameter co-variations are shown in Figs. S1-S3. To determine which types of parameters co-vary the most, we computed the rescaled standard deviation of the Pearson's correlation coefficients of all combinations of parameter types (Fig. 3D). This analysis reveals that V gene usage co-varies with itself, D and J usages are also correlated with each other, as well as N1 length with N2 length, and the insertion biases at N1 and N2 with each other.

C. Learning models of thymic selection with SONIA

After VDJ recombination, new T cells go through an initial selection process in the thymus before being released as naive T cells to the periphery. Positive thymic selection selects for functionally useful receptors, while negative selection removes T cells that recognize selfpeptides to avoid auto-immunity. Thymic selection skews the statistics of the repertoire of TRB sequences in quantifiable ways. This can be seen by comparing the length distribution of the Complementarity Determining Region 3 (CDR3, running from a conserved cysteine near the end of V segment through a conserved phenylalanine near the beginning of the J segment) of productive sequences drawn from the generation model to observed sequences (Fig 4A). We observe a substantial narrowing of the distribution post-selection, eliminating sequences much longer or shorter than 14-15 amino acids [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF].

To characterize these differences more systematically, we use a statistical model of selection to account for differences between the repertoire generated from the ...CASSPNGAWTYEQYF... V J

CDR3, L = 15

Length-Position model:

. . .

Left+Right model:

...CASSPNGAWTYEQYF... ) and left (log 10 q i,right (aa)) displayed as sequence logos for 6 individuals. The first row shows the sequence logos for the amino acid usage from the generation model alone (consistent across individuals), with the usual convention that the total height of the logo is equal to the Shannon entropy of amino acid usage at this position, and the relative height of each letter is proportional to its usage. (F,G) Distributions of selection factors for V and J genes, qV J , over the population (averaged over one of them, as selection factors as defined for the joint usage of V and J).

raw VDJ recombination (pre-selection) and the observed repertoire of productive sequences (post-selection). Since selection acts on the functionality of a receptor we restrict ourselves to productive amino acid sequence statistics. Mathematically, we require that the post-selection distribution, P post = Q(σ)P gen (σ), agrees with the statistics of productive sequences in the frequency of a select set of features, f ∈ F, while remaining as close as possible to P gen (where distance is measured by the Kullback-Leibler divergence, σ P post (σ) ln(P post (σ)/P gen (σ))). This is done by choosing the sequence-specific selection factors Q(σ) which can be shown to take the form (see Methods):

Q(σ) = 1 Z f ∈F (σ) q f , (1)
where F(σ) ⊂ F is the subset of features present in sequence σ. Solving for the factors q f that match the frequencies of features in the data is equivalent to maximum likelihood estimation (MLE).

Features may be the presence of a given amino-acid at a given position, the use of a particular V or J gene, a particular CDR3 length, or any combination thereof. For example, some of the features of the TRB designated by (CASSGRQGVATQYF, TRBV06-05, TRBJ02-05) are 'CDR3 length 14', 'S in position 2 from the left', 'Y in position -2 from the right' and 'V gene is TRBV06-05'.

To facilitate the definition and learning of such selection models, we introduce the software package SONIA. SONIA allows for a flexible definition of model features and infers the selection factors q f using MLE. The input to SONIA is a list of selected amino acid sequences and, if needed, their V and J gene choice. By default SONIA uses P gen as provided by an IGoR inferred model (using OLGA as a generation engine [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF]), but it can also take as an input a custom sample of pre-selection sequences. This can be useful for identifying selection pressures during immune challenges using different choices of pre and post-selection repertoires (see Methods for details).

We applied SONIA using two models corresponding to two choices of feature sets. In the LengthPosition model [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF], features include all possible choices of combinations of V and J genes, all possible CDR3 lengths, as well as amino acids usage at each position and length (Fig. 4B,top). This choice allows for great flexibility at the cost of many parameters. The LengthPosition model replicates the results of Ref. [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF] (Fig. 4C).

The number of parameters can be reduced by noting that selection pressures on amino acids near the 5' (left) or 3' (right) end of the CDR3 appear to depend only on their relative position to that end, regardless of CDR3 length (Fig. 4C). The Left+Right model exploits that regularity by defining features of amino-acid usage at positions relative to the 5' end of the CDR3 (denoted by a positive index), or to its 3' end (denoted by a negative index). This model has much fewer parameters, since features are defined for left and right positions regardless of CDR3 length, and can be written as a special parametrization of the LengthPosition model, in which each amino acid contributes to the selection factor through the product of a left and a right factor (Fig. 4B, bottom).

To evaluate the accuracy of the Left+Right model, we computed its predictions for the frequencies of amino acid usages at each position and length (Fig. 4D, see also Fig. S4 for overall amino-acid usage). These statistics are by construction matched by the LengthPosition model but not necessarily by the Left+Right model, and thus provide a good test of the validity of the parameter reduction it affords. While predictions from VDJ generation model (blue dots) do not reproduce the empirical frequencies well, highlighting the need of a selection model, both the LengthPosition (red dots) and the Left+Right (black dots) models match the data well. As the Left+Right model captures the observed behavior with fewer parameters, we will work with this model for the remainder of the paper.

Fig. 4E displays the selective pressures on the CDR3 amino acid composition (q i,left and q i,right ) from the left and right positions across a choice of 6 (out of 651) individuals, in the form of sequence logos. These selective factors are mostly conserved across individuals. Fig. 4F and Fig. 4G show the selection factors for the V and J genes (q V J ) averaged over one of the two segments. Again, the pattern is mostly concordant across the population, but with some substantial differences for a few genes that have greater variability. Thus, much as in the generation process, individual variability in the selection process is moderate and concentrates on gene usage rather than CDR3 statistics.

D. Population variability

To quantify more precisely the variability of the generation and selection processes across 651 individuals, we computed the distributions of log 10 P gen , log 10 Q, and log 10 P post for each individual (see Methods). Figs. 5A-C show the results as a density map over the entire population, indicating strong consistency between individuals. The distributions over sequences from the model (obtained by sampling from P post using importance sampling, black curve) agree very well with those obtained from the data (red). By contrast, sequences generated from P gen , without selection factors (blue) fail to reproduce the data.

The shift to high Q values from the pre-to the postselection model is present by construction in the distribution of the Q (Fig. 5B), because the post-selection ensemble should be enriched in high selection factors. However, a similar shift to higher probabilities from pre-to postselection is indicative of a correlation between the generation probability, P gen , and the selection factor, Q (Table I, Fig. S5). This correlation suggests that evolution has shaped VDJ recombination to favor sequences that are likely to pass thymic selection, as previously argued [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF]. Fig. 5D summarizes the distributions of probabilities P in different probability ensembles of decreasing diversity: raw VDJ recombination scenarios (black), generated nucleotide sequences (green), pre-selection productive amino acid sequences (blue, same as the blue curves in Fig. 5A), and post-selection productive amino acid sequences (red, the mean of which is the black curve in Fig. 5C). The negative of the mean of log 10 (P ) is, up to a ln(2)/ ln [START_REF] Six | The past, present and future of immune repertoire biology -the rise of next-generation repertoire analysis[END_REF] factor, equal to the Shannon entropy of the distribution expressed in bits, -log 2 (P ) P . Fig. 5E shows the distribution of these entropies across the pop-ulation. The width of the distributions of log 10 P is strongly correlated with their means across individuals and also from pre-selection to post-selection (Fig. 5F), suggesting again a single parameter driving individual variability, possibly the average number of N insertions.

We also plot the P post distributions of TRBs from the VDJdb database that are known to be specific to human viruses [START_REF] Bagaev | VDJdb in 2019: database extension, new analysis infrastructure and a T-cell receptor motif compendium[END_REF] (Fig. 6). There does not appear to be a substantial shift in the post-selection probability of these viral-specific sequences as compared to productive TRBs from blood. A similar absence of bias was previ- ously reported for the distribution of generation probabilities [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF], suggesting that the VDJ recombination process is not explicitly skewed towards generating these viralspecific sequences. Our results further show that thymic selection also does not seem to be biased to select for sequences specific to these viral epitopes.

E. Quantifying overall variability and its contribution due to generation and selection

The overall variability in the TRB repertoire can be characterized both between and within individuals in the population, by calculating the variance of the distribution of log 10 P post , which gives a measure of the typical fold-variation. Since log 10 P post (σ) = log 10 P gen (σ) + log 10 Q(σ), this variance can be decomposed as:

Var(log 10 P post ) =Var(log 10 P gen ) + Var(log 10 Q) + 2Cov(log 10 P gen , log 10 Q).

(

) 2
To quantify the range of repertoire variability within an individual we calculate the variances and covariance of log 10 P i gen , log 10 Q i , log 10 P i post over the data sequences, and synthetic sequences for each individual. Table I summarizes the average of these variances over the 651 individuals. 80% of the variation comes from the generation process, with the remainder mostly stemming from a strong correlation between selection and generation, as previously discussed (Fig. 5A andC, SI Fig. 5).

Variations in the probabilities of given sequences across individuals (averaged over sequences, see Methods for details) are much lower (Table II), highlighting the high level of consistency in the population. The total variance of 0.091 in log 10 P post corresponds to relative variations of 10 ± √ 0.091 ∈ (-50%, +100%) in the probability of sequences. While those differences are substantial in absolute terms, they are 1.6% of the variance over sequences within an individual (≈ 5.4, see Table I). Much of this variance again stems from VDJ generation.

To further characterize variability, we learned 'consensus' or 'universal' models from sequences sampled randomly from each individual. To this end we inferred a consensus VDJ generation model (P univ gen ) from outof-frame sequences, and a consensus Left+Right SONIA model (Q univ and P univ post ) from the productive sequences (Methods). We then compared each individual model to the universal model using the Jensen-Shannon divergence, an information-theoretic measure of distance between probability distributions expressed in bits and directly comparable to entropies (Methods). The distributions of JSD(P i gen , P univ gen ) and JSD(P i post , P univ post ) over the cohort highlight the consistency of these models with most individuals having < 0.3 bits JSD from both P univ gen and P univ post (Fig. 7). This should be compared to the associated entropies of > 30 bits for either distribution (Fig. 5E).

III. DISCUSSION

By applying distinct computational procedures to the nonproductive and productive sequences of the TCR repertoires of a large cohort of 651 donors, we were able to learn individual-specific models of repertoires, separating the processes of generation and thymic selection. This allowed us to quantify precisely the variability of each process within the population.

We found that the TRB generation process varied only moderately between individuals, with two main drivers: gene usage and average length of untemplated insertions. Because insertions contribute a lot to the generation probability, the latter is the main driver of variability in the distribution of P gen itself. V,D, and J gene usage variability may be due to variations in the regulatory signals, both genetic and epigenetic, that control the operation of the Recombination-Activating Gene (RAG) protein that initiates the recombination process [START_REF] Schatz | V(D)J recombination: mechanisms of initiation[END_REF]. It may also be due to variations in gene copy numbers of the gene segment, as was observed in the related case of the IgH locus [START_REF] Luo | Estimating Copy Number and Allelic Variation at the Immunoglobulin Heavy Chain Locus Using Short Reads[END_REF]. Variations in the mean number of insertions could be attributed to differences in expression of TdT as well as other proteins involved in the non-homologous end joining pathway [START_REF] Lieber | The Mechanism of Double-Strand DNA Break Repair by the Nonhomologous DNA End-Joining Pathway[END_REF].

We found that the inferred selection models were also variable between individuals, but the magnitude of these variations remains limited, which may be surprising considering that different individual's repertoires are subject to different selective pressures due to diverse HLA backgrounds. Overall, the ratio of the total variability across individuals to the variability across sequences within an individual (measured by the variance of the logarithm of the sequence probability) was only 1.6%, of which about 85% came from variations in the generation process, and 15% from the selection process. We also found that sequences that were previously identified to be specific to human viruses did not differ in their generation or selection probability from generic sequences from blood, finding no evidence in our models for an evolutionary mechanism to favor such viral-specific sequences (as suggested in [START_REF] Thomas | Selected before selection: A case for inherent antigen bias in the T cell receptor repertoire[END_REF]), neither in the process of VDJ recombination, nor through thymic selection.

Thymic selection of naive T cells was found to be well captured by a model where selection acts independently on each amino acid, regardless of the sequence context. The variability in the inferred parameters for the selection models is not large in the population, identifying reproducible features in different individuals. This suggests that the main statistical effects of thymic selection captured by our model are mostly universal, probably driven by positive selection for amino acids that makes a folding functional receptor. The effect of HLA specific positive and negative selection, on the other hand, might not be well captured by this kind of a model, which focuses on finding broad sequence features rather than specific sequences to harness more statistical power, although variations in the V and J selection factors may reflect HLA types. Our approach thus complements the strategy of looking for associations of particular TCRs with HLA type, which was previously applied to the same dataset [START_REF] Dewitt | Human t cell receptor occurrence patterns encode immune history, genetic background, and receptor specificity[END_REF]. An obvious limitation of this and other studies of that dataset is that it comprises a restricted subset of the human population.

While in this study we used SONIA for the purpose of comparing peripheral to pre-selection repertoires, the software is written to be flexible in several ways. First, if can be used to infer selection factors between any two repertoires (observed or generated), by inferring selection factors that match the statistics of the two samples. Second, SONIA can go beyond selection pressures on single amino acids, allowing features of pairs or motifs of amino acids. Finally, SONIA can be applied to other chains than TRB, notably the alpha chain of the TCR (TRA) as well as immunoglobulin IgH. SONIA's flexibility opens up the possibility of using SONIA to find statistical correlations in various biological or clinical contexts. SONIA could be applied to samples that are known to have responded to some perturbations, for example after vaccination or infection [START_REF] Thomas | Tracking global changes induced in the CD4 T-cell receptor repertoire by immunization with a complex antigen using short stretches of CDR3 protein sequence[END_REF][START_REF] Pogorelyy | Precise tracking of vaccineresponding t cell clones reveals convergent and personalized response in identical twins[END_REF]. In such a context clone sizes may be crucial to identify the underlying changes. To facilitate this, SONIA can also infer selection factors from read-count weighted repertoires. A major challenge in the field of immune repertoire profiling remains to decipher the specificity of the TCR-pMHC interaction. Vaccine design, immunotherapy and therapy for autoimmune conditions would all greatly benefit from the ability to find or design TCRs with known specificity. In the last couple of decades experimental methods have been developed for identifying TCRs specific to given antigens [START_REF] Wolfl | Activation-induced expression of cd137 permits detection, isolation, and expansion of the full repertoire of cd8+ t cells responding to antigen without requiring knowledge of epitope specificities[END_REF][START_REF] Altman | Phenotypic analysis of antigenspecific t lymphocytes[END_REF][START_REF] Glanville | Identifying specificity groups in the T cell receptor repertoire[END_REF][START_REF] Dash | Quantifiable predictive features define epitope-specific T cell receptor repertoires[END_REF]. Based on accumulated TCR binding data [START_REF] Bagaev | VDJdb in 2019: database extension, new analysis infrastructure and a T-cell receptor motif compendium[END_REF], computational methods have been proposed recently that can find clusters of similarly reactive TCRs [START_REF] Pogorelyy | Precise tracking of vaccineresponding t cell clones reveals convergent and personalized response in identical twins[END_REF][START_REF] Glanville | Identifying specificity groups in the T cell receptor repertoire[END_REF][START_REF] Dash | Quantifiable predictive features define epitope-specific T cell receptor repertoires[END_REF][START_REF] Pasetto | Tumor-and Neoantigen-Reactive T-cell Receptors Can Be Identified Based on Their Frequency in Fresh Tumor[END_REF], or to predict TCR specificity to a given epitope using machinelearning techniques [START_REF] Jurtz | NetTCR: sequence-based prediction of TCR binding to peptide-MHC complexes using convolutional neural networks[END_REF][START_REF] Sidhom | DeepTCR: a deep learning framework revealing structural concepts within TCR Repertoire[END_REF][START_REF] Springer | Prediction of specific TCR-peptide binding from large dictionaries of TCR-peptide pairs[END_REF][START_REF] Jokinen | TCRGP : Determining epitope specificity of T cell receptors[END_REF]. SONIA could be used to learn flexible models of these antigen-specific TCR subsets and to study their organization. It could also be applied to identify specific selective pressures in particular subsets, defined by HLA specificity, pathogenic history, clinical status, T-cell phenotype (naive, effector, memory, CD4, CD8, regulatory T cells), or to differentiate distinct samples from the same individual, such as blood, tissue, or tumor samples.

IV. METHODS

A. Data

The data used for the inference of both the VDJ generation models and the subsequent selection models are the Adaptive Biotechnologies sequenced TRB repertoires of Emerson et. al. [START_REF] Emerson | Immunosequencing identifies signatures of cytomegalovirus exposure history and HLAmediated effects on the T cell repertoire[END_REF]. An initial quality control pass was done over the 664 individuals to ensure at least 10,000 unique out of frame sequences to be used to infer the VDJ generation model. 651 individuals passed this threshold and all were used in the subsequent analyses.

All analyses were done on unique nucleotide reads, discarding any cell count information. This is done to ensure that each sequence is reflective of a single recombination event, which is an important restriction when modeling VDJ recombination and thymic selection. For some selection modeling purposes (e.g. modeling antigen cell counts may be incorporated. In practice, amino-acid sequences are reduced to the choice of V and J, and the full amino acid CDR3 sequence.

Sequences were determined to be productive and used in the selection analysis if they had a non-zero P gen . Beyond being an in-frame sequence without stop codons, this requires that a sequences retains the conserved residues defining the CDR3 region (Cysteine on the 5 end, Phenylalanine or Valine on the 3 end) as well as aligning to non-pseudo V and J genes.

B. Generation model

The generation model is defined at the level of the recombination scenarios in order to reflect the underlying biology of VDJ recombination. Each recombination scenario is defined by the gene choice (V , D, and J); deletions/palindromic insertions for each gene (d V , d D , d D , and d J ); and the sequence of non-templated nucleotides at each junction (m 1 , . . . , m V D and n 1 , . . . , n DJ ). The probability of a recombination scenario is given in the factorized form:

P scenario = P V (V )P delV (d V |V )P DJ (D, J) × P delD (d D , d D |D)P delJ (d J |J) × P insVD ( VD )p 0 (m 1 ) V D i=2 S VD (m i |m i-1 ) × P insDJ ( DJ )q 0 (n DJ ) DJ -1 i=1 S DJ (n i |n i+1 ) . (3)
This model factorization, originally from Murugan et al, has been shown to capture the relevant correlations between the different recombination events in TRB [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF].

The probability of a nucleotide sequence x is given by:

P nt gen (x) = scenario→x P scenario , (4)
and the probability of a productive amino-acid sequence is:

P gen (σ) = 1 F x→σ P nt gen (x), (5)
where F = scenario|prod P scenario is the total probability that a random recombination event is productive (inno stop codons, preserves conserved residues, and does not use pseudo-genes as germline gene choices). F can be computed directly from a generative model using OLGA [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF].

C. Selection model

To minimize the Kullback-Leibler distance between P post and P gen while enforcing the constraints σ:f ∈F (σ) P post (σ) ≡ P post (f ) = P data (f ) for each f , we extremize the following Lagrangian:

σ P post (σ)   ln P post (σ) P gen (σ) - f ∈F (σ) λ f -µ   , (6)
where λ f are Lagrange multipliers constraining the frequencies of f , while µ ensures the normalization of P post . This extremization yields the form of P post :

P post (σ) = P gen (σ) exp   f ∈F (σ) λ f   . (7)
Defining q f = e λ f , and Z = e -µ , we obtain Eq. 1. Given that form, the Lagrange multipliers must be adjusted to satisfy the constraints. Doing so is equivalent to maximizing the likelihood of the data under the model:

L = 1 N σ∈data ln P post (σ|{λ f }), ( 8
)
where N is the number of data sequences. This can be shown by noting that the gradient of the log-likelihood,

∂L ∂λ f = 1 N   σ∈data:f ∈F (σ) 1   -P post (f ) = P data (f ) -P post (f ), (9)
cancels when the constraints are satisfied.

D. SONIA implementation

SONIA is a python software built to define and infer feature-defined selection models. SONIA has built in procedures for defining and identifying sequence features of CDR3 sequences. SONIA also ships with the prepackaged selection models of LengthPosition and Left+Right features. With a feature model defined, SONIA takes as an input a list of productive amino acid CDR3s, along with any aligned V/J genes. This list of observed CDR3s can be either reduced to unique sequences (useful when learning thymic selection and the background statistics are based on unique sequences) or sequences taken with their clonality to account for a non-flat clone size distribution. As an optional input, SONIA can read in baseline CDR3 and aligned V/J genes to use as the background that the selection model is learned from. Alternatively, OLGA's sequence generation machinery [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF] is built into SONIA so a generation model can be specified and background sequences automatically generated.

SONIA has built-in methods to compute the feature marginals over the data sequences, background sequences, and the selection model. These marginals are use to fit the selection model iteratively using Tensor-Flow keras [START_REF] Chollet | Keras[END_REF][START_REF] Abadi | TensorFlow: Large-scale machine learning on heterogeneous systems[END_REF] with the Kullback-Leibler divergence as a loss function. We checked the convergence of the algorithm and its satisfying of the constraints after convergence (Fig. S6) An inferred SONIA model can be used to compute overall selection factors Q of any sequence. In combination with OLGA, SONIA can compute P post and to generate selected sequences through rejection sampling.

E. Distributions of probabilities

We produced the distributions of P gen , Q, and P post shown in Fig. 5A-C by comparing the productive data sequences of each individual to a synthetic sample of productive sequences generated from P i gen of that individual using OLGA [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF]. The number of generated sequences for each individual were matched to the number of productive data sequences. For each dataset, we calculated P i gen using OLGA, and Q i and P i post using SO-NIA's Left+Right model. The Q-weighted curves are determined by weighting each generated sequence by its selection factor Q i and then renormalizing.

For Fig. 5D, we used 300,000 scenarios, nucleotide, and amino acid sequences were generated from each individual's VDJ generation model. Again, we used OLGA to compute the various generation probabilities P i , where P i is P i scenario , P i,nt gen , or P i gen . Entropy was estimated as -log 2 P i over the respective generated sample. For the post-selection ensemble (P i post ), the distributions were weighted by Q i computed by SONIA, and the entropy was calculated as -Q(σ) log 2 [P gen (σ)Q(σ)] over the generated amino acid sequences.

F. Inference and Probability computation

Overall workflow is summarized in Fig. 1B. VDJ generation models were all inferred using IGoR [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]. Amino acid P gen distributions were all computed using OLGA [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF] according to the specified IGoR model parameters. All generated sequences were drawn from the corresponding VDJ generation model using OLGA. Lastly, selection models were all inferred, and evaluated using SONIA. The code for all processes is available on GitHub: IGoR: https://github.com/qmarcou/IGoR OLGA: https://github.com/zsethna/OLGA SONIA: https://github.com/statbiophys/SONIA

G. Quantifying variability

To produce the variances and covariances of Table I we took the productive data sequences from each individual along with an equivalent number of synthetic sequences drawn from the individual's VDJ generation model. For each sequence we computed P univ gen , Q univ , and P univ post using the consensus models. The variance and covariance of each quantity was computed over both the data sequences and generated sequences for each individual. These variances and covariances were then averaged over the individual cohort to yield the numbers in Table I. Error bars are the standard deviation over the cohort.

For Table II, we learned a consensus VDJ generation model P univ gen from nonproductive sequences sampled ran-domly from all individuals. 300,000 productive sequences were drawn from P univ gen to serve as a generated sequence pool. For data sequences we used 326,000 productive sequences sampled randomly from all individuals. We calculated for each sequence σ the individual specific P i gen , Q i , and P i post for each individual, then calculated the variances and covariances over i. Finally we averaged the results over the sequences σ from each pool.

The Jensen-Shannon divergence between two distribution P 1 and P 2 is defined as:

JSD(P 1 , P 2 ) = 1 2 σ P 1 (σ) ln P 1 (σ) P (σ) + P 2 (σ) ln P 2 (σ) P (σ) , (10)
with P = (P 1 + P 2 )/2. V genes (ordered as on chromosome) V genes (ordered as on chromosome) V genes (ordered as on chromosome)

FIG. 1 :

 1 FIG. 1: Analysis pipeline. (A)We analyzed data from Tcell receptor beta (TRB) repertoires of 651 donors collected by Emerson et al.[START_REF] Emerson | Immunosequencing identifies signatures of cytomegalovirus exposure history and HLAmediated effects on the T cell repertoire[END_REF]. For each person i = A, B, C, ... we define a personalized TRB generation model P i gen , and a personalized thymic selection model Q i (σ), as both processes are expected to vary across individuals as a function of their genetic background, in particular their HLA type. The generation model allows us to evaluate the probability of generating each receptor sequence σ in each individual i. Q i (σ) tells us how likely a given receptor amino-acid sequence σ is to pass thymic selection in a given individual. Combined together, the two models give the probability of a given TRB amino acid sequence in the repertoire of a given person P i post (σ) = Q i (σ)P i gen (σ). (B) To learn these models, for each individual we separated sequences into productive and nonproductive sequences. Nonproductive sequences are free of selection effects and were used to learn the generation model, Pgen, using the IGoR software[START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]. Most productive sequences are subject to selection and were used to learn the selection model, Q, by matching the statistics of the data with those of sequences generated synthetically with Pgen (using the OLGA software[START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF]) and weighted by Q. Once the model is learned, the probabilities of amino-acid TRB sequences preand post-selection can be calculated using OLGA and SONIA.

FIG. 2 :

 2 FIG. 2: Distribution of the individual P i gen model parameters over 651 individuals. All plots are violin plots with the mean and standard deviation shown by error bars. (A) Insertion length distributions of the N1 and N2 junctions. (B) Markov transition probabilities for the inserted nucleotide identities at the N1 (red) and N2 (blue) junctions. The N2 transition probabilities are organized in a reverse complementary fashion to the N1 transition probabilities. (C) V gene family usages. (D) Joint D and J gene usages. (E) Deletion profiles for individual J genes.

FIG. 3 :

 3 FIG. 3: Correlations between model parameters across individuals. (A) Mean number of N1 (VD junction) versus N2 (DJ junction) insertions (each point is an individual). (B) Variance (across sequences) versus mean number of insertions at both junctions (each point is an individual). (C) Distribution of Pearson correlation coefficients between any two usage probabilities P (V ) or P (J) across individuals. (D) Rescaled standard deviation of Pearson coefficients of parameter combinations over various recombination events. Values are rescaled by the standard deviation of the shuffled distribution (≈ 0.39 in all cases).

FIG. 4 :

 4 FIG. 4: Thymic selection models of 651 individuals. (A) Length distribution of the Complementarity Determining Region 3 (CDR3) of TRB before (as predicted by the Pgen model, in blue), and after (data, in red) thymic selection. Violin plots show variability across individuals. (B) Schematic of the two SONIA model architectures used in this article. Both models have selection factors for the joint choice of V and J, qV J , and for the CDR3 length L, qL. The LengthPosition model has selection factors defined for each amino acid at each position i and length L, qi,L. The Left+Right model factorizes those factors into two contributions depending on the position of the amino acid from the left and right, respectively. (C) Amino-acid selection factors qi,L of the LengthPosition model as a function of position i and L for each of the 20 amino acids. These factors are consistent with previous reports on a smaller cohort[START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF]. (D) Model prediction for the frequencies of all features of the LengthPosition model (V,J joint usage, CDR3 length, and amino acid usage at each position and length). The Left+Right model reproduces all the probabilities despite not having learned them directly. (E) Model parameters of the Left+Right model, for right (log 10 q i,left (aa)) and left (log 10 q i,right (aa)) displayed as sequence logos for 6 individuals. The first row shows the sequence logos for the amino acid usage from the generation model alone (consistent across individuals), with the usual convention that the total height of the logo is equal to the Shannon entropy of amino acid usage at this position, and the relative height of each letter is proportional to its usage. (F,G) Distributions of selection factors for V and J genes, qV J , over the population (averaged over one of them, as selection factors as defined for the joint usage of V and J).

FIG. 5 :

 5 FIG. 5: (A-C) Distributions of Pgen, Q, and Ppost calculated over many sequences for each individual. Shown are the postselection productive TRBs from each individual (red), and pre-selection sequences generated from the individual's VDJ generation model P i gen (blue). The distributions for all individuals are visualized using a density map indicating the local density of probability distribution curves over the cohort. (D) Density maps of the model distributions for the VDJ recombination scenarios, Pscenario, the nucleotide sequences, P nt gen , the productive amino-acid sequences upon generation, Pgen, and post-selection amino-acid sequences Ppost, over the population. The same convention for the density map is used. Error bars for (A-D) are the standard deviation over the population. (E) Distributions of the Shannon entropies of Pscenario, P nt gen , Pgen, and Ppost over the population. (F) Mean vs variance of log 10 Pgen and log 10 Ppost over both productive and generated sequences for each individual. The linear relation suggests a single parameter explaining variability in the population.

FIG. 6 :FIG. 7 :

 67 FIG. 6: Distribution of TRB sequences from the VDJdb database specific to human viruses [19] compared to the productive sequences from the blood of 651 individuals. A) log 10 (Q univ ) distribution for each individual's productive data sequences (gray heatmap) and for viral-specific TCRs from the VDJdb database. B) log 10 (P univ post ) distributions. The VDJdb log 10 (P univ post ) distributions are Gaussian-smoothed for clarity. Q univ and P univ post are 'universal' models learned from sequences randomly drawn from all individuals.

  FIG. S1: Rescaled Pearson coefficients for length insertion distributions. (A) N1-N1 correlations. (B) N2-N2 correlations. (C) N1-N2 correlations. The N1 and N2 distributions are highly correlated over the 651 individual cohort. Rescaling is done by normalizing by the standard deviation of correlation coefficients obtained by shuffling individuals for the two features independently.

  FIG. S3: (A) RescaledPearson coefficients for V-J correlations across the 651 individual cohort. (B) Rescaled Pearson coefficients for V-V correlations. V genes are ordered by position on the chromosome. While large V-J and V-V correlations exist, no obvious chromosomal structure emerges. Rescaling as in S1.

  FIG.S5: Scatter plots of log 10 (Q univ ) vs log 10 (P univ gen ) for (A) generated sequences drawn from P univ gen and (B) data sequences used to infer log 10 (Q univ ). The color scale indicates the local probability density of the points (on a log scale). This visualizes the correlation of Pgen and Q as described in Tab. I. Q univ and P univ post are 'universal' models learned from sequences randomly drawn from all individuals.

  FIG. S6: Convergence of the universal Left+Right model Q univ . (A) L1 convergence, per learning epoch, of the marginals (or frequencies) between the data features and the model features. (B) Scatter plot of the feature marginals. The x-axis shows the frequencies of features of the data, while the y-axis show the model prediction for the generation model (red) and for Q-weighted Left+Right model (blue). The L1 distance in (A) measures the mean distance between the blue dots and the diagonal.

TABLE I :

 I Intra-individual variation

	Seqs	Quantity	Intra-indiv Var	%
	Gen seqs	Var log 10 Pgen	7.05 ± 0.49	78.2%
		Var log 10 Q	0.407 ± 0.016 4.52%
		Cov (log 10 Pgen, log 10 Q) 0.778 ± 0.069 2×8.64%
		Var log 10 Ppost	9.01 ± 0.63	100%
	Data seqs	Var log 10 Pgen	4.67 ± 0.53	86.9%
		Var log 10 Q	0.185 ± 0.009 3.44%
		Cov (log 10 Pgen, log 10 Q) 0.258 ± 0.033 2×4.81%
		Var log 10 Ppost	5.37 ± 0.56	100%

TABLE II :

 II Inter-individual variation

	Seqs	Quantity	Inter-indiv Var	%
	Gen seqs	Var log 10 Pgen	0.121	95.5%
		Var log 10 Q	0.0175	13.7%
		Cov (log 10 Pgen, log 10 Q)	-0.00586	2× -4.62%
		Var log 10 Ppost	0.127	100%
	Data seqs	Var log 10 Pgen	0.0792	87.0%
		Var log 10 Q	0.0132	14.5%
		Cov (log 10 Pgen, log 10 Q)	-0.000697	2×-0.766%
		Var log 10 Ppost	0.0910	100%

Acknowledgements. The work of TM and AMW was supported in part by grant ERCCOG n. 724208. The authors have no conflicts of interest.


domain: phys sdv

text:
Giulio Isacchini

Zachary Sethna

Yuval Elhanati

Armita Nourmohammad

Aleksandra M Walczak

Thierry Mora

On generative models of T-cell receptor sequences

T-cell receptors (TCR) are key proteins of the adaptive immune system, generated randomly in each individual, whose diversity underlies our ability to recognize infections and malignancies. Modeling the distribution of TCR sequences is of key importance for immunology and medical applications. Here, we compare two inference methods trained on high-throughput sequencing data: a knowledge-guided approach, which accounts for the details of sequence generation, supplemented by a physics-inspired model of selection; and a knowledge-free Variational Auto-Encoder based on deep artificial neural networks. We show that the knowledge-guided model outperforms the deep network approach at predicting TCR probabilities, while being more interpretable, at a lower computational cost.

I. INTRODUCTION

Deep learning methods are proving a very useful approach in many areas of physics and the natural sciences [START_REF] Ching | Opportunities and obstacles for deep learning in biology and medicine[END_REF][START_REF] Carleo | Machine learning and the physical sciences[END_REF]. These algorithms are successful in identifying hidden patterns in large amounts of data, often helping make progress in situations where traditional analyses reach their limits [START_REF] Yamins | Using goal-driven deep learning models to understand sensory cortex[END_REF][START_REF] Senior | Improved protein structure prediction using potentials from deep learning[END_REF][START_REF] Yang | Machine-learningguided directed evolution for protein engineering[END_REF]. Despite the black box aspect of how the algorithm works and the lack of interpretability of the model features, machine learning is undoubtedly useful, especially in cases where the natural system of interest escapes our intuition or knowledge. However, as we show here on the example of immune repertoires, introducing physical or biological intuition into data-driven models can outperform basic uninformed machined learning approaches.

The adaptive immune system is made up of a large ensemble of diverse lymphocyte receptors that recognize different pathogens. The receptors expressed on the surface of T cells (T cell receptors -TCR) are generated by randomly assembling genomic templates for three genes (variable-V, diverse-D and junction-J) that make up of the so-called β chain, and two (V and J) for the α chain. Additionally to this combinatoric diversity, nontemplated nucleotides are added at the junctions between these templates and nucleotides are deleted. Such recombined DNA forms the newly generated TCR that later undergoes thymic selection that tests for its ability to form a receptor protein and bind, albeit not too strongly, proteins that are natural to the host organism [START_REF] Yates | Theories and quantification of thymic selection[END_REF]. TCR that pass thymic selection are released into the periphery and form the naive repertoire (i.e. non-stimulated by foreign antigens). Due to the random addition and deletion * These authors contributed equally. of nucleotides, receptors sequences have different lengths and some are even out of frame or have stop codons, in which case they are called nonproductive. Conversely, sequences with no frameshift nor stop codon are conventionally called productive. High-throughput immune repertoire sequencing experiments sample blood from individual hosts, sort out TCRs, and sequence this subset [START_REF] Heather | Highthroughput sequencing of the T-cell receptor repertoire: pitfalls and opportunities[END_REF][START_REF] Minervina | TCR and BCR repertoire profiling in adaptive immunity[END_REF][START_REF] Bradley | Using T Cell Receptor Repertoires to Understand the Principles of Adaptive Immune Recognition[END_REF]. Analysis of this kind of data makes it possible to characterize the statistics of both generated and naive repertoires.

TCR sequences differ from classical protein families, which are grouped by function and across species [START_REF] El-Gebali | The Pfam protein families database in[END_REF]. Those families are believed to have evolved over long time scales under a shared selective pressure that shapes their statistics. For such families, physics-inspired statistical inference methods have helped to predict contacts between amino-acids in the protein [START_REF] Cocco | Inverse Statistical Physics of Protein Sequences: A Key Issues Review[END_REF], define sectors of co-evolving residues [START_REF] Halabi | Protein Sectors: Evolutionary Units of Three-Dimensional Structure[END_REF], or find interaction partners [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. Deep [START_REF] Riesselman | Deep generative models of genetic variation capture the effects of mutations[END_REF] and non-deep [START_REF] Tubiana | Learning protein constitutive motifs from sequence data[END_REF] machine learning approaches have also been successfully applied. By contrast, TCR generation is fairly well understood mechanistically. Previously we developed a statistical inference technique that uses biological knowledge of the underlying assembly processes to learn the statistics of generation and calculate the generation probability of each TCR sequence [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF][START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]. Since thymic selection involves many specific interactions with antigen-presenting cells, modeling it from first principles is more difficult. Nevertheless, simple models of selection based on the assumption of an additive fitness [START_REF] Berg | Selection of DNA binding sites by regulatory proteins. Statistical-mechanical theory and application to operators and promoters[END_REF] have been shown to well recapitulate some key statistics of these ensembles [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF][START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF]. However, a direct test of the performance of this method for the abundance of specific sequences in large cohorts is still lacking.

Recently, Davidsen et al. [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] described an elegant approach for learning the distribution of T-cell receptor beta sequences (TCRβ or simply TCR in the follow-ing), based on a Variational Auto-Encoder (VAE). The method makes it possible to generate new sequences with the same statistics as real repertoires, and to evaluate the frequency of individual sequences, which agree with the data with good accuracy. Its main strength is that it does not take any information about the origin of these sequences through VDJ recombination and thymic and peripheral selection. Yet it manages to extract statistical regularities imprinted by these processes.

Here we compare the VAE method [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] with the previously proposed model of generation and selection, called SONIA [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF][START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF]. We compare their performances for predicting the distribution of TCR sequences in controlled conditions, training and validating on the same datasets. Contrary to the claims of the original VAE paper [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF], we show that that knowledge guided models perform as well as the variational auto-encoder or even better, at a lower computational cost.

II. MODEL DEFINITIONS A. Knowledge-guided model

To predict the probability distribution of TCR sequences, we build a generative model that proceeds in two steps: initial generation, and selection.

First, a recombination model for the probability of generation of a sequence σ, denoted by P gen (σ), is learned from failed, nonproductive rearrangements, which are free of selection biases [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF][START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]. This model describes in detail the probabilities of V, D, and J usages, and of deletion and insertion profiles. Calling E the collective variable describing the recombination scenario, the model predicts its probability P scenario (E). Its parameters are learned through Expectation-Maximization using the IGoR software [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF].

Although the model is trained on non-productive sequences, it can be used to predict the probability of any sequence. Denoting σ(E) the amino-acid sequence produced by scenario E, we define the generation probability of a productive amino-acid sequence σ as:

P gen (σ) = 1 F E P scenario (E)I[σ(E) = σ], (1)
where I(•) is the indicator function, and F = E P scenario (E)I(σ(E) is productive) is the probability that a random recombination scenario results in a productive sequence. More precisely, σ is defined by the choice of V and J genes (σ V and σ J ), as well as the amino-acid sequence of the Complementarity Determining Region 3 (CDR3) that lies between V and J, σ 1 , . . . , σ L . The sum in (1) involves a large number of terms due to the degeneracy of both the genetic code and the recombination process, but it can be done using a recursive technique akin to transfer matrices implemented in the OLGA software [START_REF] Sethna | OLGA: fast computation of generation probabili-ties of B-and T-cell receptor amino acid sequences and motifs[END_REF].

Second, a model of selection, called SONIA [START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF], is learned on top of the generation probability P gen to describe the distribution of productive sequences,

P SONIA (σ) = Q(σ)P gen (σ), (2)
where

Q(σ) = 1 Z exp h V JL (σ V , σ J , L) + L i=1 h i,L (σ i ) (3)
is a selection factor calculated through additive "fields" h acting on the sequence elements, similarly to additive position-weight matrix models first introduced for DNA binding sites [START_REF] Berg | Selection of DNA binding sites by regulatory proteins. Statistical-mechanical theory and application to operators and promoters[END_REF]. Within this framework, we can define three models according to the parametrization of h. In the first two models, the VJL field is decomposed as

h V JL (σ V , σ J , L) = h V J (σ V , σ J ) + h L (L). A first model in which h i,L is left unconstrained is called the "Length-Position" (LP)
model. This choice corresponds to the original model of [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF], in which the selective pressure on each aminoacid may depend on the sequence length L. However, observations [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF] suggest that these factors are to some extend independent of L. This invariance can be incorporated by assuming that the field can be decomposed into two contributions depending on the position of the amino acid from the right and left ends of the CDR3: h i,L = h i,right + h L-i+1,left . The resulting "Left+Right" (LR) model has much fewer parameters and is less likely to overfit the data. For these two models, parameters are learned by maximizing the log-likelihood with an L 2 regularization using gradient ascent, as specified in Ref. [START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF].

In addition, because no software implementation of the selection model was provided with the original article [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF], Davidsen et al. [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] compared their VAE approach to a reduced version of this selection model (not examined in [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF]), which they call OLGA.Q. In that model, only VJ usage and CDR3 length were included: h i,L = 0. Its parameters h V JL were fitted by maximizing the likelihood analytically.

B. Variational Auto-Encoder

A VAE is an auto-encoder whose structure can be used as a generative probabilistic model. A good introduction can be found in Ref. [START_REF] Kingma | An Introduction to Variational Autoencoders[END_REF]. In short, a VAE consists of a probabilistic encoder q(z|σ) and a probabilistic decoder p(σ|z), converting the sequence into a continuous multidimensional latent variable z and back. The goal of the encoder is to make the probabilistic mapping from σ to itself through q and p as faithful as possible, while at the same time making the distribution of the latent variable z as close as possible to a simple distribution, i.e. multivariate Gaussian with unit covariance.

Both p and q are parametrized by deep neural networks, whose parameters are optimized for these two objectives, using stochastic gradient descent. Once the model is learned, new sequences can be generated by drawing z from p 0 (z), and σ from p(σ|z), so that σ is distributed according to P VAE (σ) = dz p(σ|z)p 0 (z).

In practice, the predicted probability of a given sequence P VAE (σ) is evaluated using Monte-Carlo importance sampling. In Ref. [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF], a variant of the traditional auto-encoder detailed in [START_REF] Higgins | beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework[END_REF] was used. Here we focus on the version of the VAE called basic in that paper.

III. MODEL COMPARISON A. Datasets and model training

The data consists of TCRβ sequence repertoires of 666 individuals [START_REF] Emerson | Immunosequencing identifies signatures of cytomegalovirus exposure history and HLAmediated effects on the T cell repertoire[END_REF]. We use the exact same procedure, dataset, and subsamples as in [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] for reproducibility.

For each individual, read counts are first discarded as they stem from clonal expansions. To train an initial P gen model on which SONIA is built and trained, we used 2 • 10 5 nonproductive sequences drawn randomly from all donors. For all models, unique amino-acid sequences were first separated into a training and a testing dataset of equal sizes. All models were then trained on 2 • 10 5 or 10 6 TCRβ sequences randomly sampled from the training dataset with replacement, according to their frequency in the cohort, counting each unique nucleotide sequence in each patient. Their performance was assessed by their ability to predict the frequency of sequences from the testing set, P data (σ).

B. Predicting sequence frequencies

We used two measures of performance: Pearson's ρ 2 between the logarithms of the frequencies as in [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF], and the Kullback-Leibler divergence: D KL = log 2 [P data (σ)/P model (σ)] (model = VAE or SONIA), where the average • is taken over 10 4 sequences from the testing set, sampled according to their relative frequencies within that set. We excluded ∼ 0.3% of sequences for which P gen = 0, probably due to sequencing errors. Note that, if not for the L 2 regularization, maximizing the loglikelihood would be equivalent to minimizing the D KL . The scale of D KL may be compared to the total entropy of the ensemble,σ P SONIA (σ) log 2 P SONIA (σ) ≈ 31 bits [START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF].

Fig. 1 shows the predicted frequencies of the Left+Right SONIA model and the VAE model, both trained on the same 2 • 10 5 sequences, and compares them to data. The performances of all models and both datasets are reported in Table I. SONIA models perform generally better than the VAE, especially the Left+Right model which is the best model according to both measures of performance. Note that the Length-Position model of Ref. [START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF], also performs as well as the VAE. Davidsen et al. [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] did not compare their model to it owing to the absence of a readily available implementation.

Strikingly, even the basic model of generation with no selection (h = 0), P gen , performs comparably to the VAE, and sometimes better according to the ρ 2 measure, despite the model being trained on nonproductive sequences. Accordingly, the OLGA.Q model, which adds a minimal layer of selection on top of P gen , also performs very well. These results differ substantially from the ρ 2 = 0.26 -0.27 reported in [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] for OLGA.Q. In [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF], the default model for P gen was not actually trained on the dataset of interest, but rather used with its default parameters learned from a different dataset, which explains the poor reported performance.

We can also compare the two models by asking whether the distribution of frequencies are well reproduced by one another, using another TCR dataset from [START_REF] Neuter | Memory CD4+ T cell receptor repertoire data mining as a tool for identifying cytomegalovirus serostatus[END_REF] to allow for a direct comparison to the results of Ref. [START_REF] Davidsen | Deep generative models for T cell receptor protein sequences[END_REF] (Fig. 2). Both the VAE and SONIA agree with the data in their distribution of P model . VAE-generated sequences have the same distribution of P SONIA as SONIA-generated sequences, with a slight under-estimation of the distribution peak, and an excess of low-frequency sequences (Fig. 2a). The converse is true when looking at the distribution of P VAE for SONIA-versus VAE-generated sequences (Fig. 2b). This suggests that the VAE and SO-NIA capture some features of the sequence statistics that are distinct from one another.

C. Computational times

SONIA is an order of magnitude faster than the VAE, which uses Monte-Carlo sampling to calculate predicted frequencies. The average computing time for P SONIA (σ) is 14 ms per sequence on a laptop computer and 3 ms on a 16-core computer, versus 0.18 s for P VAE (σ) (no parallelization possible).

SONIA was also faster to train. It took 33 minutes to train a SONIA model on 10 6 sequences using a 30-core computer, to which one should add 31 minutes to train an IGoR model on 2 • 10 5 nonproductive sequences. For the same amount of data and on the same machine, the VAE took 7 hours to train.

IV. CONCLUSION

In summary, both approaches, VAE and SONIA, perform equally well, with perhaps a slight advantage for the latter. SONIA is also much faster. These results suggest that, while knowledge-free approaches such as the VAE perform well, there is still value in preserving the structure implied by the VDJ recombination process as a baseline for learning complex distributions of immune repertoires. Extending the SONIA model considered here beyond a simple linear combination of features, and taking ideas from the modeling strategy of the VAE, offers interesting directions for future improvement in repertoire modeling.

In a more general context, while machine learning approaches are undoubtably a very useful tool, they can be made even more powerful when combined with models that describe the underlying physics or biology. This is the case when training data is limited, as has been reported in complex image processing of non-animate matter [START_REF] Colas | Nonlinear denoising for characterization of solid friction under low confinement pressure[END_REF]. As we show, even if data is abundant, using models to guide learning can help.

FIG. 1 :FIG. 2 :

 12 FIG. 1: Predicted TCR sequence probabilities (y-axis) versus empirical frequencies (y-axis), for (a) the SONIA Left+Right model (ρ 2 = 0.53) and (b) the VAE model (ρ 2 = 0.47). Models were trained on 2 • 10 5 sequences sampled from the training set assembled from the TCR β repertoires of 666 donors [25]. Frequencies refer to empirical frequencies in the same datasets. The SONIA model was built on top of a Pgen model trained on 2 • 10 5 non-productive sequences from the same donors.

  10 6 ρ 2 10 6 DKL 2 • 10 5 ρ 2 2 • 10 5 DKL

	VAE	0.48	1.7	0.47	2.0
	Pgen	0.48	4.5	0.51	4.5
	OLGA.Q	0.48	2.6	0.47	2.6
	SONIA LP 0.52	1.8	0.52	1.7
	SONIA LR 0.53	1.4	0.53	1.4

TABLE I :

 I Pearson's correlation coefficients ρ 2 and Kullback-Leibler divergence DKL (in bits) for the various models. Either 1 million or 200,000 sequences were used in the training dataset.

Acknowledgements. The work of GI, TM and AMW was supported by grant European Research Council COG 724208.

Code availability. All code for reproducing the figures of this comment can be found at https://github. com/statbiophys/compare_selection_models_2019/. The SONIA package upon which that code builds is available at https://github.com/statbiophys/SONIA/ [20].