Neural Network Models of TCR and BCR Selection.

Previous work has inferred biophysically informed models of V(D)J recombination underlying the generation of TCRs and BCRs (3, 30). In brief, these models are parametrized according to the probabilities by which different V, D, J genes are used and base pairs are inserted in or deleted from the CDR3 junctions to generate a receptor sequence. We infer the parameters of these models using the IGoR software (3) from unproductive receptor sequences, which are generated but due to a frameshift or insertion of stop codons, are not expressed and hence, are not subject to functional selection. The inferred models are used to characterize the generation probability of a receptor sequence $$ P_{\text{gen}}$$ and to synthetically generate an ensemble of preselection receptors (5). These generated receptors define a baseline $$ \mathcal{G}$$ for statistics of repertoires prior to any functional selection.

The amino acid sequence of an immune receptor protein determines its function. To identify sequence properties that are linked to function, we compare the statistics of sequence features $$ f$$ (e.g., V-, J-gene usage and CDR3 amino acid composition) present in a given B or T cell functional repertoire with the expected baseline of receptor generation (Fig. 1C). To do so, we encode a receptor sequence $$ \sigma $$ as a binary vector $$ x$$ whose elements $$ x_f\in \left\{\mathrm{0,1}\right\}$$ specify whether the feature $$ f$$ is present in a sequence $$ \sigma $$. The probability $$ P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}^{\theta }\left(x\right)$$ for a given receptor $$ x$$ to belong to a functional repertoire is described by modulating the receptor’s generation probability $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}\left(x\right)$$ by a selection factor $$ {\mathcal{Q}}^{\theta }\left(x\right)$$,

In soNNia, we divide the sequence features $$ f$$ into three categories: 1) (V, J) usage; 2) CDR3 length; and 3) CDR3 amino acid composition encoded by a $$ 20\times 50$$ binary matrix that specifies the identity of an amino acid and its relative position within a 25-amino acid range from both the 5′ and 3′ ends of the CDR3, equivalent to the left–right encoding of the SONIA model (4) (Fig. 1D). Inputs from each of the three categories are first propagated through their own network. Outputs from these three networks are then combined and transformed through a dense layer. This choice of architecture reduces the number of parameters in the DNN and makes the contributions of the three categories (which have different dimensions) comparable; SI Appendix, SI Text and Figs. S1, S3, and S4 have details on the architecture of the DNN.

The baseline ensemble $$ \mathcal{G}$$, which we have described as being generated from the $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ model (Fig. 1C), can in principle be replaced by any dataset, including empirical ones, at no additional computational cost, for selection inference with soNNia. This is especially useful when the goal is to only compare the selection models associated with different subrepertoires with distinct functions. We will use this feature of soNNia to learn selection coefficients of subsets relative to an empirically constructed generic functional repertoire. In that case, the inferred selection factors $$ Q$$ only reflect differential selection relative to the generic baseline. Importantly, this approach enables us to infer differential selection without having to infer a common underlying generation model $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ for the subrepertoires. After two soNNia models have been learned from two distinct datasets, their statistics may be compared by computing a sequence-dependence log-likelihood ratio $$ r\left(x\right)=\log Q_1\left(x\right)/Q_2\left(x\right)$$ predicting the preference of a sequence for a subset over the other. This log-likelihood ratio can be used as a functional classifier for receptor repertoires (Fig. 1E).

Deep Nonlinear Selection Model Best Describes Functional TCR Repertoire.

First, we systematically compare the accuracy of the (nonlinear) soNNia model with linear SONIA (4) (Fig. 1B) by inferring selection on TCR$$ \beta $$ repertoires from a large cohort of 743 individuals from ref. 33. Our goal is to characterize selection on functional receptors irrespective of their phenotype. To avoid biases caused by expansions of particular receptors in different individuals, we pool the unique nucleotide sequences of receptors from all individuals and construct a universal donor. Multiplicity of an amino acid sequence in this universal donor indicates the number of independent recombination events that have led to that receptor (in different individuals or in the same individual by convergent recombination).

We randomly split the pooled dataset into a training and a test set of equal sizes and trained both a SONIA and a soNNia selection model on the training set (Methods and Fig. 1C). Our inference is highly stable, and the selection models are reproducible when trained on subsets of the training data (SI Appendix, SI Text and Fig. S2).

To assess the performance of our selection models, we compared their inferred probabilities $$ P_{\text{post}}\left(x\right)$$ with the observed frequencies of the receptor sequences $$ P_{data}\left(x\right)$$ in the test set (Fig. 2 A and B). Prediction accuracy can be quantified through the Pearson correlation between the two log frequencies or through their Kullback–Leibler divergence $$ {\mathcal{D}}_{\mathrm{K}\mathrm{L}}\left(P_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}|P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}\right)$$ (Methods and Fig. 2C). A smaller Kullback–Leibler divergence indicates a higher accuracy of the inferred model in predicting the data. The estimated accuracy of an inferred model is limited by the correlation between the test and the training set, which provides a lower bound on the Kullback–Leibler divergence $$ {\mathcal{D}}_{KL}\simeq 0.4$$ bits and an upper bound on the Pearson correlation $$ {\rho }^2\simeq 0.8$$.

Fig. 2.

Performance of selection models on TCR repertoires. Scatterplot of observed frequency, $$ P_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}$$, vs. predicted probability $$ P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$$ for (A) linear SONIA and (B) DNN soNNia models trained on the TCR$$ \beta $$ repertoires of 743 individuals from ref. 33. The baseline is formed by sampling $$ 10^7$$ sequences from the $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ model, learned from the nonproductive sequences of the same dataset (SI Appendix). Color indicates number of sequences. (C) The soNNia model performs significantly better, as quantified by both the Kullback–Leibler divergence $$ {\mathcal{D}}_{\mathrm{K}\mathrm{L}}$$ (Methods) and the Pearson correlation coefficient $$ {\rho }^2$$, without overfitting (SI Appendix, Fig. S2).

We observe a substantial improvement of selection inference for the generalized selection model soNNia with $$ {\mathcal{D}}_{KL}\simeq 1.0$$ bits (and Pearson correlation $$ {\rho }^2\simeq 0.61$$) compared with the linear SONIA model with $$ {\mathcal{D}}_{KL}\simeq 1.6$$ bits (and Pearson correlation $$ {\rho }^2\simeq 0.48$$) (Fig. 2). Both models show a strong effect of selection, reducing the $$ {\mathcal{D}}_{KL}$$ from $$ 3.03$$ bits (and increasing the correlation $$ {\rho }^2$$ from $$ 0.43$$) for the comparison of data with the $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ model alone (Fig. 2). This result highlights the role of complex nonlinear selection factors acting on receptor features that shape a functional T cell repertoire. The features that are still inaccessible to the soNNia selection factors are likely due to the sampling of rare features, individual history of pathogenic exposures, or differences in human leukocyte antigen (HLA) complexes among individuals.

Intra- and Interchain Interactions in TCRs and BCRs.

TCRs are disulfide-linked membrane-bound proteins made of variable $$ \alpha $$ and $$ \beta $$ chains and expressed as part of a complex that interact with pathogens. Similarly, BCRs and antibodies are made up of a heavy and two major groups ($$ \kappa $$ and $$ \lambda $$) of light chains. Previous work has identified low but consistent correlations between features of $$ \alpha \beta $$-chain pairs in TCRs, with the largest contributions between $$ V_{\alpha },V_{\beta }$$ and $$ J_{\alpha },V_{\beta }$$ (34353637–38). In B cells, preferences for receptor features within immunoglobulin heavy (IgH) and light ($$ \lambda $$ or $$ \kappa $$) chains have been studied separately (39, 40), but interchain correlations have not been systematically investigated.

We first aimed to quantify dependencies between chains by reanalyzing recently published single-cell datasets: TCR$$ \alpha \beta $$ pairs of unfractionated repertoires from ref. 37 and BCR of naive cells from ref. 41 (Methods). The blue bars of Fig. 3 show the mutual information between the V and J choices and CDR3 length of each chain, for TCR$$ \alpha \beta $$ (Fig. 3A), IgH$$ \lambda $$ (Fig. 3B), and IgH$$ \kappa $$ (Fig. 3C) repertoires. Mutual information is a nonparametric measure of correlation between pairs of variables (SI Appendix).

Fig. 3.

Inference of selection on intra- and interchain receptor features. Mutual information between pairs of major intra- and interchain features ($$ V$$ and $$ J$$ gene choice and $$ L=$$ CDR3 length for each chain) for (A) TCR$$ \alpha \beta $$, (B) IgH$$ \lambda $$, and (C) IgH$$ \kappa $$ paired chains are shown. Mutual information is estimated directly from data (blue) and from receptors generated based on inferred models: generative baseline (green), linear SONIA (pink), deep-indep (red), and deep joint (yellow). For both TCRs and BCRs, only the deep-joint model (yellow), which correlates the features of both chains through a DNN, is able to recover interchain correlations. Mutual information is corrected for finite-size bias, and error bars are obtained by subsampling (SI Appendix). The diversity of the paired-chain B and T cell repertoires and the contributions of different features to this diversity are reported in SI Appendix, Table S1.

Both TCRs and BCRs have intra- and interchain correlations of sequence features, with stronger empirical mutual dependencies present within chains (Fig. 3).

To account for these dependencies between chains, we generalize the selection model of Eq. 1 to pairs, $$ x=\left(x^a,x^b\right)$$, where $$ \left(a,b\right)=\left(\alpha ,\beta \right)$$ in TCRs or $$ \left(\mathrm{H},\kappa \right)$$ or $$ \left(\mathrm{H},\lambda \right)$$ in BCRs:

Analogously to single chains, we first define a linear selection model specified by $$ Q\left(x\right)=\exp \left({\sum }_f{\theta }_f{x}_f\right)$$, where the sum now runs over features of both chains $$ a$$ and $$ b$$. Because of its multiplicative form, selection can then be decomposed as the product of selection factors for each chain: $$ Q\left(x\right)=Q^a\left(x^a\right)Q^b\left(x^b\right)$$, where $$ Q^a$$ and $$ Q^b$$ are linear models. We also define a deep independent model (deep-indep), which has the multiplicative form $$ Q\left(x\right)=Q^a\left(x^a\right)Q^b\left(x^b\right)$$ but where $$ Q^a$$ and $$ Q^b$$ are each described by DNNs that can account for complex correlations between features of the same chain, similar to the single-chain case (SI Appendix, Fig. S3). The resulting postselection distributions for both the linear and the deep-indep model factorize, $$ P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}\left(x\right)=P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}^a\left(x^a\right)P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}^b\left(x^b\right)$$, making the two chains independent. Thus, by construction neither the linear nor the deep-indep model can account for correlations between chains. Finally, we define a full soNNia model (deep joint) where $$ Q\left(x\right)$$ is a neural network combining and correlating the features of both chains (SI Appendix, Fig. S4).

We trained these three classes of models on each of the TCR$$ \alpha \beta $$ and BCR H$$ \kappa $$ and H$$ \lambda $$ paired repertoire data described earlier. We then used these models to generate synthetic data with a depth similar to the real data and calculated mutual information between pairs of features (Fig. 3). The preselection generation model [$$ Q\left(x\right)=1$$; green bars] explains part but not all of the intrachain feature dependencies, for both T and B cells, while the linear (purple), deep-indep (red), and deep-joint (yellow) models explain them very well (Fig. 3). By construction, the generation, linear, and deep-indep models do not allow for interchain correlations. Only the deep-joint model (yellow) is able to recover part of the interchain dependencies observed in the data. It even overestimates some correlations in BCRs, specifically between the CDR3 length distributions of the two chains and between the heavy-chain $$ J$$ and the light-chain CDR3 length. Thus, the deep structure of soNNia recapitulates both intrachain and interchain dependencies of feature-forming immune receptors.

The inferred selection on correlated interchain receptor features is consistent with previous analyses in TCRs (34353637–38) and is likely due to the synergy of the two chains interacting with self-antigens presented during thymic development for TCRs and preperipheral selection (including central tolerance) for BCRs or later when recognizing antigens in the periphery. Notably, the largest interchain dependencies and synergistic selection are associated with the V-gene usages of the two chains (Fig. 3), which encode a significant portion of antigen-engaging regions in both TCRs and BCRs.

Our results show that the process of selection in BCRs is restrictive, in agreement with previous findings (7), significantly increasing interchain feature correlations. Notably, the increase in correlations (difference between green and other bars) due to selection is larger in naive B cells than in unsorted (memory and naive) T cells. However, the selection strengths inferred by our models should not be directly compared with estimates of the percentage of cells passing preperipheral selection: $$ \sim 10\%$$ for B cells vs. $$ 3\text{to}5\%$$ for T cells (6). Our models identify features under selection without making reference to the number of cells carrying these features. Since the T cell pool in our analysis is a mixture of naive and memory cells, we can expect stronger selection pressures in the T cell data than in the purely naive T cells. However, previous work analyzing naive and memory TCRs separately using linear selection models did not report substantial differences between the two subsets (31).

Lastly, to quantify the diversity of immune receptor repertoires, we compared the entropy of unpaired- and paired-chain repertoires in SI Appendix, Table S1 (SI Appendix, SI Text). These entropy measures suggest a repertoire size (i.e., a typical number of amino acid sequences) of about $$ 10^9$$ receptors for TCR$$ \beta $$ (consistent with ref. 4), $$ 10^7$$ receptors for TCR$$ \alpha $$, $$ 10^{13}$$ receptors for BCR heavy-chain, and $$ 10^4$$ receptors for BCR light-chain sequences. The paired-chain entropy measures suggest repertoire sizes of $$ 10^{16}$$ for TCR$$ \alpha \beta $$ and $$ 10^{17}$$ BCR IgH$$ \lambda $$ and IgH$$ \kappa $$ receptors, which are compatible with the small correlations observed between heavy and light chains in Fig. 3 and previously reported in refs. 34353637–38.

Cell Type- and Tissue-Specific Selection on T Cells.

During maturation in the thymus, T cells differentiate into two major cell types: cytotoxic (CD8⁺) and helper (CD4⁺) T cells. CD8⁺ cells bind peptides presented on MHC class I molecules that are expressed by all cells, whereas CD4⁺ cells bind peptides presented on MHC class II molecules, which are only expressed on specialized antigen presenting cells. Differences in sequence features of CD8⁺ and CD4⁺ T cells should reflect the distinct recognition targets of these receptors. Although these differences have already been investigated in refs. 36 and 42, we still lack an understanding as how selection contributes to the differences between CD8⁺ and CD4⁺ TCRs. In addition to functional differentiation at the cell-type level, T cells also migrate and reside in different tissues, where they encounter different environments and are prone to infections by different pathogens. As a result, we expect to detect tissue-specific TCR preferences that reflect tissue-specific T cell signatures.

To characterize differential sequence features of TCRs between cell types in different tissues, we pool unique TCRs from nine individuals (from ref. 42) sorted into three cell types (CD4⁺ conventional T cells [Tconvs], CD4⁺ Tregs, and CD8⁺ T cells) and harvested from three tissues (pancreatic draining lymph nodes [pLNs], “irrelevant” nonpancreatic draining lymph nodes [iLNs], and spleen).

Training a nonlinear soNNia model (Fig. 1C) for each subset leads to overfitting issues due to limited data. To solve this problem, we train the model in two steps. First, we use the unfractionated data from ref. 43 to construct a shared baseline for all repertoire subsets. We then learn independent linear SONIA models for each repertoire subset so that the inferred $$ Q$$ factors only reflect selection relative to the baseline. We approach this problem in two ways. 1) We infer an SONIA model atop an empirical baseline set $$ \mathcal{G}$$ constructed from the unfractionated repertoire, and 2) we use the technique of transfer learning, which consists of learning a shared nonlinear soNNia model for the unfractionated repertoire and then adding an additional linear layer (similar to standard SONIA) for each subrepertoire (SI Appendix, SI Text and Fig. S5). The subrepertoire selection factors $$ Q$$ inferred by these two approaches are very similar (SI Appendix, Fig. S5), but the former method is simpler, and we use it for our main analysis in Fig. 4. For comparison, we also used the generation model $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ (trained earlier for Fig. 2) as a baseline, in which case the selection factors include selection effects that are shared among the subrepertoires. Distributions of selection factors obtained by both approaches are shown in SI Appendix, Fig. S6.

Fig. 4.

Cell type- and tissue-specific selection on TCRs. (A) Jensen–Shannon divergences ($$ D_{JS}$$) (Eq. 8) computed from models trained on different subrepertoires are shown. (B) Difference in the marginal probability for amino acid composition along the CDR3, $$ P_{\text{post}}^{\text{CD}8}\left(a\right)-P_{\text{post}}^{\text{CD}4}\left(a\right)$$, between CD8⁺ and CD4⁺ Tconv (Left) and the mean difference in the corresponding log-selection factors for amino acid usage $$ \mathrm{\Delta }\log Q=\log Q^{\text{CD}8}-\log Q^{\text{CD}4}$$ (Right) are shown (the mean is taken over the distribution $$ \left(P_{\text{post}}^{\text{CD}8}+P_{\text{post}}^{\text{CD}4}\right)/2$$). The negatively charged amino acids (Aspartate, D, and Glutamate, E) and the positively charged amino acids (Lysine, K, and Arginine, R) are indicated in red and blue, respectively. Other amino acids are shown in gray. (C) Maximum likelihood inference of the fraction of CD8⁺ TCRs in mixed repertoires of conventional CD4⁺ T cells (Tconvs) and CD8⁺ cells from spleen (Eq. 4) is shown. Each repertoire comprises $$ 5\times 10^3$$ unique TCRs. (D) Same as C but for a mixture of Tconv and Treg TCRs. (E) Mean squared error of the inferred sample fraction from C as a function of sample size $$ N$$, averaged over all fractions, using models of increasing complexity: “$$ Q_{VJL}$$” is a linear model with only features for CDR3 length and VJ usage, “linear” is linear SONIA model, and “deep” is the full soNNia model (Fig. 1C). (F) ROC for classifying individual sequences coming from CD8⁺ cells or from CD4⁺ Tconvs from spleen, using the log-likelihood ratios. Curves are generated by varying the threshold in Eq. 5. The accuracy of the classifier is compared with a traditional logistic classifier inferred on the same set of features as our selection models. The training set for the logistic classifier has $$ N=3\times 10^5$$ Tconv CD4⁺ and $$ N=8.7\times 10^4$$ CD8⁺ TCRs, and the test set has $$ N=2\times 10^4$$ CD4⁺ and $$ N=2\times 10^4$$ CD8⁺ TCR sequences.

To quantify differential selection on subrepertoires, we use Jensen–Shannon divergence $$ D_{\mathrm{J}\mathrm{S}}$$ between the distributions of receptors $$ P_{\text{post}}^r$$ and $$ P_{\text{post}}^{r^{\prime }}$$ for pairs of subrepertoires ($$ r,\hspace{0.17em}{r}^{\prime }$$) (Methods). Clustering of cell types based on Jensen–Shannon divergence shows strong differential selection preferences between the CD4⁺ and CD8⁺ receptors, with an average $$ D_{\mathrm{J}\mathrm{S}}\simeq 0.08\pm 0.01$$ bits across respective tissues and subrepertoires (Fig. 4A; SI Appendix, Fig. S7A has similar results where $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ is used as baseline). We identify differential selection between Tconv and Treg receptors within CD4⁺ cells with $$ D_{JS}\simeq 0.015\pm 0.004$$. We also detect moderate tissue specificity for CD8⁺ and Treg receptors, but no such signal can be detected for CD4⁺ Tconv cells across different tissues.

Examining the linear selection factors of the SONIA model trained atop $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ as a baseline reveals the VJ (SI Appendix, Fig. S8) and amino acid usage features (Fig. 4C) that are differentially selected in the Tconv CD4⁺ and CD8⁺ subsets (in spleen). Linear selection models are organized according to a hierarchy from the least to the most constrained model. As one adds selection factors for each feature, the Kullback–Leibler divergence between the repertoire and the baseline increases (SI Appendix). Decomposing in this way the divergence between CD4⁺ Tconv and CD8⁺ repertoires, we find that contributions to the total divergence are evenly split between amino acid features and VJ gene usage, with only a minor contribution from CDR3 length (SI Appendix, Fig. S9). It should be noted that the baseline models $$ P_{\text{gen}}$$ for these subrepertoires, inferred from their unproductive receptors, are similar (SI Appendix, Fig. S10) and do not contribute to these differential preferences.

One key difference between CD4⁺ and CD8⁺ TCRs amino acid composition is their CDR3 charge preferences. We observe an overrepresentation of positively charged (Lysine, K, and Arginine, R) and suppression of negatively charged (Aspartate, D, and Glutamate, E) amino acids in CD4⁺ TCRs compared with CD8⁺ TCRs (Fig. 4B), consistent with previous observations (44). These charge preferences arise due to differential selection on the two subtypes (Fig. 4B), indicating broad differences between amino acid features of peptide–MHC-I and peptide–MHC-II complexes, which respectively interact with CD8⁺ and CD4⁺ TCRs. For example, a statistical survey of peptides presented by different MHC classes shows that MHC-I molecules tend to present more positively charged peptides compared with MHC-II molecules—a bias that is complementary to the charge preferences of the respective TCR subtypes (45).

Decomposing Unsorted Repertoires Using Selection Models.

Knowing $$ P_{\text{post}}^r$$ models specific to subrepertoires enables us to infer the fraction of each class $$ r$$ in unsorted data. Estimating the relative fraction of CD4⁺ and CD8⁺ subtypes in a repertoire can be informative for clinical purposes (e.g., as a probe for tumor-infiltrating lymphocytes), where overabundance of CD8⁺ cells in the sample has been associated with positive prognosis in ovarian cancer (46). Given a repertoire composed of the mixture of two subrepertoires $$ r$$ and $$ r^{\prime }$$ in unknown proportions, we maximize the log-likelihood function $$ L\left(f\right)$$ based on our selection models to find the fraction $$ f$$ of a subrepertoire $$ r$$ within the mixture:

Previous work has used differential V- and J-gene usage and CDR3 length to characterize the relative fraction of CD4⁺ and CD8⁺ cells in an unfractionated repertoire (47). The log-likelihood function in Eq. 4 provides a principled approach for inferring cell-type composition using selection factors that capture the differential receptor features of each subrepertoire, including but not limited to their V and J usage and CDR3 length and amino acid preferences.

To test the accuracy of our method, we formed a synthetic mixture of previously sorted CD4⁺ [Tconv from spleen (42)] and CD8⁺ [from spleen (42)] receptors with different proportions and show that our selection-based inference can accurately recover the relative fraction of CD8⁺ in the mix (Fig 4C). Our method can also infer the proportion of Treg cells in a mixture of Tconv and Treg CD4⁺ cells from spleen (Fig. 4D), which is a much harder task since these subsets are very similar (Fig. 4A). The accuracy of the inference depends on the size of the unfractionated data, with a mean expected error that falls below 1% for datasets with size $$ 10^4$$ or larger for the CD8⁺/CD4⁺ mixture (red and orange lines in Fig. 4E).

Our method uses a theoretically grounded maximum likelihood approach, which includes all of the features captured by the soNNia model. Nonetheless, a simple linear selection model with only V- and J-gene usage and CDR3 length information (blue line in Fig. 4E), analogous to the method used in ref. 47, reliably infers the composition of the mixture repertoire. Additional information about amino acid usage in the linear SONIA model results in moderate but significant improvement (orange line in Fig. 4E). The accuracy of the inference is insensitive to the choice of the baseline model for receptor repertoires: Using the empirical baseline from ref. 43 (Fig. 4E) or $$ P_{\mathrm{g}\mathrm{e}\mathrm{n}}$$ (SI Appendix, Fig. S7D) does not substantially change the results.

The method can be extended to the decomposition of three or more subrepertoires. To illustrate this, we inferred the fractions of Tconv, Treg, and CD8⁺ cells in synthetic unfractionated repertoires from spleen, showing an accuracy of $$ 3\pm 1\%$$ in reconstructing all three fractions (SI Appendix, Fig. S11) in a mixture of size $$ 5\times 10^3$$.

Computational Sorting of CD4⁺ and CD8⁺ TCRs.

Selection models are powerful in characterizing the broad statistical differences between distinct functional subsets of immune receptors, including the CD4⁺ and CD8⁺ TCRs (Fig. 4A). A more difficult task, which we call computational sorting, is to classify individual receptors into functional classes based on their sequence features. In other words, how accurately can one classify a given receptor as a member of a functional subset (e.g., CD4⁺ or CD8⁺ TCRs)?

We use selection models inferred for distinct subrepertoires $$ r$$ and $$ r^{\prime }$$ to estimate a log-likelihood ratio $$ R\left(x\right)$$ for a given receptor $$ x$$ to belong to either of the subrepertoires,

For comparison, we also used a common approach for categorical classification and optimized a linear logistic classifier that takes receptor features (similar to the selection model) as input and classifies receptors into CD8⁺ or CD4⁺ cells. The model predicts the probability that sequence $$ x$$ belongs to subrepertoire $$ r$$ (rather than $$ r^{\prime }$$) as $$ ŷ\left(\mathbf{x}\right)=\zeta \left(R_{\mathrm{l}\mathrm{o}\mathrm{g}}\left(x\right)\right)$$, with $$ R_{\mathrm{l}\mathrm{o}\mathrm{g}}\left(x\right)={\sum }_f{w}_f{x}_f+b$$ and $$ \zeta \left(x\right)=e^x/\left(1+e^x\right)$$. We learn the model parameters $$ w_f$$ and $$ b$$ by maximizing the log likelihood of the training set:

It should be emphasized that despite comparable performances, our fully linear selection-based method provides a biologically interpretable basis for subtype classification, in contrast to black box approaches (36). For example, the relative importance of different sequence features (i.e., CDR3 length, V/J gene identity, and amino acid composition) for CD4⁺ vs. CD8⁺ classification is shown in SI Appendix, Fig. S9.

Classification of TCRs Targeting Distinct Antigenic Epitopes.

Recognition of a pathogenic epitope by a TCR is mediated through molecular interactions between the two proteins. The strength of this interaction depends on the complementarity of a TCR against an antigen presented by an MHC molecule on the T cell surface. Recent growth of data on paired TCRs and their target epitopes (26, 48) has led to the development of machine learning methods for TCR–epitope mapping (25262728–29). A TCR–epitope map is a classification problem that determines whether a TCR binds to a specific epitope. We use our selection-based classifier (Eq. 5) to address this problem. We determine the target ensemble $$ P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}^r$$ from the training set of TCRs associated with a given epitope (positive data) and the alternative ensemble $$ P_{\text{post}}^{r^{\prime }}$$ from a set of generic unfractionated TCRs (negative data). For comparison, we also perform the classification task using the linear logistic regression approach (Eq. 6) and the state-of-the-art TCRex algorithm (28), which uses a random forest model for classification.

We performed classification for the following CD8⁺-specific epitopes, presented on HLA-A*02 molecules: 1) the influenza GILGFVFTL epitope (with $$ N=\mathrm{3,107}$$ associated TCRs), 2) the cytomegalovirus (CMV) NLVPMVATV epitope ($$ N=\mathrm{4,812}$$), and 3) the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) YLQPRTFLL epitope ($$ N=315$$). The first two epitopes have the most abundant associated TCR sets in VDJdb (26, 48), and the latter is relevant for the ongoing COVID-19 pandemic. For consistency with TCRex (28), we used the pooled data from ref. 43 as the negative set and used 10 times more negative data than positive data for training. To quantify performance of each classifier, we performed a fivefold cross-validation procedure. Due to the scarcity of data, we limit our selection inference to the linear SONIA model (Fig. 1C). The ROC curves show comparable performances for the three classification methods on the three epitope-specific TCR sets (Fig. 5 A–C).

Fig. 5.

Selection-based prediction of epitope specificity for TCR. TCRs are classified based on their reactivity to three pathogenic epitopes (columns), using three classification methods: TCRex, log-likelihood ratio (Eq. 5), and linear logistic regression (Eq. 6). (A–C) ROC curves and (D–F) precision-recall curves for (A and D) influenza epitope GILGFVFTL ($$ N=\mathrm{3,107}$$ TCR), (B and E) CMV epitope NLVPMVATV ($$ N=\mathrm{4,812}$$), and (C and F) SARS-CoV-2 epitope YLQPRTFLL ($$ N=315$$) are shown. (G–I) Comparison between log-likelihood scores $$ R\left(x\right)$$ and logistic regression scores $$ R_{\mathrm{l}\mathrm{o}\mathrm{g}}\left(x\right)$$ for the three epitopes. Red points are TCRs that bind the specific epitope (positive set), and black points are TCRs from bulk sequencing (negative set). $$ r$$ is Pearson’s correlation. For all panels, we used pooled data from ref. 43 as the negative set. We used 10 times more negative data than positive data for training. Performance was quantified using fivefold cross-validation.

The TCR–epitope mapping is a highly unbalanced classification problem, where reactive receptors against a specific epitope comprise a very small fraction of the repertoire [less than $$ 10^{-5}$$ (6)]. Precision-recall curves are best suited to evaluate the performance of classification for imbalanced problems. In this case, a classifier should show a large precision (fraction of true predicted positives among all predicted positives) for a broad range of recall or sensitivity (fraction of true predicted positives among positives = true positives + false negatives). The precision-recall curves in Fig. 5 D–F show that TCRex and the logistic classifier can equally well classify the data and moderately outperform the selection-based classifier. While both the logistic classifier and TCRex are optimized for classification tasks, the selection-based classifier is a generative model trained to infer the receptor distribution of interest (positive set) and identify its distinguishing features from the baseline (negative set). As a result, selection-based classification underperforms in the low-data regime, for which fitting a reliable distribution is difficult (e.g., for the SARS-CoV-2 epitope model, with only $$ N=315$$ positive examples). By contrast, the logistic classifier finds a hyperplane that best separates the two sets and therefore, is better suited for classification tasks, and it may be trained on smaller datasets. Nonetheless, we see a strong correlation between the selection-based log-likelihood ratio $$ R\left(x\right)$$ (Eq. 5) and the estimator of the logistic classifier $$ ŷ$$ (Eq. 6), shown for positive set (red points) and the negative set (black points) in Fig. 5 G–I for the three epitopes. This result indicates that the separation hyperplane identified by the logistic classifier aligns well along the effective coordinates of selection that represent sequence features relevant for function in each epitope class.

Data Description.

In this work, we used different datasets to evaluate selection on TCR and BCR features.

1)To quantify the accuracy of the soNNia model (Fig. 2), we used the TCR$$ \beta $$ repertoires from a large cohort of 743 individuals from ref. 33. We pool the unique nucleotide sequences of receptors from all individuals and construct a universal donor totaling $$ 9\times 10^7$$ sequences. We randomly split the pooled dataset into a training set and a test set of equal sizes. We then subsampled the training set to $$ 10^7$$ to reduce the computational cost of inference.

2)To characterize selection on paired-chain receptors (Fig. 3), we analyzed TCR$$ \alpha \beta $$ pairs of unfractionated repertoires from ref. 37 (totaling $$ 5\times 10^5$$ receptors) and BCR of naive cells from ref. 41 (totaling $$ 22\times 10^3$$ and $$ 28\times 10^3$$ receptors for the H$$ \lambda $$ and H$$ \kappa $$ repertoires, respectively).

3)To characterize differential sequence features of TCRs between cell types in different tissues (Fig. 4), we pooled unique TCRs from nine healthy individuals from ref. 42, sorted into CD4⁺ Tconvs, CD4⁺ Tregs, and CD8⁺ T cells, harvested from three tissues: pLNs ($$ 2.3\times 10^5$$ Tconvs, $$ 2.9\times 10^5$$ Tregs, $$ 2.5\times 10^5$$ CD8s), iLNs ($$ 2.0\times 10^5$$ Tconvs, $$ 9.0\times 10^4$$ Tregs, $$ 1.0\times 10^5$$ CD8s), and spleen ($$ 3.2\times 10^5$$ Tconvs, $$ 1.1\times 10^5$$ Tregs, $$ 1.1\times 10^5$$ CD8s). We used the unfractionated data from ref. 43, composed of $$ 2.2\times 10^6$$ receptor sequences, to construct a baseline model for this analysis.

Quantifying Accuracy of Selection Models.

To assess the performance of our selection models, we compare their inferred probabilities $$ P_{\text{post}}\left(x\right)$$ with the observed frequencies of the receptor sequences $$ P_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}\left(x\right)$$ in the test set. Prediction accuracy can be quantified through the Pearson correlation between the two log frequencies or the Kullback–Leibler divergence between the data and the distribution predicted by the selection model $$ P_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$$:

Comparing Selection on Different Subrepertoires.

To characterize differences in subrepertoires due to selection, we evaluate the Jensen–Shannon divergence $$ D_{\mathrm{J}\mathrm{S}}\left(r,r^{\prime }\right)$$ between the distribution of pairs $$ \left(r,r^{\prime }\right)$$ of subrepertoires $$ P_{\text{post}}^r$$ and $$ P_{\text{post}}^{r^{\prime }}$$,