# A transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity

Barbara Bravi, Andrea Di Gioacchino, Jorge Fernandez-De-Cossio-Diaz, Aleksandra Walczak, Thierry Mora, Simona Cocco, Rémi Monasson

► To cite this version:

Barbara Bravi, Andrea Di Gioacchino, Jorge Fernandez-De-Cossio-Diaz, Aleksandra Walczak, Thierry Mora, et al.. A transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity. eLife, 2023, 12, pp.e.85126. 10.7554/eLife.85126. hal-04252230

HAL Id: hal-04252230

https://hal.science/hal-04252230

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_1.png" /></figure>

HAL
open science

# A transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity

Barbara Bravi, Andrea Di Gioacchino, Jorge Fernandez-De-Cossio-Diaz, Aleksandra Walczak, Thierry Mora, Simona Cocco, Rémi Monasson

► To cite this version:

Barbara Bravi, Andrea Di Gioacchino, Jorge Fernandez-De-Cossio-Diaz, Aleksandra Walczak, Thierry Mora, et al.. A transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity. eLife, 2023, 12, pp.e.85126. 10.7554/eLife.85126. hal-04252230

HAL Id: hal-04252230

https://hal.science/hal-04252230v1

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Learning the differences: a transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity

Barbara Bravi¹,², Andrea Di Gioacchino²,⁺, Jorge Fernandez-de-Cossio-Diaz²,⁺, Aleksandra M. Walczak²,*, Thierry Mora²,*, Simona Cocco²,*, and Rémi Monasson²,\*

¹Department of Mathematics, Imperial College London, London, SW7 2AZ United Kingdom

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France

⁺Equal Contribution

*Equal Contribution

## Abstract

Antigen immunogenicity and the specificity of binding of T-cell receptors to antigens are key properties underlying effective immune responses. Here we propose diffRBM, an approach based on transfer learning and Restricted Boltzmann Machines, to build sequence-based predictive models of these properties. DiffRBM is designed to learn the distinctive patterns in amino acid composition that, one the one hand, underlie the antigen's probability of triggering a response, and on the other hand the T-cell receptor's ability to bind to a given antigen. We show that the patterns learnt by diffRBM allow us to predict putative contact sites of the antigen-receptor complex. We also discriminate immunogenic and non-immunogenic antigens, antigen-specific and generic receptors, reaching performances that compare favorably to existing sequence-based predictors of antigen immunogenicity and T-cell receptor specificity. More broadly, diffRBM provides a general framework to detect, interpret and leverage selected features in biological data.

## 1 Introduction

T cells play an essential role in the immune response to pathogens and malignancies. Killer T cells are activated following the binding of their surface receptors (T-cell receptors or TCRs) to short portions of pathogen-related proteins (peptide antigens) that are presented by class I major histocompatibility complexes (MHCs) forming the peptide-MHC epitope (pMHC).

Only a fraction of peptides presented by the MHC are immunogenic, meaning that they possess biochemical properties that can promote a T-cell response [Sette et al., 1994]. Accurate prediction of immunogenicity is crucial to the successful identification of microbial antigens and cancer neoantigens (antigens carrying cancer-related mutations) that help develop vaccines and immune-based cancer therapies. A very recent systematic assessment [Buckley et al., 2022] of the available models to identify immunogenic targets from pathogens and cancers reports suboptimal overall performances, with none of the models able to substantially improve beyond pure MHC-presentation prediction when evaluated on immunogenic peptides from a new virus (SARS-CoV-2). The largest-scale validation [Wells et al., 2020] of existing computational pipelines for neoantigen discovery has highlighted a general lack of consensus among their predictions and a rather low average success rate, with only 6% of the predicted neoantigens validated as truly immunogenic.

We also know that a given pMHC epitope elicits the response of only specific small subsets of the human T-cell repertoire. Predicting the molecular composition of TCRs that have the potential to be reactive to a given epitope is a difficult computational problem that is yet not fully solved despite
# Learning the differences: a transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity

Barbara Bravi¹,², Andrea Di Gioacchino²,⁺, Jorge Fernandez-de-Cossio-Diaz²,⁺, Aleksandra M. Walczak²,*, Thierry Mora²,*, Simona Cocco²,*, and Rémi Monasson²,\*

¹Department of Mathematics, Imperial College London, London, SW7 2AZ United Kingdom

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France

+Equal Contribution

*Equal Contribution

## Abstract

Antigen immunogenicity and the specificity of binding of T-cell receptors to antigens are key properties underlying effective immune responses. Here we propose diffRBM, an approach based on transfer learning and Restricted Boltzmann Machines, to build sequence-based predictive models of these properties. DiffRBM is designed to learn the distinctive patterns in amino acid composition that, one the one hand, underlie the antigen's probability of triggering a response, and on the other hand the T-cell receptor's ability to bind to a given antigen. We show that the patterns learnt by diffRBM allow us to predict putative contact sites of the antigen-receptor complex. We also discriminate immunogenic and non-immunogenic antigens, antigen-specific and generic receptors, reaching performances that compare favorably to existing sequence-based predictors of antigen immunogenicity and T-cell receptor specificity. More broadly, diffRBM provides a general framework to detect, interpret and leverage selected features in biological data.

## 1 Introduction

T cells play an essential role in the immune response to pathogens and malignancies. Killer T cells are activated following the binding of their surface receptors (T-cell receptors or TCRs) to short portions of pathogen-related proteins (peptide antigens) that are presented by class I major histocompatibility complexes (MHCs) forming the peptide-MHC epitope (pMHC).

Only a fraction of peptides presented by the MHC are immunogenic, meaning that they possess biochemical properties that can promote a T-cell response [Sette et al., 1994]. Accurate prediction of immunogenicity is crucial to the successful identification of microbial antigens and cancer neoantigens (antigens carrying cancer-related mutations) that help develop vaccines and immune-based cancer therapies. A very recent systematic assessment [Buckley et al., 2022] of the available models to identify immunogenic targets from pathogens and cancers reports suboptimal overall performances, with none of the models able to substantially improve beyond pure MHC-presentation prediction when evaluated on immunogenic peptides from a new virus (SARS-CoV-2). The largest-scale validation [Wells et al., 2020] of existing computational pipelines for neoantigen discovery has highlighted a general lack of consensus among their predictions and a rather low average success rate, with only 6% of the predicted neoantigens validated as truly immunogenic.

We also know that a given pMHC epitope elicits the response of only specific small subsets of the human T-cell repertoire. Predicting the molecular composition of TCRs that have the potential to be reactive to a given epitope is a difficult computational problem that is yet not fully solved despite
numerous recent advances [Gielis et al., 2019; Jokinen et al., 2019; Springer et al., 2020; Montemurro et al., 2021; Weber et al., 2021]. Improvements in predicting immunogenicity and TCR specificity would have direct consequences for medical applications, including the study of an individual's infection history from their T-cell repertoire and personalized adoptive T-cell therapy for cancer treatment.

Both antigen immunogenicity and epitope-specificity of T-cell receptors have a molecular-level component. They result from specific physico-chemical constraints on the sequence composition of antigens and T-cell receptors. Immunogenic antigens and epitope-specific receptors display an enrichment in specific patterns of amino acid composition. For example, several works have shown enrichment in hydrophobic [Chowell et al., 2015; Riley et al., 2019] and aromatic [Schmidt et al., 2021] residues in immunogenic peptides, compared to all presented peptides (which are predominantly non-binders of TCRs). TCRs specifically responding to the same peptide are also characterized by convergent amino acid motifs [Cinelli et al., 2017; Dash et al., 2017; Glanville et al., 2017; Pogorelyy et al., 2019; Huang et al., 2020; Shomuradova et al., 2020; Mayer-Blackwell et al., 2021; Minervina et al., 2021; Goncharov et al., 2022], whose retrieval is the focus of several clustering approaches [Dash et al., 2017; Glanville et al., 2017; Meysman et al., 2019; Pogorelyy and Shugay, 2019; Thakkar and Bailey-Kellogg, 2019; Mayer-Blackwell et al., 2021; Valkiers et al., 2021]. These patterns are observed in addition to the sequence patterns that are more broadly shared across antigens and TCRs. These general patterns reflect baseline constraints ensuring viability and function (ensuring that TCRs are structurally stable and have the basic binding properties allowing them to pass thymic selection, or ensuring that antigens have high binding affinity to the presenting HLA protein). An outstanding question is how to disentangle sequence pattern enrichment underlying immunogenicity and TCR epitope specificity from the stronger statistical signatures stemming from these baseline constraints. This separation could generate insight into the molecular basis of antigen immunogenicity and epitope specificity and could enable predicting them, at least to some extent, from the sequence alone.

To tackle this question, we here introduce a strategy of 'differential learning' within the unsupervised learning scheme of Restricted Boltzmann Machines [Hinton, 2002; Hinton and Salakhutdinov, 2006], which we call *diffRBM* (differential Restricted Boltzmann Machine). *DiffRBM* relies on a transfer learning procedure, where we first learn general background distributions of antigen or TCR sequences, and then capture, from few sequence data, the distinctive features that confer them, respectively, immunogenicity or epitope specificity. We next interpret the sequence patterns learnt by *diffRBM*, and assess its performance compared to existing computational tools.

## 2 A differential learning approach: *diffRBM*

In this paper, we apply *diffRBM* to build models of peptide immunogenicity and T-cell binding to specific peptides. The modeling setting to which *diffRBM* is applicable is more general. It could be applied to any data that have some distinctive features compared to a much larger pool of data endowed with the baseline properties. We will refer to these two different sets of data as 'selected' and 'background' datasets (Figure 1A-B). In our example, they correspond to sets of immunogenic and only presented antigens in the case of the model of immunogenicity, or antigen-specific and bulk-repertoire TCRs in the case of the model of TCR epitope specificity (Figure 1C-D). Scenarios where it is useful to predict distinctive properties of a subset of data that are important for their functional selection against a background are numerous in biology, with the increasing deployment of SELEX [Ellington and Szostak, 1990; Tuerk and Gold, 1990; Sola et al., 2020] and directed evolution protocols [Jäckel et al., 2008; Packer and Liu, 2015; Arnold, 2018].

Our differential approach is akin to the machine learning technique known as transfer learning, whereby a model learnt for one task is transferred to the second task in such a way that the information embedded in the first model facilitates the learning of the second model. Transfer learning has been recently applied in the context of immune receptors, to model antigen-specific T-cell receptors [Isacchini et al., 2021; Lu et al., 2021] and antibodies [Akbar et al., 2022], but not in the context of modeling antigen immunogenicity. The benefit of a transfer learning approach is two-fold. First, it allows one to leverage the background data availability to control for overfitting. Using the background RBM, which is
numerous recent advances [Gielis et al., 2019; Jokinen et al., 2019; Springer et al., 2020; Montemurro et al., 2021; Weber et al., 2021]. Improvements in predicting immunogenicity and TCR specificity would have direct consequences for medical applications, including the study of an individual's infection history from their T-cell repertoire and personalized adoptive T-cell therapy for cancer treatment.

Both antigen immunogenicity and epitope-specificity of T-cell receptors have a molecular-level component. They result from specific physico-chemical constraints on the sequence composition of antigens and T-cell receptors. Immunogenic antigens and epitope-specific receptors display an enrichment in specific patterns of amino acid composition. For example, several works have shown enrichment in hydrophobic [Chowell et al., 2015; Riley et al., 2019] and aromatic [Schmidt et al., 2021] residues in immunogenic peptides, compared to all presented peptides (which are predominantly non-binders of TCRs). TCRs specifically responding to the same peptide are also characterized by convergent amino acid motifs [Cinelli et al., 2017; Dash et al., 2017; Glanville et al., 2017; Pogorelyy et al., 2019; Huang et al., 2020; Shomuradova et al., 2020; Mayer-Blackwell et al., 2021; Minervina et al., 2021; Goncharov et al., 2022], whose retrieval is the focus of several clustering approaches [Dash et al., 2017; Glanville et al., 2017; Meysman et al., 2019; Pogorelyy and Shugay, 2019; Thakkar and Bailey-Kellogg, 2019; Mayer-Blackwell et al., 2021; Valkiers et al., 2021]. These patterns are observed in addition to the sequence patterns that are more broadly shared across antigens and TCRs. These general patterns reflect baseline constraints ensuring viability and function (ensuring that TCRs are structurally stable and have the basic binding properties allowing them to pass thymic selection, or ensuring that antigens have high binding affinity to the presenting HLA protein). An outstanding question is how to disentangle sequence pattern enrichment underlying immunogenicity and TCR epitope specificity from the stronger statistical signatures stemming from these baseline constraints. This separation could generate insight into the molecular basis of antigen immunogenicity and epitope specificity and could enable predicting them, at least to some extent, from the sequence alone.

To tackle this question, we here introduce a strategy of 'differential learning' within the unsupervised learning scheme of Restricted Boltzmann Machines [Hinton, 2002; Hinton and Salakhutdinov, 2006], which we call *diffRBM* (differential Restricted Boltzmann Machine). *DiffRBM* relies on a transfer learning procedure, where we first learn general background distributions of antigen or TCR sequences, and then capture, from few sequence data, the distinctive features that confer them, respectively, immunogenicity or epitope specificity. We next interpret the sequence patterns learnt by *diffRBM*, and assess its performance compared to existing computational tools.

## 2 A differential learning approach: *diffRBM*

In this paper, we apply *diffRBM* to build models of peptide immunogenicity and T-cell binding to specific peptides. The modeling setting to which *diffRBM* is applicable is more general. It could be applied to any data that have some distinctive features compared to a much larger pool of data endowed with the baseline properties. We will refer to these two different sets of data as ‘selected’ and ‘background’ datasets (Figure 1A-B). In our example, they correspond to sets of immunogenic and only presented antigens in the case of the model of immunogenicity, or antigen-specific and bulk-repertoire TCRs in the case of the model of TCR epitope specificity (Figure 1C-D). Scenarios where it is useful to predict distinctive properties of a subset of data that are important for their functional selection against a background are numerous in biology, with the increasing deployment of SELEX [Ellington and Szostak, 1990; Tuerk and Gold, 1990; Sola et al., 2020] and directed evolution protocols [Jäckel et al., 2008; Packer and Liu, 2015; Arnold, 2018].

Our differential approach is akin to the machine learning technique known as transfer learning, whereby a model learnt for one task is transferred to the second task in such a way that the information embedded in the first model facilitates the learning of the second model. Transfer learning has been recently applied in the context of immune receptors, to model antigen-specific T-cell receptors [Isacchini et al., 2021; Lu et al., 2021] and antibodies [Akbar et al., 2022], but not in the context of modeling antigen immunogenicity. The benefit of a transfer learning approach is two-fold. First, it allows one to leverage the background data availability to control for overfitting. Using the background RBM, which is
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Cartoon of the differential RBM (diffRBM) learning approach. A: The parameters of background RBM (gray) are learnt from the 'background' sequence dataset. B: The diffRBM units (gold) are learnt from a small subset of 'selected' sequences. C: We consider the application of diffRBM to modeling peptide immunogenicity or T-cell receptor (TCR) antigen specificity, whereby the background dataset consists, respectively, of all antigens presented by a given Human Leukocyte Antigen class I complex (HLA) or of generic TCRs from the bulk repertoire. D: The selected sequences correspond to HLA-specific antigens validated to be immunogenic or to TCRs that are antigen-specific responders. The inferred parameters associated to the diffRBM units allow one to identify putative contact positions in the peptide-HLA-TCR structure (E) and more generally to assign scores that distinguish the selected from the background sequences (F). E is an example of a peptide-HLA-TCR structure for the CMV peptide NLVPMVATV (PDB-ID:3GSN), where the contact points along the peptide and the TCR are highlighted in different colors (image obtained with Mol* [Sehnal et al., 2018]).</figcaption></figure>

pre-trained on the larger background dataset as a starting point for learning the diffRBM units, reduces the amount of parameters and data features to infer from the selected dataset (typically much smaller than the background dataset). A second advantage, particularly relevant to the interpretable modeling of biological data, is the fact that the diffRBM units focus on the differences relative to the background, disentangling in this way the data features that make them a selected subset with distinctive properties. As a result, scores based on the diffRBM units can distinguish selected data from the background data (Figure 1F). In the context of modeling antigen immunogenicity and TCR specificity, the diffRBM units and their parameters can also identify antigen-TCR contact points in the three-dimensional molecular structure (Figure 1E).

# 3 Results

## 3.1 DiffRBM model of antigen immunogenicity

Several approaches have been proposed to predict antigen immunogenicity [Calis et al., 2013, Trolle and Nielsen, 2014, Chowell et al., 2015, Ogishi and Yotsuyanagi, 2019, Gao et al., 2020, Schmidt et al., 2021, Lin et al., 2021, Buckley et al., 2022], some of them specifically developed for computational pipelines of neoantigen discovery [Łuksza et al., 2017, Smith et al., 2019, Riley et al., 2019, Schenck et al., 2019, Schaap-Johansen et al., 2021]. The model [Calis et al., 2013] is currently the most used resource for MHC class I immunogenicity prediction, being both implemented by the IEDB tool (http://tools.iedb.org/immunogenicity/) for immunogenicity predictions and integrated in the T cell-antigen interaction prediction by the NetTepi server [Trolle and Nielsen, 2014]. It relies on the simple assumption of position-independence of the main biophysical properties underlying immunogenicity. Some of these approaches are based on a preliminary choice of peptide positions [Schmidt et al., 2021] or properties (such as hydrophobicity [Chowell et al., 2015]), known or assumed to be important for recognition by TCRs. The diffRBM approach we propose
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Cartoon of the differential RBM (diffRBM) learning approach. A: The parameters of background RBM (gray) are learnt from the 'background' sequence dataset. B: The diffRBM units (gold) are learnt from a small subset of 'selected' sequences. C: We consider the application of diffRBM to modeling peptide immunogenicity or T-cell receptor (TCR) antigen specificity, whereby the background dataset consists, respectively, of all antigens presented by a given Human Leukocyte Antigen class I complex (HLA) or of generic TCRs from the bulk repertoire. D: The selected sequences correspond to HLA-specific antigens validated to be immunogenic or to TCRs that are antigen-specific responders. The inferred parameters associated to the diffRBM units allow one to identify putative contact positions in the peptide-HLA-TCR structure (E) and more generally to assign scores that distinguish the selected from the background sequences (F). E is an example of a peptide-HLA-TCR structure for the CMV peptide NLVPMVATV (PDB-ID:3GSN), where the contact points along the peptide and the TCR are highlighted in different colors (image obtained with Mol* [Sehnal et al., 2018]).</figcaption></figure>

pre-trained on the larger background dataset as a starting point for learning the diffRBM units, reduces the amount of parameters and data features to infer from the selected dataset (typically much smaller than the background dataset). A second advantage, particularly relevant to the interpretable modeling of biological data, is the fact that the diffRBM units focus on the differences relative to the background, disentangling in this way the data features that make them a selected subset with distinctive properties. As a result, scores based on the diffRBM units can distinguish selected data from the background data (Figure 1F). In the context of modeling antigen immunogenicity and TCR specificity, the diffRBM units and their parameters can also identify antigen-TCR contact points in the three-dimensional molecular structure (Figure 1E).

# 3 Results

## 3.1 DiffRBM model of antigen immunogenicity

Several approaches have been proposed to predict antigen immunogenicity [Calis et al., 2013, Trolle and Nielsen, 2014, Chowell et al., 2015, Ogishi and Yotsuyanagi, 2019, Gao et al., 2020, Schmidt et al., 2021, Lin et al., 2021, Buckley et al., 2022], some of them specifically developed for computational pipelines of neoantigen discovery [Łuksza et al., 2017, Smith et al., 2019, Riley et al., 2019, Schenck et al., 2019, Schaap-Johansen et al., 2021]. The model [Calis et al., 2013] is currently the most used resource for MHC class I immunogenicity prediction, being both implemented by the IEDB tool (http://tools.iedb.org/immunogenicity/) for immunogenicity predictions and integrated in the T cell-antigen interaction prediction by the NetTepi server [Trolle and Nielsen, 2014]. It relies on the simple assumption of position-independence of the main biophysical properties underlying immunogenicity. Some of these approaches are based on a preliminary choice of peptide positions [Schmidt et al., 2021] or properties (such as hydrophobicity [Chowell et al., 2015]), known or assumed to be important for recognition by TCRs. The diffRBM approach we propose
is a sequence-based method that does not require an *ad hoc* selection of residue positions or properties and does not rely on the assumption of peptide position independence.

We collected from the Immune Epitope Database (IEDB) [Vita et al., 2019] sets of peptides that elicited a T-cell reaction in T cell assays (referred to as ‘immunogenic’) and sets of peptides that were not T-cell-reactive (‘non-immunogenic’). The peptide-presenting MHCs are specialized proteins coded by highly polymorphic human genes called Human Leukocyte Antigen (HLA) gene. We selected only peptides presented by 3 HLA-I alleles (HLA-A\*02:01, HLA-B\*35:01, HLA-B\*07:02). We chose these HLA-I alleles since they are associated to at least 200 immunogenic peptides in IEDB and at least one TCR-pMHC structure in the Protein Data Bank (Methods). We trained diffRBM for each set of HLA-restricted immunogenic peptides in two steps, by training first a background RBM on samples of antigens presented by one specific HLA via the RBM-based algorithm RBM-MHC [Bravi et al., 2021b], and next by training the diffRBM units on the set of immunogenic antigens of the same HLA type (Figures 2A, S2). Background RBM can predict scores of presentation on the specific HLA under consideration, while the diffRBM units predict scores of immunogenicity.

Differently from some existing approaches to modeling immunogenicity [Calis et al., 2013, Chowell et al., 2015, Schmidt et al., 2021], we train *HLA-specific* models. Preliminary inspection of the datasets revealed that patterns of amino acid enrichment differ between immunogenic peptides presented by different HLAs, apart from some general trend in terms of dominant amino acid properties [Calis et al., 2013, Chowell et al., 2015, Riley et al., 2019, Schmidt et al., 2021, Łuksza et al., 2022]. This is true also when we restrict to the peptide positions known to be relevant for immunogenicity [Schmidt et al., 2021, Rudolph et al., 2006, Calis et al., 2012, Schmidt et al., 2021, Łuksza et al., 2022] (Figure 3A). For instance, the extent to which enrichment in hydrophobicity can discriminate immunogenic from non-immunogenic peptides was observed to vary across HLAs [Buckley et al., 2022], since presentation by certain HLAs alone might be associated to hydrophobic residues, supporting HLA-specific strategies to model immunogenicity. On the practical side, cross-HLA imbalances in the training sets, due to the substantial variation in the availability of immunogenic peptides across HLAs, were found to skew predictions toward the most characterized HLAs, in particular toward HLA-A\*02:01 [Buckley et al., 2022].

Equipped with these single-alleles models of immunogenicity, we can perform two tasks, the prediction of peptide sites in contact with the TCR (Figure 1E) and the classification of immunogenic peptides against non-immunogenic ones (Figure 1F).

## 3.2 Validation of model predictions against TCR-pMHC structures

Since the diffRBM units are inferred to capture the distinctive patterns of peptide immunogenicity, we hypothesized that the associated inferred parameters can be informative about the actual structural properties of the pMHC-TCR complex. To test this hypothesis, we collected a set of resolved crystal structures publicly available in the Protein Data Bank [Berman et al., 2000] describing peptides in complex with the presenting HLA molecule and a cognate TCR. For each of these TCR-pMHC complexes, we estimated the peptide positions in contact with the TCR and the HLA (Methods, Table S1). Figures 2B-C show the frequency of contacts at each peptide position with the HLA (B) and the TCR (C) in peptides presented by HLA-A\*02:01 (the HLA allele to which the large majority of structures is available, 41 over 46 structures, see Methods). These contact frequency distributions highlight that positions 2 and 9 (and 1 and 3 to a lower extent) are the anchor sites for the binding of the peptide to the HLA-I protein (Figure 2B), while central positions (4-8) tend to be in contact to the TCR (Figure 2C), consistently with the analyses of structures reported in [Rudolph et al., 2006, Calis et al., 2012, Schmidt et al., 2021, Milighetti et al., 2021]. Previous measures of TCR functional avidity with mutant peptides indicate that amino acid changes at the peptide central positions impact the most T-cell activation [Hoof et al., 2010, Schmidt et al., 2021, Łuksza et al., 2022], suggesting that these positions are important for TCR response.

**Anchor sites of peptide-HLA binding can be inferred from the amino acid statistics.** We first focus on HLA-peptide binding. Anchor sites for the bond with the HLA constrain the amino acid
is a sequence-based method that does not require an *ad hoc* selection of residue positions or properties and does not rely on the assumption of peptide position independence.

We collected from the Immune Epitope Database (IEDB) [Vita et al., 2019] sets of peptides that elicited a T-cell reaction in T cell assays (referred to as ‘immunogenic’) and sets of peptides that were not T-cell-reactive (‘non-immunogenic’). The peptide-presenting MHCs are specialized proteins coded by highly polymorphic human genes called Human Leukocyte Antigen (HLA) gene. We selected only peptides presented by 3 HLA-I alleles (HLA-A\*02:01, HLA-B\*35:01, HLA-B\*07:02). We chose these HLA-I alleles since they are associated to at least 200 immunogenic peptides in IEDB and at least one TCR-pMHC structure in the Protein Data Bank (Methods). We trained diffRBM for each set of HLA-restricted immunogenic peptides in two steps, by training first a background RBM on samples of antigens presented by one specific HLA via the RBM-based algorithm RBM-MHC [Bravi et al., 2021b], and next by training the diffRBM units on the set of immunogenic antigens of the same HLA type (Figures 2A, S2). Background RBM can predict scores of presentation on the specific HLA under consideration, while the diffRBM units predict scores of immunogenicity.

Differently from some existing approaches to modeling immunogenicity [Calis et al., 2013, Chowell et al., 2015, Schmidt et al., 2021], we train *HLA-specific* models. Preliminary inspection of the datasets revealed that patterns of amino acid enrichment differ between immunogenic peptides presented by different HLAs, apart from some general trend in terms of dominant amino acid properties [Calis et al., 2013, Chowell et al., 2015, Riley et al., 2019, Schmidt et al., 2021, Łuksza et al., 2022]. This is true also when we restrict to the peptide positions known to be relevant for immunogenicity [Schmidt et al., 2021, Rudolph et al., 2006, Calis et al., 2012, Schmidt et al., 2021, Łuksza et al., 2022] (Figure 3A). For instance, the extent to which enrichment in hydrophobicity can discriminate immunogenic from non-immunogenic peptides was observed to vary across HLAs [Buckley et al., 2022], since presentation by certain HLAs alone might be associated to hydrophobic residues, supporting HLA-specific strategies to model immunogenicity. On the practical side, cross-HLA imbalances in the training sets, due to the substantial variation in the availability of immunogenic peptides across HLAs, were found to skew predictions toward the most characterized HLAs, in particular toward HLA-A\*02:01 [Buckley et al., 2022].

Equipped with these single-alleles models of immunogenicity, we can perform two tasks, the prediction of peptide sites in contact with the TCR (Figure 1E) and the classification of immunogenic peptides against non-immunogenic ones (Figure 1F).

## 3.2 Validation of model predictions against TCR-pMHC structures

Since the diffRBM units are inferred to capture the distinctive patterns of peptide immunogenicity, we hypothesized that the associated inferred parameters can be informative about the actual structural properties of the pMHC-TCR complex. To test this hypothesis, we collected a set of resolved crystal structures publicly available in the Protein Data Bank [Berman et al., 2000] describing peptides in complex with the presenting HLA molecule and a cognate TCR. For each of these TCR-pMHC complexes, we estimated the peptide positions in contact with the TCR and the HLA (Methods, Table S1). Figures 2B-C show the frequency of contacts at each peptide position with the HLA (B) and the TCR (C) in peptides presented by HLA-A\*02:01 (the HLA allele to which the large majority of structures is available, 41 over 46 structures, see Methods). These contact frequency distributions highlight that positions 2 and 9 (and 1 and 3 to a lower extent) are the anchor sites for the binding of the peptide to the HLA-I protein (Figure 2B), while central positions (4-8) tend to be in contact to the TCR (Figure 2C), consistently with the analyses of structures reported in [Rudolph et al., 2006, Calis et al., 2012, Schmidt et al., 2021, Milighetti et al., 2021]. Previous measures of TCR functional avidity with mutant peptides indicate that amino acid changes at the peptide central positions impact the most T-cell activation [Hoof et al., 2010, Schmidt et al., 2021, Łuksza et al., 2022], suggesting that these positions are important for TCR response.

**Anchor sites of peptide-HLA binding can be inferred from the amino acid statistics.** We first focus on HLA-peptide binding. Anchor sites for the bond with the HLA constrain the amino acid
usage at those positions across peptides presented by that same HLA, increasing the frequency of the amino acids required for binding (e.g., I, L, V at positions 2 and 9 of HLA-A*02:01 ligands). As a result, the sites occupied by the high-frequency amino acids at those positions are typically anchor sites (Figure 2B). We can translate this simple observation into a recipe to predict peptide contacts with the HLA. For each peptide sequence from the TCR-pMHC-A*02:01 structures, we ranked sequence positions by the frequency of the amino acids estimated from the background dataset of HLA-A*02:01-presented peptides. We took the top ranking positions as putative anchor sites and we verified this prediction against the true contacts by calculating an average Positive Predictive Value (PPV), see the inset of Figure 2B and Methods. This average PPV gives the fraction of ranked positions that correspond to the true contacts averaged over HLA-A*02:01-specific peptides from the available structures as a function of the number of ranked positions (i.e. predicted contacts). The PPV is equal to its maximal value, 1, for the first ranked position and stays close to it for the successive ones (Figure 2B inset, Figure S5).

We can also formulate a model-based prediction of contact sites. We took the background RBM parameters, which are learnt to reproduce the amino acid statistics in HLA-A*02:01-presented peptides, and we used them to define peptide site-specific quantities that we call ‘single-site factors’. These correspond essentially to the background RBM log probabilities of single sites conditional on the rest of the sequence (Methods). Large single-site factors detect the peptide-HLA binding anchor sites, similarly to the amino acid frequency (Figure 2B), hence anchor sites can be predicted as the top ranking sequence positions according to the single-site factors. The average PPV of this prediction is comparable to the amino acid frequency-based prediction (Figure 2B inset).

**The diffRBM single-site factors flag up peptide positions important for immunogenicity.**
Also in sets of immunogenic peptides we expect the statistics at the contact positions with the TCR to reflect the constraint of being in contact, when compared to the statistics of all presented peptides. These constraints are captured by the parameters linked to the diffRBM units, as they are learnt in such a way as to pick up the statistical differences with respect to all presented peptides. In analogy to the prediction of HLA-peptide binding via background RBM, we introduce what we call ‘diffRBM single-site factors’ to predict the single residue potential to establish a contact with the TCR. The diffRBM single-site factors are defined from the parameters of the diffRBM units in such a way as to give approximately log odds-ratios between the full RBM and the background RBM probabilities of a certain residue conditional on the rest of the sequence (equation 17 in Methods). Once evaluated on a peptide, these model-dependent terms provide a measure of the predicted contribution to immunogenicity of the amino acid at each peptide position. Since they are based on the diffRBM units, which are connected simultaneously to all sequence sites, the diffRBM single-site factors deliver effectively single-site predictions while accounting for the sequence context given by all the other sites (hence, for the same position and amino acid, their values can vary between sequences).

Figure 2C shows the average magnitude of the diffRBM single-site factors evaluated on the 41 peptides binding to HLA-A*A02:01, which identifies positions 4-8 as the most relevant for immunogenicity. The pattern of positional contact frequency with the TCR in the same figure supports the structural interpretation of the model’s prediction in terms of binding between peptide and TCR.

The diffRBM units’ prediction recovers the pattern of positions important for immunogenicity without restricting *a priori* the input sequences to a subset of peptide positions already known or assumed to be involved in TCR binding. This is in contrast with existing approaches [Calis et al., 2013; Schmidt et al., 2021] that take advantage of considerations on peptide-HLA binding to choose *a priori* what positions to retain in the formulation of the immunogenicity model. [Calis et al., 2013] excludes the positions that are most commonly the anchor sites in the binding to the HLA (1,2 and 9). [Schmidt et al., 2021] learns a classifier of immunogenicity using only the HLA-specific set of positions (for example 5-8 for HLA-A*02:01) estimated to have a minimal impact on HLA affinity based on the statistical information content of HLA-I binding motifs.

**The diffRBM single-site factors predict peptide contact positions with the TCR.** We ranked sequence positions by the single-site factors' magnitude for each peptide in the TCR-pMHC complexes, using the HLA-specific model corresponding to the peptide's HLA type, and we took the highest ranking
usage at those positions across peptides presented by that same HLA, increasing the frequency of the amino acids required for binding (e.g., I, L, V at positions 2 and 9 of HLA-A\*02:01 ligands). As a result, the sites occupied by the high-frequency amino acids at those positions are typically anchor sites (Figure 2B). We can translate this simple observation into a recipe to predict peptide contacts with the HLA. For each peptide sequence from the TCR-pMHC-A\*02:01 structures, we ranked sequence positions by the frequency of the amino acids estimated from the background dataset of HLA-A\*02:01-presented peptides. We took the top ranking positions as putative anchor sites and we verified this prediction against the true contacts by calculating an average Positive Predictive Value (PPV), see the inset of Figure 2B and Methods. This average PPV gives the fraction of ranked positions that correspond to the true contacts averaged over HLA-A\*02:01-specific peptides from the available structures as a function of the number of ranked positions (i.e. predicted contacts). The PPV is equal to its maximal value, 1, for the first ranked position and stays close to it for the successive ones (Figure 2B inset, Figure S5).

We can also formulate a model-based prediction of contact sites. We took the background RBM parameters, which are learnt to reproduce the amino acid statistics in HLA-A\*02:01-presented peptides, and we used them to define peptide site-specific quantities that we call ‘single-site factors’. These correspond essentially to the background RBM log probabilities of single sites conditional on the rest of the sequence (Methods). Large single-site factors detect the peptide-HLA binding anchor sites, similarly to the amino acid frequency (Figure 2B), hence anchor sites can be predicted as the top ranking sequence positions according to the single-site factors. The average PPV of this prediction is comparable to the amino acid frequency-based prediction (Figure 2B inset).

**The diffRBM single-site factors flag up peptide positions important for immunogenicity.**
Also in sets of immunogenic peptides we expect the statistics at the contact positions with the TCR to reflect the constraint of being in contact, when compared to the statistics of all presented peptides. These constraints are captured by the parameters linked to the diffRBM units, as they are learnt in such a way as to pick up the statistical differences with respect to all presented peptides. In analogy to the prediction of HLA-peptide binding via background RBM, we introduce what we call ‘diffRBM single-site factors’ to predict the single residue potential to establish a contact with the TCR. The diffRBM single-site factors are defined from the parameters of the diffRBM units in such a way as to give approximately log odds-ratios between the full RBM and the background RBM probabilities of a certain residue conditional on the rest of the sequence (equation 17 in Methods). Once evaluated on a peptide, these model-dependent terms provide a measure of the predicted contribution to immunogenicity of the amino acid at each peptide position. Since they are based on the diffRBM units, which are connected simultaneously to all sequence sites, the diffRBM single-site factors deliver effectively single-site predictions while accounting for the sequence context given by all the other sites (hence, for the same position and amino acid, their values can vary between sequences).

Figure 2C shows the average magnitude of the diffRBM single-site factors evaluated on the 41 peptides binding to HLA-A\*A02:01, which identifies positions 4-8 as the most relevant for immunogenicity. The pattern of positional contact frequency with the TCR in the same figure supports the structural interpretation of the model’s prediction in terms of binding between peptide and TCR.

The diffRBM units’ prediction recovers the pattern of positions important for immunogenicity without restricting *a priori* the input sequences to a subset of peptide positions already known or assumed to be involved in TCR binding. This is in contrast with existing approaches [Calis et al., 2013; Schmidt et al., 2021] that take advantage of considerations on peptide-HLA binding to choose *a priori* what positions to retain in the formulation of the immunogenicity model. [Calis et al., 2013] excludes the positions that are most commonly the anchor sites in the binding to the HLA (1,2 and 9). [Schmidt et al., 2021] learns a classifier of immunogenicity using only the HLA-specific set of positions (for example 5-8 for HLA-A\*02:01) estimated to have a minimal impact on HLA affinity based on the statistical information content of HLA-I binding motifs.

**The diffRBM single-site factors predict peptide contact positions with the TCR.** We ranked sequence positions by the single-site factors' magnitude for each peptide in the TCR-pMHC complexes, using the HLA-specific model corresponding to the peptide's HLA type, and we took the highest ranking
positions as predicted contact points. The peptide-averaged PPV for this prediction as a function of the ranked positions (Figure 2D) indicates a model's predictive power substantially higher than the random expectation (Methods).

Differential amino acid usage between immunogenic and non-immunogenic peptides at single positions, each considered independently of each other, has been essential to inspect what biochemical properties contribute to immunogenicity [Calis et al., 2013, Chowell et al., 2015] and has been leveraged to build predictors of immunogenicity [Calis et al., 2013]. We compared the prediction by the diffRBM units to such predictions based on differential amino acid usage (Figure 2C-D). Here the single-site factors are given by the log ratio between the position-specific amino acid frequency in immunogenic peptides of a given HLA type and the one in the set of either all presented peptides or the non-immunogenic peptides with the same HLA type (Methods). Thus these single-site factors describe, under an independent-site assumption, the enrichment in amino acid usage in immunogenic peptides relative to all presented peptides or non-immunogenic ones respectively. The diffRBM units outperform these predictions based on amino acid frequency ratios, as quantified by the average PPV (Figures 2D, S4). For the HLA-A\*02:01 peptides, we also found that the magnitude of the diffRBM single-site factors across positions correlates with the pattern of contact frequency better than predictions based on amino acid frequency ratios (Figure 2C).

Finally, Figure 2D shows that a prediction based only on the frequency of an amino acid in the dataset of immunogenic peptides, similarly to the dataset of all presented peptides, tends to identify as putative contacts the anchor sites of binding to the HLA, leading to completely incorrect predictions concerning the binding sites to the TCR.

## 3.3 DiffRBM encodes molecular-level features of immunogenicity

Next we assessed whether our model predictions of the residues' contribution to immunogenicity, based on single-site factors, are consistent with previous findings [Piepenbrink et al., 2013, Chowell et al., 2015, Riley et al., 2019, Schmidt et al., 2021, Łuksza et al., 2022]. Inspecting the amino acid enrichment with respect to non-immunogenic antigens at the common peptide-TCR contact points (positions 4-8) reveals similarities but also differences compared to other accounts of immunogenicity (Figure 3A). Some of the trends visible in Figure 3A resemble the ones found in previous studies, for instance: the bias, in immunogenic sequences, towards hydrophobic amino acids [Chowell et al., 2015, Riley et al., 2019, Łuksza et al., 2022, Buckley et al., 2022] (especially valine, denoted by symbol V); the abundance of glutamic acid (E) [Calis et al., 2013] and, to a more moderate extent, of tryptophan (W) [Schmidt et al., 2021, Łuksza et al., 2022]; the depletion of the small polar amino acid serine (S) and the positively charged amino acid lysine (K), consistently with the observations respectively in [Chowell et al., 2015] and [Schmidt et al., 2021]. We noted also a number of discrepancies, presumably due to the use of datasets of immunogenic and non-immunogenic peptides that were updated with respect to previously available ones and restricted to sequences of only three selected HLA types. The most striking is the under-representation of the aromatic amino acids phenylalanine (F) and tyrosine (Y), particularly severe in HLA-A\*02:01, in contrast with the experiment-based observations in [Schmidt et al., 2021, Piepenbrink et al., 2013]. (Note however that also other analyses performed on IEDB data, like [Calis et al., 2013, Chowell et al., 2015], did not flag up a significant enrichment in Y).

**The diffRBM single-site factors recover the positive contribution to immunogenicity of key residues.** We then verified whether these apparent discrepancies in terms of site-specific amino acid abundances in our IEDB-derived dataset persisted at the level of model predictions. We considered a few combinations of amino acids and positions along the peptide that were suggested to play a crucial role in T-cell reactivity and binding (in the context of HLA-A\*02:01 epitopes) based on structural or functional analyses. For instance, W at position 6 and F at position 7 [Schmidt et al., 2021], as well as Y at position 8 [Piepenbrink et al., 2013], were observed to form a variety of stabilizing interactions with the TCR. Testing functional avidity of TCRs against peptides harbouring single-point mutations, [Schmidt et al., 2021] detected that F and W at position 5 triggered the strongest activation signal, while [Łuksza et al., 2022] found that the substitution of methionine (M) at position 5 systematically abrogated TCR
positions as predicted contact points. The peptide-averaged PPV for this prediction as a function of the ranked positions (Figure 2D) indicates a model's predictive power substantially higher than the random expectation (Methods).

Differential amino acid usage between immunogenic and non-immunogenic peptides at single positions, each considered independently of each other, has been essential to inspect what biochemical properties contribute to immunogenicity [Calis et al., 2013, Chowell et al., 2015] and has been leveraged to build predictors of immunogenicity [Calis et al., 2013]. We compared the prediction by the diffRBM units to such predictions based on differential amino acid usage (Figure 2C-D). Here the single-site factors are given by the log ratio between the position-specific amino acid frequency in immunogenic peptides of a given HLA type and the one in the set of either all presented peptides or the non-immunogenic peptides with the same HLA type (Methods). Thus these single-site factors describe, under an independent-site assumption, the enrichment in amino acid usage in immunogenic peptides relative to all presented peptides or non-immunogenic ones respectively. The diffRBM units outperform these predictions based on amino acid frequency ratios, as quantified by the average PPV (Figures 2D, S4). For the HLA-A\*02:01 peptides, we also found that the magnitude of the diffRBM single-site factors across positions correlates with the pattern of contact frequency better than predictions based on amino acid frequency ratios (Figure 2C).

Finally, Figure 2D shows that a prediction based only on the frequency of an amino acid in the dataset of immunogenic peptides, similarly to the dataset of all presented peptides, tends to identify as putative contacts the anchor sites of binding to the HLA, leading to completely incorrect predictions concerning the binding sites to the TCR.

## 3.3 DiffRBM encodes molecular-level features of immunogenicity

Next we assessed whether our model predictions of the residues' contribution to immunogenicity, based on single-site factors, are consistent with previous findings [Piepenbrink et al., 2013, Chowell et al., 2015, Riley et al., 2019, Schmidt et al., 2021, Łuksza et al., 2022]. Inspecting the amino acid enrichment with respect to non-immunogenic antigens at the common peptide-TCR contact points (positions 4-8) reveals similarities but also differences compared to other accounts of immunogenicity (Figure 3A). Some of the trends visible in Figure 3A resemble the ones found in previous studies, for instance: the bias, in immunogenic sequences, towards hydrophobic amino acids [Chowell et al., 2015, Riley et al., 2019, Łuksza et al., 2022, Buckley et al., 2022] (especially valine, denoted by symbol V); the abundance of glutamic acid (E) [Calis et al., 2013] and, to a more moderate extent, of tryptophan (W) [Schmidt et al., 2021, Łuksza et al., 2022]; the depletion of the small polar amino acid serine (S) and the positively charged amino acid lysine (K), consistently with the observations respectively in [Chowell et al., 2015] and [Schmidt et al., 2021]. We noted also a number of discrepancies, presumably due to the use of datasets of immunogenic and non-immunogenic peptides that were updated with respect to previously available ones and restricted to sequences of only three selected HLA types. The most striking is the under-representation of the aromatic amino acids phenylalanine (F) and tyrosine (Y), particularly severe in HLA-A\*02:01, in contrast with the experiment-based observations in [Schmidt et al., 2021, Piepenbrink et al., 2013]. (Note however that also other analyses performed on IEDB data, like [Calis et al., 2013, Chowell et al., 2015], did not flag up a significant enrichment in Y).

**The diffRBM single-site factors recover the positive contribution to immunogenicity of key residues.** We then verified whether these apparent discrepancies in terms of site-specific amino acid abundances in our IEDB-derived dataset persisted at the level of model predictions. We considered a few combinations of amino acids and positions along the peptide that were suggested to play a crucial role in T-cell reactivity and binding (in the context of HLA-A\*02:01 epitopes) based on structural or functional analyses. For instance, W at position 6 and F at position 7 [Schmidt et al., 2021], as well as Y at position 8 [Piepenbrink et al., 2013], were observed to form a variety of stabilizing interactions with the TCR. Testing functional avidity of TCRs against peptides harbouring single-point mutations, [Schmidt et al., 2021] detected that F and W at position 5 triggered the strongest activation signal, while [Łuksza et al., 2022] found that the substitution of methionine (M) at position 5 systematically abrogated TCR
<figure><img src="image_3.png" /></figure>

B HLA-A*02:01-specific peptides

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

D Peptide-residues in contact with TCR

<figure><img src="image_6.png" /><figcaption>**Figure 2:** DiffRBM model of immunogenicity and structural interpretation of its parameters. A: The diffRBM units are learnt from HLA-specific peptides annotated as immunogenic. B: From the background dataset (HLA-A*02:01-presented antigens), we estimate the expected amino acid frequency at each position; here we plot its average value for the 41 HLA-A*02:01-restricted peptides in the available crystal structures (black line) and we compare it to the magnitude of the single-site factors from background RBM (HLA-A*02:01-specific presentation model), averaged over the 41 peptides. The trend of the amino acid frequency and single-site factors resembles the one of the frequency of contacts with the HLA across the 41 structures for each peptide position. Right inset: Average Positive Predictive Value (PPV) for the prediction of peptide positions in contact with the HLA as a function of the number of ranked positions. The PPV is averaged over the 41 HLA-A*02:01-restricted peptides. The random expectation (blue dotted line) gives the average PPV over a uniformly random prediction of each set of contacts (Methods). C: Same as B, where peptide contact positions are the ones with the TCR and the single-site factors are estimated from either the diffRBM units modeling immunogenicity or amino acid frequency ratios (ratio of position-specific amino acid frequency in immunogenic peptides to non-immunogenic peptides or all presented peptides, see legend in D). ρ denotes the coefficient of correlation between the contact frequency distribution and the magnitude of single-site factors. Contact positions with the HLA (respectively TCR) are defined as the peptide positions whose distance to the residues of the HLA (respectively TCR) in the resolved structure is within 3.5 (respectively 4) Å. D: PPV for the prediction of peptide positions in contact with the TCR as a function of the peptide positions ranked by the models' single-site factors. For each peptide, predictions are made using the corresponding HLA-specific model of immunogenicity. The average is over 46 structures (4 for HLA-B*35:01, 41 for HLA-A*02:01, 1 for HLA-B*07:02). All average PPVs are reweighed by a term accounting for the sequence similarity between peptide entries (Methods, Figures S4A-B, S5A-B).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2:** DiffRBM model of immunogenicity and structural interpretation of its parameters. A: The diffRBM units are learnt from HLA-specific peptides annotated as immunogenic. B: From the background dataset (HLA-A*02:01-presented antigens), we estimate the expected amino acid frequency at each position; here we plot its average value for the 41 HLA-A*02:01-restricted peptides in the available crystal structures (black line) and we compare it to the magnitude of the single-site factors from background RBM (HLA-A*02:01-specific presentation model), averaged over the 41 peptides. The trend of the amino acid frequency and single-site factors resembles the one of the frequency of contacts with the HLA across the 41 structures for each peptide position. Right inset: Average Positive Predictive Value (PPV) for the prediction of peptide positions in contact with the HLA as a function of the number of ranked positions. The PPV is averaged over the 41 HLA-A*02:01-restricted peptides. The random expectation (blue dotted line) gives the average PPV over a uniformly random prediction of each set of contacts (Methods). C: Same as B, where peptide contact positions are the ones with the TCR and the single-site factors are estimated from either the diffRBM units modeling immunogenicity or amino acid frequency ratios (ratio of position-specific amino acid frequency in immunogenic peptides to non-immunogenic peptides or all presented peptides, see legend in D). ρ denotes the coefficient of correlation between the contact frequency distribution and the magnitude of single-site factors. Contact positions with the HLA (respectively TCR) are defined as the peptide positions whose distance to the residues of the HLA (respectively TCR) in the resolved structure is within 3.5 (respectively 4) Å. D: PPV for the prediction of peptide positions in contact with the TCR as a function of the peptide positions ranked by the models' single-site factors. For each peptide, predictions are made using the corresponding HLA-specific model of immunogenicity. The average is over 46 structures (4 for HLA-B*35:01, 41 for HLA-A*02:01, 1 for HLA-B*07:02). All average PPVs are reweighed by a term accounting for the sequence similarity between peptide entries (Methods, Figures S4A-B, S5A-B).</figcaption></figure>
response.

For each of these combinations, we identified the sequences possessing the particular amino acid at the combined position and calculated the sequence-dependent diffRBM single-site factors for that position (Figure 3B). For all amino acid/position pairs, the distribution of these single-site factors is skewed toward positive values (Figures 3B, S6A), meaning that the diffRBM units predict a positive contribution to immunogenicity, and the degree of such contribution varies sequence by sequence in relation to the other amino acids. In contrast, measuring purely the amino acid frequency ratio at that position under the independent-site assumption predicts a contribution to immunogenicity that is significantly smaller (Figure 3B). In other words, the diffRBM's ability to capture the sequence context reconciles its predictions with previous findings on residues that are key to immunogenicity. This result emphasizes that the understanding of the biochemical basis of immunogenicity can benefit from modeling approaches able to assess the role of single residues in a sequence context-dependent manner.

**Model predictions are in agreement with data on TCR reactivity to mutant peptides.**
To further corroborate the connection between the model's predictions and T-cell reactivity assays, we considered the data from [Luksza et al., 2022] on the TCR response to one of the highly immunogenic peptides to which our HLA-A\*02:01-specific diffRBM model of immunogenicity can be applied (NLVPMVATV from the human cytomegalovirus). These data measure TCR reactivity to all possible single-site mutants of NLVPMVATV for 3 TCRs specific to it (Methods, Figure 3C). Some of these single-site mutations do not cause a complete loss of TCR response (we refer to them as 'non-lethal mutations'), meaning that TCR reactivity can be recovered by increasing the peptide concentration (Figure 3C). We can estimate a 'cost' for such mutations in terms of the TCR cross-reactivity between NLVPMVATV and its mutants, measured in [Luksza et al., 2022] as the log ratio between the half maximal effective concentration for TCR activation after the peptide has been mutated ($EC_{50}^{MT}$) and the one before the mutation ($EC_{50}^{WT}$), see Figure 3C. Positive costs reflect the increase in concentration needed to recover TCR reactivity and the corresponding decrease in peptide immunogenicity. Other mutations completely destroy peptide-TCR binding, and TCR reactivity cannot be restored even at highest peptide concentrations considered (we refer to them as 'lethal mutations', see Figure 3C).

Lethal mutations tend to occur at the typical TCR-contact positions along the peptide (Figure 3D), hence the diffRBM units, through the corresponding single-site factors, should account for the drop in TCR response upon these mutations. We found that the diffRBM units predict mainly positive costs for these mutations, which matches qualitatively the observed loss of immunogenicity (Figure 3E).

We observed that the magnitude of non-lethal mutations across all TCR-mutant pairs is concentrated at the peptide positions 2 and 9, which are the anchor sites for binding the HLA in HLA-A\*02:01-specific peptides (Figure 3F). We therefore hypothesized that an increase in peptide concentration is needed first to compensate for decreased HLA presentability, hence a model of presentation alone should be capable of predicting their effect. We found that the prediction by our model of presentation by HLA-A\*02:01 (background RBM) correlates significantly to the experimental mutation costs, and the degree of correlation is consistent across all 3 TCRs, given that these mutations disrupt predominantly peptide presentation on the HLA (Figures 3G, S6B). In contrast, the diffRBM units cannot predict these mutational effects concentrated at the anchor sites for presentation (Figure S6B), since its parameters capture the distinctive molecular composition of immunogenic peptides at the central positions (Figure 2C).

## 3.4 DiffRBM discriminates immunogenic vs non-immunogenic peptides

The diffRBM units learn distinctive sequence patterns of immunogenicity, having the background RBM captured the sequence constraints associated to presentability. Such patterns should contribute to distinguish immunogenic from non-immunogenic peptides. We therefore assigned scores of immunogenicity based on the diffRBM units to held-out test sets of peptides (Figure 4A) and we measured the score's ability to discriminate HLA-specific immunogenic peptides from non-immunogenic ones sharing the same HLA specificity in terms of the area under the receiver operating characteristic curve (AUC), see inset of Figure 4A and Methods.
response.

For each of these combinations, we identified the sequences possessing the particular amino acid at the combined position and calculated the sequence-dependent diffRBM single-site factors for that position (Figure 3B). For all amino acid/position pairs, the distribution of these single-site factors is skewed toward positive values (Figures 3B, S6A), meaning that the diffRBM units predict a positive contribution to immunogenicity, and the degree of such contribution varies sequence by sequence in relation to the other amino acids. In contrast, measuring purely the amino acid frequency ratio at that position under the independent-site assumption predicts a contribution to immunogenicity that is significantly smaller (Figure 3B). In other words, the diffRBM's ability to capture the sequence context reconciles its predictions with previous findings on residues that are key to immunogenicity. This result emphasizes that the understanding of the biochemical basis of immunogenicity can benefit from modeling approaches able to assess the role of single residues in a sequence context-dependent manner.

**Model predictions are in agreement with data on TCR reactivity to mutant peptides.**

To further corroborate the connection between the model's predictions and T-cell reactivity assays, we considered the data from [Luksza et al., 2022] on the TCR response to one of the highly immunogenic peptides to which our HLA-A\*02:01-specific diffRBM model of immunogenicity can be applied (NLVPMVATV from the human cytomegalovirus). These data measure TCR reactivity to all possible single-site mutants of NLVPMVATV for 3 TCRs specific to it (Methods, Figure 3C). Some of these single-site mutations do not cause a complete loss of TCR response (we refer to them as 'non-lethal mutations'), meaning that TCR reactivity can be recovered by increasing the peptide concentration (Figure 3C). We can estimate a 'cost' for such mutations in terms of the TCR cross-reactivity between NLVPMVATV and its mutants, measured in [Luksza et al., 2022] as the log ratio between the half maximal effective concentration for TCR activation after the peptide has been mutated ($EC_{50}^{MT}$) and the one before the mutation ($EC_{50}^{WT}$), see Figure 3C. Positive costs reflect the increase in concentration needed to recover TCR reactivity and the corresponding decrease in peptide immunogenicity. Other mutations completely destroy peptide-TCR binding, and TCR reactivity cannot be restored even at highest peptide concentrations considered (we refer to them as 'lethal mutations', see Figure 3C).

Lethal mutations tend to occur at the typical TCR-contact positions along the peptide (Figure 3D), hence the diffRBM units, through the corresponding single-site factors, should account for the drop in TCR response upon these mutations. We found that the diffRBM units predict mainly positive costs for these mutations, which matches qualitatively the observed loss of immunogenicity (Figure 3E).

We observed that the magnitude of non-lethal mutations across all TCR-mutant pairs is concentrated at the peptide positions 2 and 9, which are the anchor sites for binding the HLA in HLA-A\*02:01-specific peptides (Figure 3F). We therefore hypothesized that an increase in peptide concentration is needed first to compensate for decreased HLA presentability, hence a model of presentation alone should be capable of predicting their effect. We found that the prediction by our model of presentation by HLA-A\*02:01 (background RBM) correlates significantly to the experimental mutation costs, and the degree of correlation is consistent across all 3 TCRs, given that these mutations disrupt predominantly peptide presentation on the HLA (Figures 3G, S6B). In contrast, the diffRBM units cannot predict these mutational effects concentrated at the anchor sites for presentation (Figure S6B), since its parameters capture the distinctive molecular composition of immunogenic peptides at the central positions (Figure 2C).

## 3.4 DiffRBM discriminates immunogenic vs non-immunogenic peptides

The diffRBM units learn distinctive sequence patterns of immunogenicity, having the background RBM captured the sequence constraints associated to presentability. Such patterns should contribute to distinguish immunogenic from non-immunogenic peptides. We therefore assigned scores of immunogenicity based on the diffRBM units to held-out test sets of peptides (Figure 4A) and we measured the score's ability to discriminate HLA-specific immunogenic peptides from non-immunogenic ones sharing the same the HLA specificity in terms of the area under the receiver operating characteristic curve (AUC), see inset of Figure 4A and Methods.
A

Amino acid enrichment in immunogenic peptides (positions 4-8)

<figure><img src="image_7.png" /></figure>

B

Predictions for immunogenicity-related residues

<figure><img src="image_8.png" /></figure>

Data on TCR cross-reactivity to NLVPMVATV mutants

From Łuksza et al. Nature 2022

<figure><img src="image_9.png" /><figcaption>**Figure 3: DiffRBM units encode molecular-level features of immunogenicity.** A: Enrichment in amino acid usage in immunogenic peptides relative to non-immunogenic ones calculated as a log frequency ratio across the central positions (4-8) for each HLA type. The color code used in the *x*-axis labels indicates amino acid properties: red = negatively charged (E,D), blue = positively charged (H, K, R), purple = non-charged polar hydrophilic (N, T, S, Q), golden = aromatic (F, W, Y), black = aliphatic hydrophobic (I, L, M, V), green = cysteine (C), grey = tiny (A, G, P) or other. B: Distribution of the diffRBM single-site factors evaluated across the HLA-A\*02:01-specific immunogenic sequences with W at position 5. For all amino acid/position combinations that were found to be key to immunogenicity, the diffRBM single-site factors assume predominantly positive values (see also Figure S6A). This means that diffRBM predicts a positive contribution to immunogenicity of these key residues, in agreement with observations. Taking into account the sequence context of these residues is essential here: the single-site factors given by the immunogenic *vs* non-immunogenic amino acid frequency ratio, which do not include the sequence context (Methods), predict a much lower contribution to immunogenicity, as indicated by the p-values of their difference with respect to the average of the diffRBM single-site factors distribution. The p-value is estimated from the one-sample T-test and the black dashed line marks the reference p-value=0.01. C: Cartoon illustration of the TCR activation curves obtained in [Łuksza et al., 2022] for the wild-type (WT) peptide NLVPMVATV and its single-site mutants (MT). *EC*<sub>50</sub><sup>WT</sup> and *EC*<sub>50</sub><sup>MT</sup> indicate the half-point peptide concentration of the TCR activation curve for respectively WT and MT (Methods). D: The total count of lethal mutations' costs (214 out of the 513 TCR-mutant combinations) is plotted per mutated peptide position. E: The distribution of predictions of the cost of lethal mutations by diffRBM units is significantly shifted towards positive values (Methods). F: Sum of the costs of non-lethal mutations (299 out of the 514 TCR-mutant combinations) per mutated peptide position. G: Experimental mutation costs *vs* the mutation costs predicted by background RBM for the non-lethal mutations in the context of one TCR (TCR1). The Spearmann correlation coefficient *r* is comparable across the 3 TCRs, with a p-value ≤ 10<sup>-6</sup> for all of them (Figure S6B).</figcaption></figure>
A

Amino acid enrichment in immunogenic peptides (positions 4-8)

<figure><img src="image_5.png" /></figure>

B

Predictions for immunogenicity-related residues

<figure><img src="image_6.png" /></figure>

Data on TCR cross-reactivity to NLVPMVATV mutants

From Łuksza et al. Nature 2022

<figure><img src="image_7.png" /><figcaption>**Figure 3: DiffRBM units encode molecular-level features of immunogenicity.** A: Enrichment in amino acid usage in immunogenic peptides relative to non-immunogenic ones calculated as a log frequency ratio across the central positions (4-8) for each HLA type. The color code used in the *x*-axis labels indicates amino acid properties: red = negatively charged (E,D), blue = positively charged (H, K, R), purple = non-charged polar hydrophilic (N, T, S, Q), golden = aromatic (F, W, Y), black = aliphatic hydrophobic (I, L, M, V), green = cysteine (C), grey = tiny (A, G, P) or other. B: Distribution of the diffRBM single-site factors evaluated across the HLA-A\*02:01-specific immunogenic sequences with W at position 5. For all amino acid/position combinations that were found to be key to immunogenicity, the diffRBM single-site factors assume predominantly positive values (see also Figure S6A). This means that diffRBM predicts a positive contribution to immunogenicity of these key residues, in agreement with observations. Taking into account the sequence context of these residues is essential here: the single-site factors given by the immunogenic *vs* non-immunogenic amino acid frequency ratio, which do not include the sequence context (Methods), predict a much lower contribution to immunogenicity, as indicated by the p-values of their difference with respect to the average of the diffRBM single-site factors distribution. The p-value is estimated from the one-sample T-test and the black dashed line marks the reference p-value=0.01. C: Cartoon illustration of the TCR activation curves obtained in [Łuksza et al., 2022] for the wild-type (WT) peptide NLVPMVATV and its single-site mutants (MT). *EC*<sub>50</sub><sup>WT</sup> and *EC*<sub>50</sub><sup>MT</sup> indicate the half-point peptide concentration of the TCR activation curve for respectively WT and MT (Methods). D: The total count of lethal mutations' costs (214 out of the 513 TCR-mutant combinations) is plotted per mutated peptide position. E: The distribution of predictions of the cost of lethal mutations by diffRBM units is significantly shifted towards positive values (Methods). F: Sum of the costs of non-lethal mutations (299 out of the 514 TCR-mutant combinations) per mutated peptide position. G: Experimental mutation costs *vs* the mutation costs predicted by background RBM for the non-lethal mutations in the context of one TCR (TCR1). The Spearmann correlation coefficient *r* is comparable across the 3 TCRs, with a p-value ≤ 10<sup>-6</sup> for all of them (Figure S6B).</figcaption></figure>
In Figure 4B, we compare the diffRBM units' AUC for the HLA-A\*02:01-specific model to the AUC of the full RBM, that includes the background model in its predictions and hence yields a joint representation of the enriched patterns and of the background constraints. The AUC of discrimination progressively decreases, as expected, from the model part that disentangles enriched patterns from the background constraints (diffRBM units), to the full RBM capturing both, to the model fit to the background constraints only (background RBM). (In fact, there is no reason why background RBM should predict anything at all unless presentability and immunogenicity are correlated). This trend indicates that learning the background constraints that are shared by immunogenic and non-immunogenic peptides along with the sequence pattern enrichment distinctive of immunogenicity can act as a confounding factor when we look for the features that characterize and distinguish immunogenic peptides. Also simpler differential models relying on the independent-site assumption, while returning lower AUCs than RBM-based models, exhibit a decrease in AUC between the differential part and the full model (Figure S7). The AUC values and trend remain stable when we compare scores assigned to immunogenic peptides and to peptides from the human proteome (Figure S8A) and when we score immunogenic and non-immunogenic peptides from the same organism having trained the models on the immunogenic peptides from all the other organisms (Figure S9).

Conversely, the diffRBM units are not designed to capture the background constraints (here associated to presentability). As a result, they cannot successfully discriminate presentable antigens from generic peptides that are predominantly non-presentable (like peptides randomly drawn from the human proteome), in contrast to background RBM and the full RBM (Figure S8B). The same trends of discrimination performance are consistently found across the 3 HLA alleles considered (Figures 4C-D, S7, S8A-B).

The prediction by the background model (trained on peptides probed for their binding to a given HLA-I) provides no clear signal that already the binding affinity to the HLA can discriminate immunogenic peptides (average AUC=0.55 for HLA-A\*02:01, 0.50 for HLA-B\*07:02, 0.53 for HLA-B\*35:01, see Figure 4). To further check this prediction, we scored HLA-A\*02:01 peptides by their binding affinity to the HLA through NetMHCpan4.1 [Reynisson et al., 2020] and found a comparable difference in score distributions between immunogenic and non-immunogenic peptides (AUC = 0.54 for HLA-A\*02:01, 0.48 for HLA-B\*07:02, 0.53 for HLA-B\*35:01). Our observation is in line with a recent large-scale mapping of killer T-cell recognition of candidate neopeptides at high HLA affinity [Kristensen et al., 2022], which did not find a significantly different distribution of HLA-binding NetMHCpan scores between immunogenic vs non-immunogenic neopeptides. Other studies, however, have suggested that immunogenic peptides bind more strongly to the HLA compared to non-immunogenic ones, both in the case of viral epitopes [Croft et al., 2019, Buckley et al., 2022] and neo-epitopes [Bjerregaard et al., 2017, Buckley et al., 2022]. More work is needed in the future to clarify the link between binding affinity to the HLA and immunogenicity.

**A deeper classifier reaches optimal performance, but diffRBM stays comparable.** To perform a task of classification of immunogenic *vs* non-immunogenic peptides, it is more effective to leverage the information in both sets of immunogenic and non-immunogenic peptides. In our framework of differential learning, this can be done by training a second set of differential units on the non-immunogenic peptides, which learns differences in their amino acid statistics with respect to the background. The difference of the scores of the diffRBM units trained on immunogenic peptides and the ones trained on non-immunogenic peptides plays as well the role of a score of immunogenicity, expected to be positive when evaluated on immunogenic peptides and negative on non-immunogenic peptides (Methods). We use it to classify immunogenic *vs* non-immunogenic peptides, resorting to the AUC to measure its classification performance (Figure 4B-D). Such performance is, by construction, the same as the difference of the full RBM scores (Methods, Figure S10).

Exploiting the availability of both immunogenic and non-immunogenic peptides, we also trained a deep neural-network-based classifier. The classifier architecture was optimized among several ones of different depth and width (Methods, Figure S11, Table S2) and reaches a performance quantified by a cross-HLA average AUC = 0.71 ± 0.02 (see 'optimal classifier' in Figure 4B-D; here and in the following uncertainties are estimated over several training-test partitions, see Methods). As a term of reference we considered a linear classifier, whose performance (AUC = 0.66 ± 0.02) shows that there is a gain in
In Figure 4B, we compare the diffRBM units' AUC for the HLA-A\*02:01-specific model to the AUC of the full RBM, that includes the background model in its predictions and hence yields a joint representation of the enriched patterns and of the background constraints. The AUC of discrimination progressively decreases, as expected, from the model part that disentangles enriched patterns from the background constraints (diffRBM units), to the full RBM capturing both, to the model fit to the background constraints only (background RBM). (In fact, there is no reason why background RBM should predict anything at all unless presentability and immunogenicity are correlated). This trend indicates that learning the background constraints that are shared by immunogenic and non-immunogenic peptides along with the sequence pattern enrichment distinctive of immunogenicity can act as a confounding factor when we look for the features that characterize and distinguish immunogenic peptides. Also simpler differential models relying on the independent-site assumption, while returning lower AUCs than RBM-based models, exhibit a decrease in AUC between the differential part and the full model (Figure S7). The AUC values and trend remain stable when we compare scores assigned to immunogenic peptides and to peptides from the human proteome (Figure S8A) and when we score immunogenic and non-immunogenic peptides from the same organism having trained the models on the immunogenic peptides from all the other organisms (Figure S9).

Conversely, the diffRBM units are not designed to capture the background constraints (here associated to presentability). As a result, they cannot successfully discriminate presentable antigens from generic peptides that are predominantly non-presentable (like peptides randomly drawn from the human proteome), in contrast to background RBM and the full RBM (Figure S8B). The same trends of discrimination performance are consistently found across the 3 HLA alleles considered (Figures 4C-D, S7, S8A-B).

The prediction by the background model (trained on peptides probed for their binding to a given HLA-I) provides no clear signal that already the binding affinity to the HLA can discriminate immunogenic peptides (average AUC=0.55 for HLA-A\*02:01, 0.50 for HLA-B\*07:02, 0.53 for HLA-B\*35:01, see Figure 4). To further check this prediction, we scored HLA-A\*02:01 peptides by their binding affinity to the HLA through NetMHCpan4.1 [Reynisson et al., 2020] and found a comparable difference in score distributions between immunogenic and non-immunogenic peptides (AUC = 0.54 for HLA-A\*02:01, 0.48 for HLA-B\*07:02, 0.53 for HLA-B\*35:01). Our observation is in line with a recent large-scale mapping of killer T-cell recognition of candidate neopeptides at high HLA affinity [Kristensen et al., 2022], which did not find a significantly different distribution of HLA-binding NetMHCpan scores between immunogenic vs non-immunogenic neopeptides. Other studies, however, have suggested that immunogenic peptides bind more strongly to the HLA compared to non-immunogenic ones, both in the case of viral epitopes [Croft et al., 2019, Buckley et al., 2022] and neo-epitopes [Bjerregaard et al., 2017, Buckley et al., 2022]. More work is needed in the future to clarify the link between binding affinity to the HLA and immunogenicity.

**A deeper classifier reaches optimal performance, but diffRBM stays comparable.** To perform a task of classification of immunogenic *vs* non-immunogenic peptides, it is more effective to leverage the information in both sets of immunogenic and non-immunogenic peptides. In our framework of differential learning, this can be done by training a second set of differential units on the non-immunogenic peptides, which learns differences in their amino acid statistics with respect to the background. The difference of the scores of the diffRBM units trained on immunogenic peptides and the ones trained on non-immunogenic peptides plays as well the role of a score of immunogenicity, expected to be positive when evaluated on immunogenic peptides and negative on non-immunogenic peptides (Methods). We use it to classify immunogenic *vs* non-immunogenic peptides, resorting to the AUC to measure its classification performance (Figure 4B-D). Such performance is, by construction, the same as the difference of the full RBM scores (Methods, Figure S10).

Exploiting the availability of both immunogenic and non-immunogenic peptides, we also trained a deep neural-network-based classifier. The classifier architecture was optimized among several ones of different depth and width (Methods, Figure S11, Table S2) and reaches a performance quantified by a cross-HLA average AUC = 0.71 ± 0.02 (see 'optimal classifier' in Figure 4B-D; here and in the following uncertainties are estimated over several training-test partitions, see Methods). As a term of reference we considered a linear classifier, whose performance (AUC = 0.66 ± 0.02) shows that there is a gain in
resorting to a more flexible architecture in terms of depth and width. The AUC of the optimal classifier sets the maximal predictive performance that can be achieved, with the datasets under consideration, by a supervised method that is trained to discriminate immunogenic and non-immunogenic antigens. Its value (AUC = 0.71 ± 0.02) indicates that the predictability of immunogenicity from peptide sequences is limited, both by data availability and by the fact that sequence patterns along the peptide are not the only determinant of a positive T-cell response. In the future, more exhaustive assessments of peptide immunogenicity should account for the composition of cognate TCRs, peptide expression levels as well as the regulatory dynamics underlying T-cell response in physiological conditions. The performance of the diffRBM units' scores (AUC = 0.69 ± 0.02) is slightly lower than the one of the optimal classifier but higher than the one of the linear classifier (Figure 4B-D). In this comparison, it is important to remark that the deep classifier is optimized for the classification task, while classification through the difference of diffRBM units' scores is based on unsupervised models trained to learn the distribution of interest (the one of immunogenic or non-immunogenic antigens) and identify its distinguishing features from the baseline distribution of presented antigens.

Finally, we compared the discrimination performance of diffRBM to two established methods for immunogenicity prediction, the IEDB immunogenicity tool and PRIME [Schmidt et al., 2021]. Since these methods are not designed to be re-trained on custom datasets such as ours, the comparison cannot be fully consistent in terms of using the *same* training and test sets, as done for all the methods in Figure 4. The IEDB tool for immunogenicity prediction is based on the model by [Calis et al., 2013] and can be downloaded from http://tools.iedb.org/immunogenicity/. By applying it to our set of immunogenic and non-immunogenic peptides, we obtained AUC=0.54 for the HLA-A*02:01-specific peptides, AUC=0.60 for HLA-B*07:02 and AUC=0.57 for HLA-B*35:01, which are all lower than the diffRBM average AUC values (respectively 0.66, 0.65, 0.67, see Figure 4). By applying PRIME [Schmidt et al., 2021] to our set of immunogenic peptides, we found again poorer performance compared to the diffRBM units (average cross-HLA AUC = 0.51, see Methods), despite not having excluded from our test sets the immunogenic peptides that are either in PRIME's or in the IEDB tool's training sets. Differences in the peptide datasets used may contribute to explain this result. For instance, our dataset consists of epitopes with experimentally validated positive T-cell responses from IEDB, mainly of microbial origin, while PRIME's training set [Schmidt et al., 2021] was constructed in such a way as to contain a high proportion of neoantigens.

As an additional point of comparison, Riley et al. [Riley et al., 2019] propose a neural network trained in a supervised way to classify immunogenic against non-immunogenic peptides using sequence as well as structural features of the peptide/HLA complex. When they check the performance of a peptide sequence-only predictor (using nonamer HLA-A*02:01-restricted peptides), they obtain AUC = 0.61 on the training set and AUC = 0.50 on the test. This further confirms that the performance of our diffRBM approach, despite not being directly trained as a discriminator of immunogenic *vs* non-immunogenic peptides, compares favorably to existing sequence-based tools of prediction.

## 3.5 DiffRBM model for T-cell specific binding

The idea of a differential learning of unsupervised models of TCR statistics, designed to describe the process of thymic selection acting on top of the receptor generation via V(D)J recombination, has first been put forward by [Elhanati et al., 2014]. This effort has culminated in the release of the software SONIA [Sethna et al., 2020] and, most recently, of its deeper version soNNia [Isacchini et al., 2021]. Depending on the TCR dataset used for training, these unsupervised models are applicable also as tools to classify TCRs by type (helper *vs* killer) and by epitope specificity [Isacchini et al., 2021]. The computational prediction of the epitope specificity of a given receptor, in particular, is extremely relevant for therapeutic design, and has been attempted through several algorithms [Cinelli et al., 2017, Gielis et al., 2019, Jokinen et al., 2019, Davidsen et al., 2019, Springer et al., 2020, Luu et al., 2021, Sidhom et al., 2021, Chronister et al., 2021, Weber et al., 2021, Lin et al., 2021, Montemurro et al., 2021, Zhang et al., 2021], spanning a variety of learning methods from random forests [Gielis et al., 2019] to deep neural network-like architectures (convolutional neural networks [Montemurro et al., 2021, Sidhom et al., 2021, Moris et al., 2021, Luu et al., 2021], bimodal attention networks [Weber et al., 2021], long short-term
resorting to a more flexible architecture in terms of depth and width. The AUC of the optimal classifier sets the maximal predictive performance that can be achieved, with the datasets under consideration, by a supervised method that is trained to discriminate immunogenic and non-immunogenic antigens. Its value (AUC = 0.71 ± 0.02) indicates that the predictability of immunogenicity from peptide sequences is limited, both by data availability and by the fact that sequence patterns along the peptide are not the only determinant of a positive T-cell response. In the future, more exhaustive assessments of peptide immunogenicity should account for the composition of cognate TCRs, peptide expression levels as well as the regulatory dynamics underlying T-cell response in physiological conditions. The performance of the diffRBM units' scores (AUC = 0.69 ± 0.02) is slightly lower than the one of the optimal classifier but higher than the one of the linear classifier (Figure 4B-D). In this comparison, it is important to remark that the deep classifier is optimized for the classification task, while classification through the difference of diffRBM units' scores is based on unsupervised models trained to learn the distribution of interest (the one of immunogenic or non-immunogenic antigens) and identify its distinguishing features from the baseline distribution of presented antigens.

Finally, we compared the discrimination performance of diffRBM to two established methods for immunogenicity prediction, the IEDB immunogenicity tool and PRIME [Schmidt et al., 2021]. Since these methods are not designed to be re-trained on custom datasets such as ours, the comparison cannot be fully consistent in terms of using the *same* training and test sets, as done for all the methods in Figure 4. The IEDB tool for immunogenicity prediction is based on the model by [Calis et al., 2013] and can be downloaded from http://tools.iedb.org/immunogenicity/. By applying it to our set of immunogenic and non-immunogenic peptides, we obtained AUC=0.54 for the HLA-A\*02:01-specific peptides, AUC=0.60 for HLA-B\*07:02 and AUC=0.57 for HLA-B\*35:01, which are all lower than the diffRBM average AUC values (respectively 0.66, 0.65, 0.67, see Figure 4). By applying PRIME [Schmidt et al., 2021] to our set of immunogenic peptides, we found again poorer performance compared to the diffRBM units (average cross-HLA AUC = 0.51, see Methods), despite not having excluded from our test sets the immunogenic peptides that are either in PRIME's or in the IEDB tool's training sets. Differences in the peptide datasets used may contribute to explain this result. For instance, our dataset consists of epitopes with experimentally validated positive T-cell responses from IEDB, mainly of microbial origin, while PRIME's training set [Schmidt et al., 2021] was constructed in such a way as to contain a high proportion of neoantigens.

As an additional point of comparison, Riley et al. [Riley et al., 2019] propose a neural network trained in a supervised way to classify immunogenic against non-immunogenic peptides using sequence as well as structural features of the peptide/HLA complex. When they check the performance of a peptide sequence-only predictor (using nonamer HLA-A\*02:01-restricted peptides), they obtain AUC = 0.61 on the training set and AUC = 0.50 on the test. This further confirms that the performance of our diffRBM approach, despite not being directly trained as a discriminator of immunogenic *vs* non-immunogenic peptides, compares favorably to existing sequence-based tools of prediction.

## 3.5 DiffRBM model for T-cell specific binding

The idea of a differential learning of unsupervised models of TCR statistics, designed to describe the process of thymic selection acting on top of the receptor generation via V(D)J recombination, has first been put forward by [Elhanati et al., 2014]. This effort has culminated in the release of the software SONIA [Sethna et al., 2020] and, most recently, of its deeper version soNNia [Isacchini et al., 2021]. Depending on the TCR dataset used for training, these unsupervised models are applicable also as tools to classify TCRs by type (helper *vs* killer) and by epitope specificity [Isacchini et al., 2021]. The computational prediction of the epitope specificity of a given receptor, in particular, is extremely relevant for therapeutic design, and has been attempted through several algorithms [Cinelli et al., 2017, Gielis et al., 2019, Jokinen et al., 2019, Davidsen et al., 2019, Springer et al., 2020, Luu et al., 2021, Sidhom et al., 2021, Chronister et al., 2021, Weber et al., 2021, Lin et al., 2021, Montemurro et al., 2021, Zhang et al., 2021], spanning a variety of learning methods from random forests [Gielis et al., 2019] to deep neural network-like architectures (convolutional neural networks [Montemurro et al., 2021, Sidhom et al., 2021, Moris et al., 2021, Luu et al., 2021], bimodal attention networks [Weber et al., 2021], long short-term
Performance at discriminating immunogenic vs non-immunogenic

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 4: Performance at discriminating immunogenic vs non-immunogenic peptides.** A: For a given HLA (e.g. HLA-A*02:01), we assign diffRBM units' scores to held-out sets of immunogenic and non-immunogenic antigens and we measure the discrimination performance via the Area Under the Curve (AUC), see Methods. B: Performance of diffRBM units, full RBM, background RBM and other methods for the HLA-A*02:01-specific dataset as measured by the AUC. We have divided the methods into unsupervised and supervised ones (color code). Unsupervised methods are meant as the methods trained on immunogenic (or presented) only peptides to learn an unsupervised model. Supervised methods are trained on both immunogenic and non-immunogenic to optimally discriminate them. The approach 'diffRBM units (difference)' is intermediate: it exploits the annotation of peptides as immunogenic and non-immunogenic but it is not optimized for the discrimination task. C,D: Same plot as B for the immunogenicity models built for HLA-B*07:02 (C) and HLA-B*35:01 (D). AUC values shown are the average over 50 partitions into training and test sets and error bars give the corresponding standard deviation (Methods).</figcaption></figure>
Performance at discriminating immunogenic vs non-immunogenic

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 4: Performance at discriminating immunogenic vs non-immunogenic peptides.** A: For a given HLA (e.g. HLA-A*02:01), we assign diffRBM units' scores to held-out sets of immunogenic and non-immunogenic antigens and we measure the discrimination performance via the Area Under the Curve (AUC), see Methods. B: Performance of diffRBM units, full RBM, background RBM and other methods for the HLA-A*02:01-specific dataset as measured by the AUC. We have divided the methods into unsupervised and supervised ones (color code). Unsupervised methods are meant as the methods trained on immunogenic (or presented) only peptides to learn an unsupervised model. Supervised methods are trained on both immunogenic and non-immunogenic to optimally discriminate them. The approach 'diffRBM units (difference)' is intermediate: it exploits the annotation of peptides as immunogenic and non-immunogenic but it is not optimized for the discrimination task. C,D: Same plot as B for the immunogenicity models built for HLA-B*07:02 (C) and HLA-B*35:01 (D). AUC values shown are the average over 50 partitions into training and test sets and error bars give the corresponding standard deviation (Methods).</figcaption></figure>
memory networks and autoencoders [Springer et al., 2020, Lu et al., 2021]). These models are typically calibrated towards achieving high predictive power as classifiers of TCR specificity, while less attention is paid to the interpretability of their predictions apart from a few exceptions [Papadopoulou et al., 2022]. Our diffRBM approach differs from and complements these existing methods because, similarly to SONIA, it focuses on the unsupervised discovery of salient amino acid patterns that can allow for biological insight into the molecular properties determining binding specificity.

To train diffRBM models of epitope specificity, we first collected from the VDJdb database [Shugay et al., 2018, Bagaev et al., 2020] datasets of TCRs specific to 4 epitopes (Methods): the M1₅₈ peptide from the influenza virus (with sequence GILGFVFTL), the pp65₄₉₅ from the human cytomegalovirus (CMV, with sequence NLPMVATV), the BMLF1₂₈₀ peptide from the Epstein-Barr virus (EBV, with sequence GLCTLVAML), the peptide from the Spike protein S₂₆₉ from Sars-Cov-2 (with sequence YLQPRTFLL). We limited the search to sequences of the *β* chain of the TCR (TCR*β*), where the sites of binding to the antigen are concentrated in a region called complementarity determining region 3 (CDR3*β*).

The background dataset, in this case, is meant as a typical bulk TCR*β* repertoire in normal conditions (Figure S13). In particular we take, to train background RBM, the repertoire of a hypothetical universal donor that was constructed by [Isacchini et al., 2021] from the TCR*β* repertoires of the large scale study [Emerson et al., 2017] (Methods). The sequence features captured by background RBM concern germline-encoded amino acid usage related to stability and binding constraints (as is the case for the two conserved residues cysteine and phenylalanine delimiting the CDR3*β* region), as well as additional biases in amino acid usage stemming from VDJ recombination and thymic selection.

After having trained the background RBM, we train a set of diffRBM units on each set of epitope-specific CDR3*β* (Figures 5A, S13). By design, these diffRBM units capture antigen-driven convergent sequence features that have been documented in connection to epitope specificity [Dash et al., 2017, Glanville et al., 2017, Meysman et al., 2019, Pogorelyy and Shugay, 2019, Thakkar and Bailey-Kellogg, 2019, Mayer-Blackwell et al., 2021, Valkiers et al., 2021]. As such, similarly to the model of antigen immunogenicity, the diffRBM units can predict contact sites along the CDR3*β* (Figure 1E) and classify specific receptor against generic, predominantly non-specific ones (Figure 1F).

## 3.6 DiffRBM predicts CDR3β contact positions with the peptide

It has been discussed that convergent features in receptors responding to same antigens have structural interpretation in terms of interactions across the peptide-TCR interface [Dash et al., 2017, Glanville et al., 2017] and that the TCR contact regions are dominated by CDR3β residues [Glanville et al., 2017, Ostmeyer et al., 2019, Milighetti et al., 2021] (see also Figure 1E). Starting from the available TCR-pMHC structures already analyzed, we looked at the peptide-TCR contacts lying in the CDR3β region of the TCR, focusing on the structures involving the 4 epitopes under consideration (12 structures in total, see Methods). Figure 5B shows the distribution of contacts along the CDR3β positions, which consists, as already observed in [Glanville et al., 2017, Ostmeyer et al., 2019, Milighetti et al., 2021], of stretches of 3-5 contiguous amino acids in the central part of the CDR3 (6-8 positions from the left and right anchor).

We estimated the diffRBM single-site factors for the CDR3β sequences from the available structures. Their average value across CDR3β sequences concentrates on the CDR3β central positions and well correlates with the contact frequency distribution (Figure 5B). We next ranked the positions by the single-site factor magnitude and we took top ranking positions as positions of predicted contact. Figure 5C shows the PPV averaged over the 12 available structures (Methods). Similarly to the predictions for peptides (Figure 2D), the predictive power of the diffRBM units is superior to the position-specific amino acid frequency ratio between the antigen-specific and the bulk-repertoire set of receptors. The PPV trend stays robust varying the distance cutoff from 4 Å to 5 Å (the value chosen to determine peptide-TCR contacts in other work [Calis et al., 2012, Glanville et al., 2017, Ostmeyer et al., 2019, Milighetti et al., 2021]), see Figure S15D.
memory networks and autoencoders [Springer et al., 2020, Lu et al., 2021]). These models are typically calibrated towards achieving high predictive power as classifiers of TCR specificity, while less attention is paid to the interpretability of their predictions apart from a few exceptions [Papadopoulou et al., 2022]. Our diffRBM approach differs from and complements these existing methods because, similarly to SONIA, it focuses on the unsupervised discovery of salient amino acid patterns that can allow for biological insight into the molecular properties determining binding specificity.

To train diffRBM models of epitope specificity, we first collected from the VDJdb database [Shugay et al., 2018, Bagaev et al., 2020] datasets of TCRs specific to 4 epitopes (Methods): the M1₅₈ peptide from the influenza virus (with sequence GILGFVFTL), the pp65₄₉₅ from the human cytomegalovirus (CMV, with sequence NLVPMVATV), the BMLF1₂₈₀ peptide from the Epstein-Barr virus (EBV, with sequence GLCTLVAML), the peptide from the Spike protein S₂₆₉ from Sars-Cov-2 (with sequence YLQPRTFLL). We limited the search to sequences of the *β* chain of the TCR (TCR*β*), where the sites of binding to the antigen are concentrated in a region called complementarity determining region 3 (CDR3*β*).

The background dataset, in this case, is meant as a typical bulk TCR*β* repertoire in normal conditions (Figure S13). In particular we take, to train background RBM, the repertoire of a hypothetical universal donor that was constructed by [Isacchini et al., 2021] from the TCR*β* repertoires of the large scale study [Emerson et al., 2017] (Methods). The sequence features captured by background RBM concern germline-encoded amino acid usage related to stability and binding constraints (as is the case for the two conserved residues cysteine and phenylalanine delimiting the CDR3*β* region), as well as additional biases in amino acid usage stemming from VDJ recombination and thymic selection.

After having trained the background RBM, we train a set of diffRBM units on each set of epitope-specific CDR3*β* (Figures 5A, S13). By design, these diffRBM units capture antigen-driven convergent sequence features that have been documented in connection to epitope specificity [Dash et al., 2017, Glanville et al., 2017, Meysman et al., 2019, Pogorelyy and Shugay, 2019, Thakkar and Bailey-Kellogg, 2019, Mayer-Blackwell et al., 2021, Valkiers et al., 2021]. As such, similarly to the model of antigen immunogenicity, the diffRBM units can predict contact sites along the CDR3*β* (Figure 1E) and classify specific receptor against generic, predominantly non-specific ones (Figure 1F).

## 3.6 DiffRBM predicts CDR3β contact positions with the peptide

It has been discussed that convergent features in receptors responding to same antigens have structural interpretation in terms of interactions across the peptide-TCR interface [Dash et al., 2017, Glanville et al., 2017] and that the TCR contact regions are dominated by CDR3β residues [Glanville et al., 2017, Ostmeyer et al., 2019, Milighetti et al., 2021] (see also Figure 1E). Starting from the available TCR-pMHC structures already analyzed, we looked at the peptide-TCR contacts lying in the CDR3β region of the TCR, focusing on the structures involving the 4 epitopes under consideration (12 structures in total, see Methods). Figure 5B shows the distribution of contacts along the CDR3β positions, which consists, as already observed in [Glanville et al., 2017, Ostmeyer et al., 2019, Milighetti et al., 2021], of stretches of 3-5 contiguous amino acids in the central part of the CDR3 (6-8 positions from the left and right anchor).

We estimated the diffRBM single-site factors for the CDR3β sequences from the available structures. Their average value across CDR3β sequences concentrates on the CDR3β central positions and well correlates with the contact frequency distribution (Figure 5B). We next ranked the positions by the single-site factor magnitude and we took top ranking positions as positions of predicted contact. Figure 5C shows the PPV averaged over the 12 available structures (Methods). Similarly to the predictions for peptides (Figure 2D), the predictive power of the diffRBM units is superior to the position-specific amino acid frequency ratio between the antigen-specific and the bulk-repertoire set of receptors. The PPV trend stays robust varying the distance cutoff from 4 Å to 5 Å (the value chosen to determine peptide-TCR contacts in other work [Calis et al., 2012, Glanville et al., 2017, Ostmeyer et al., 2019, Milighetti et al., 2021]), see Figure S15D.
<figure><img src="image_13.png" /><figcaption>Figure 5: DiffRBM model of TCR epitope specificity and structural interpretation of its parameters. A: The diffRBM units are learnt from the CDR3β sequences of antigen-specific T-cell receptors. B: Frequency of contact with the peptide of each CDR3β position calculated across 12 structures (2 for YLQPRTFLL, 3 for NLVPMVATV, 1 for GLCTLVAML, 6 for GILGFVFTL). A position along the CDR3β is estimated as the distance either from the left or from the right anchor site. Contacts are estimated from the structures choosing a 4 Å distance cutoff between CDR3β and peptide. We plot, along with the contact frequency distribution, the magnitude of single-site factors based on the diffRBM units or the amino acid frequency ratio (of peptide-specific sequences relative to bulk-repertoire sequences) averaged over the 12 CDR3β. C: Average PPV for the prediction of CDR3β positions in contact with the peptide as a function of the number of ranked positions (the average is calculated over the 12 structures listed in the caption of B). For each CDR3β sequence, predictions are obtained from the single-site factors of the corresponding peptide-specific models of TCR binding (diffRBM or the model based on the amino acid frequency ratio). The PPV is reweighed by a term accounting for the CDR3β sequence similarity (Methods, Figure S15A-B).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 5: DiffRBM model of TCR epitope specificity and structural interpretation of its parameters.** A: The diffRBM units are learnt from the CDR3β sequences of antigen-specific T-cell receptors. B: Frequency of contact with the peptide of each CDR3β position calculated across 12 structures (2 for YLQPRTFLL, 3 for NLVPMVATV, 1 for GLCTLVAML, 6 for GILGFVFTL). A position along the CDR3β is estimated as the distance either from the left or from the right anchor site. Contacts are estimated from the structures choosing a 4 Å distance cutoff between CDR3β and peptide. We plot, along with the contact frequency distribution, the magnitude of single-site factors based on the diffRBM units or the amino acid frequency ratio (of peptide-specific sequences relative to bulk-repertoire sequences) averaged over the 12 CDR3β. C: Average PPV for the prediction of CDR3β positions in contact with the peptide as a function of the number of ranked positions (the average is calculated over the 12 structures listed in the caption of B). For each CDR3β sequence, predictions are obtained from the single-site factors of the corresponding peptide-specific models of TCR binding (diffRBM or the model based on the amino acid frequency ratio). The PPV is reweighed by a term accounting for the CDR3β sequence similarity (Methods, Figure S15A-B).</figcaption></figure>
## 3.7 DiffRBM discriminates specific receptors

We tested the power to discriminate receptors specific to the 4 epitopes under consideration (GILGFVFTL, NLVPMVATV, GLCTLVAML, YLQPRTFLL) from generic sequences drawn from a bulk repertoire (background dataset). These can be seen as a proxy for non-specific receptors, since the large majority of them is not expected to respond to a specific epitope. We measured the performance at discriminating specific from generic receptors by the AUC (Figure 6A) and the results for the 4 epitope-specific models are reported in Figures 6B-E, S16. We observed the same trend described for the models of immuno-genicity, whereby the AUC of discrimination for the diffRBM units is consistently higher than the one for the full RBM. Also in this context, singling out the sequence features associated to epitope specificity, as the differential units do, enhances the model's predictive performance compared to the full RBM, where the information on those features is added to the background constraints. Any discrimination power is lost when using the background RBM, as it should, since it has no information on epitope specificity.

To further check the robustness of our model's performance, we also repeated the training and testing constructing the background model from a different set of bulk-repertoire TCR$\beta$ sequences from healthy donors (the dataset from [Dean et al., 2015]). All the performance metrics are almost unaffected by this change of the background dataset (Figure S17). The average AUC attained by the diffRBM units across the 4 peptide-specific models is 0.83 ± 0.01 with the background dataset from [Emerson et al., 2017] and 0.82 ± 0.01 with the background dataset from [Dean et al., 2015] (uncertainties are estimated over several training-test partitions, see Methods).

**DiffRBM reaches state-of-the-art performance.** For the sake of comparison of diffRBM to other tools, we considered another generative model of antigen-specific repertoires, SONIA [Elhanati et al., 2014, Sethna et al., 2020], and a series of methods that are trained to discriminate target-specific from unspecific receptors (labeled as 'supervised' in Figure 6B-E). First, as baseline supervised method to compare with, we considered a *k*-Nearest Neighbours (*k*-NN) classifier (Methods). In [Weber et al., 2021] a *k*-NN classifier served as a baseline method to determine the predictive power achievable from TCR sequence similarity alone and, despite its simplicity, it was found to outperform some existing methods. Next, we compared diffRBM to state-of-the-art predictors of TCR specificity, namely TCRex [Gielis et al., 2019] and NetTCR-2.0 [Montemurro et al., 2021]. TCRex is a random forest based classifier of TCRs by epitope specificity, while NetTCR belongs to a family of methods [Springer et al., 2020, Weber et al., 2021, Milighetti et al., 2021] that are trained on TCR-antigen pairs to predict a score of binding, but can also be applied to the prediction of whether an unseen TCR is likely to be specific to a given peptide target.

To make a fair comparison of performance, we re-trained all these methods on our training sets and evaluated the AUC of discrimination of specific receptors on our test sets. For TCRex and NetTCR-2.0, we considered their version for the $\beta$-chain only (Methods). The performance of the diffRBM units is fully comparable to TCRex (both with an average AUC across peptides of 0.83 ± 0.01) and only slightly higher than SONIA and *k*-NN (AUC = 0.80 ± 0.01 for both), see Figure 6B-E. We observed a moderately lower performance of NetTCR-2.0 (AUC= 0.73 ± 0.01), which could be partially due to the fact that the CDR3$\beta$ version of NetTCR-2.0 does not take as input the V and J genes. To check the role of including information on V and J, we also re-trained diffRBM using only the CDR3$\beta$ amino acid sequence (Figure S19). The performance of the diffRBM units is reduced (AUC = 0.79 ± 0.01), but still higher than NetTCR-2.0. We note, however, that the strength of methods like NetTCR [Montemurro et al., 2021, Springer et al., 2020, Weber et al., 2021, Milighetti et al., 2021] relies on being trained to be 'pan-specific', that is, predictive across potentially any peptide and not only for specific peptides. An approach like diffRBM (similarly to SONIA) returns peptide-specific models and hence needs to be trained on TCRs known to recognize a certain peptide for any prediction of specificity to that peptide. Furthermore, the performance by diffRBM decreases compared to other methods when we look at a different task, the one of discriminating a set of antigen-specific receptors from receptors with another antigen specificity [Meysman et al., 2022]. This is due to the fact that our approach, in contrast to other methods [Gielis et al., 2019, Montemurro et al., 2021, Weber et al., 2021], is unsupervised, meaning that it is not trained using receptors of different antigen specificity as negatives against which the positives
## 3.7 DiffRBM discriminates specific receptors

We tested the power to discriminate receptors specific to the 4 epitopes under consideration (GILGFVFTL, NLVPMVATV, GLCTLVAML, YLQPRTFLL) from generic sequences drawn from a bulk repertoire (background dataset). These can be seen as a proxy for non-specific receptors, since the large majority of them is not expected to respond to a specific epitope. We measured the performance at discriminating specific from generic receptors by the AUC (Figure 6A) and the results for the 4 epitope-specific models are reported in Figures 6B-E, S16. We observed the same trend described for the models of immuno-genicity, whereby the AUC of discrimination for the diffRBM units is consistently higher than the one for the full RBM. Also in this context, singling out the sequence features associated to epitope specificity, as the differential units do, enhances the model's predictive performance compared to the full RBM, where the information on those features is added to the background constraints. Any discrimination power is lost when using the background RBM, as it should, since it has no information on epitope specificity.

To further check the robustness of our model's performance, we also repeated the training and testing constructing the background model from a different set of bulk-repertoire TCR$\beta$ sequences from healthy donors (the dataset from [Dean et al., 2015]). All the performance metrics are almost unaffected by this change of the background dataset (Figure S17). The average AUC attained by the diffRBM units across the 4 peptide-specific models is 0.83 ± 0.01 with the background dataset from [Emerson et al., 2017] and 0.82 ± 0.01 with the background dataset from [Dean et al., 2015] (uncertainties are estimated over several training-test partitions, see Methods).

**DiffRBM reaches state-of-the-art performance.** For the sake of comparison of diffRBM to other tools, we considered another generative model of antigen-specific repertoires, SONIA [Elhanati et al., 2014, Sethna et al., 2020], and a series of methods that are trained to discriminate target-specific from unspecific receptors (labeled as 'supervised' in Figure 6B-E). First, as baseline supervised method to compare with, we considered a *k*-Nearest Neighbours (*k*-NN) classifier (Methods). In [Weber et al., 2021] a *k*-NN classifier served as a baseline method to determine the predictive power achievable from TCR sequence similarity alone and, despite its simplicity, it was found to outperform some existing methods. Next, we compared diffRBM to state-of-the-art predictors of TCR specificity, namely TCRex [Gielis et al., 2019] and NetTCR-2.0 [Montemurro et al., 2021]. TCRex is a random forest based classifier of TCRs by epitope specificity, while NetTCR belongs to a family of methods [Springer et al., 2020, Weber et al., 2021, Milighetti et al., 2021] that are trained on TCR-antigen pairs to predict a score of binding, but can also be applied to the prediction of whether an unseen TCR is likely to be specific to a given peptide target.

To make a fair comparison of performance, we re-trained all these methods on our training sets and evaluated the AUC of discrimination of specific receptors on our test sets. For TCRex and NetTCR-2.0, we considered their version for the $\beta$-chain only (Methods). The performance of the diffRBM units is fully comparable to TCRex (both with an average AUC across peptides of 0.83 ± 0.01) and only slightly higher than SONIA and *k*-NN (AUC = 0.80 ± 0.01 for both), see Figure 6B-E. We observed a moderately lower performance of NetTCR-2.0 (AUC= 0.73 ± 0.01), which could be partially due to the fact that the CDR3$\beta$ version of NetTCR-2.0 does not take as input the V and J genes. To check the role of including information on V and J, we also re-trained diffRBM using only the CDR3$\beta$ amino acid sequence (Figure S19). The performance of the diffRBM units is reduced (AUC = 0.79 ± 0.01), but still higher than NetTCR-2.0. We note, however, that the strength of methods like NetTCR [Montemurro et al., 2021, Springer et al., 2020, Weber et al., 2021, Milighetti et al., 2021] relies on being trained to be 'pan-specific', that is, predictive across potentially any peptide and not only for specific peptides. An approach like diffRBM (similarly to SONIA) returns peptide-specific models and hence needs to be trained on TCRs known to recognize a certain peptide for any prediction of specificity to that peptide. Furthermore, the performance by diffRBM decreases compared to other methods when we look at a different task, the one of discriminating a set of antigen-specific receptors from receptors with another antigen specificity [Meysman et al., 2022]. This is due to the fact that our approach, in contrast to other methods [Gielis et al., 2019, Montemurro et al., 2021, Weber et al., 2021], is unsupervised, meaning that it is not trained using receptors of different antigen specificity as negatives against which the positives
<figure><img src="image_14.png" /><figcaption>**Figure 6: Performance at discriminating antigen-specific from generic T-cell receptors.** A: For a given epitope model (e.g. the Influenza epitope GILGFVFTL), we assign diffRBM units' scores to held-out sets of antigen-specific CDR3β and generic CDR3β from the bulk repertoire, and we measure the discrimination performance via the Area Under the Curve (AUC), see Methods. B: AUC of the diffRBM units, full RBM, background RBM and other supervised and unsupervised methods trained and tested on CDR3β sequences specific to the Influenza epitope GILGFVFTL. C-E: The performance assessment illustrated in A-B is repeated for the models of specific response to the CMV epitope NLVPMVATV (C) EBV epitope GLCTLVAML (D), and the Sars-Cov-2 epitope YLQPRTFLL (E). AUC values shown are the average over 50 partitions into training and test sets and error bars give the corresponding standard deviation (Methods).</figcaption></figure>

should be discriminated.

Overall, these comparisons support the conclusion that our approach, when used for classifying receptors that are epitope binders vs generic non-binders from a bulk repertoire, reaches the performance of state-of-the-art TCR specificity predictors [Gielis et al., 2019, Montemurro et al., 2021, Isacchini et al., 2021].

# 4 Discussion

We have introduced a framework, based on the probabilistic graphical model of Restricted Boltzmann Machines, aimed at learning what we called ‘differential’ models, which capture differences of selected subsets of sequences with respect to the full set of sequences of similar type (the ‘background’). More specifically, the ‘differential’ part of the model (the diffRBM units) is designed to learn such differences (Figures 1, S1). First, we applied this framework to the problem of modeling antigen immunogenicity, the ability to trigger a positive T-cell response. We showed that the diffRBM units are able to encode several features relevant to immunogenicity. More precisely, we used the parameters associated to the diffRBM units to estimate which positions are more likely to be in contact with the TCR and we showed that they can distinguish immunogenic from non-immunogenic peptides better than a full RBM, their performance being comparable to a supervised model (the ‘optimal classifier’) that we built and optimized for classification.

To work with samples of sufficient size, the diffRBM models of immunogenicity discussed here were trained on IEDB peptide data, which are predominantly pathogenic. The diffRBM approach is in principle applicable also to cancer neoantigen discovery, but the particular models may need to be re-trained on neoantigen-specific training sets, given that substantial variation in method performance and in the immunogenicity-related predicted features between pathogenic and cancer data has been
<figure><img src="image_12.png" /><figcaption>**Figure 6: Performance at discriminating antigen-specific from generic T-cell receptors.** A: For a given epitope model (e.g. the Influenza epitope GILGFVFTL), we assign diffRBM units' scores to held-out sets of antigen-specific CDR3β and generic CDR3β from the bulk repertoire, and we measure the discrimination performance via the Area Under the Curve (AUC), see Methods. B: AUC of the diffRBM units, full RBM, background RBM and other supervised and unsupervised methods trained and tested on CDR3β sequences specific to the Influenza epitope GILGFVFTL. C-E: The performance assessment illustrated in A-B is repeated for the models of specific response to the CMV epitope NLVPMVATV (C) EBV epitope GLCTLVAML (D), and the Sars-Cov-2 epitope YLQPRTFLL (E). AUC values shown are the average over 50 partitions into training and test sets and error bars give the corresponding standard deviation (Methods).</figcaption></figure>

should be discriminated.

Overall, these comparisons support the conclusion that our approach, when used for classifying receptors that are epitope binders vs generic non-binders from a bulk repertoire, reaches the performance of state-of-the-art TCR specificity predictors [Gielis et al., 2019, Montemurro et al., 2021, Isacchini et al., 2021].

# 4 Discussion

We have introduced a framework, based on the probabilistic graphical model of Restricted Boltzmann Machines, aimed at learning what we called ‘differential’ models, which capture differences of selected subsets of sequences with respect to the full set of sequences of similar type (the ‘background’). More specifically, the ‘differential’ part of the model (the diffRBM units) is designed to learn such differences (Figures 1, S1). First, we applied this framework to the problem of modeling antigen immunogenicity, the ability to trigger a positive T-cell response. We showed that the diffRBM units are able to encode several features relevant to immunogenicity. More precisely, we used the parameters associated to the diffRBM units to estimate which positions are more likely to be in contact with the TCR and we showed that they can distinguish immunogenic from non-immunogenic peptides better than a full RBM, their performance being comparable to a supervised model (the ‘optimal classifier’) that we built and optimized for classification.

To work with samples of sufficient size, the diffRBM models of immunogenicity discussed here were trained on IEDB peptide data, which are predominantly pathogenic. The diffRBM approach is in principle applicable also to cancer neoantigen discovery, but the particular models may need to be re-trained on neoantigen-specific training sets, given that substantial variation in method performance and in the immunogenicity-related predicted features between pathogenic and cancer data has been
evidenced [Buckley et al., 2022]. We did an additional analysis as a preliminary check on the diffRBM performance in the cancer setting: we assigned scores to the HLA-A*02:01-presented neoantigens from the TESLA dataset [Wells et al., 2020] (11 immunogenic and 227 non-immunogenic peptides) and we obtained an AUC of immunogenic vs non-immunogenic discrimination differing significantly from the random expectation (0.6 for the diffRBM units' scores and 0.69 for the difference of diffRBM units' scores).

We also learnt differential models from datasets of epitope-specific T-cell receptors to capture the sequence features underlying their ability to bind to a given epitope, and we have successfully tested the models' power to discriminate epitope-specific from generic receptors and to identify CDR3β residues in contact with the antigen. We looked at the task of distinguishing a small number of antigen-specific receptors from the bulk repertoire (forming the background), a task which reflects many typical scenarios of application. Indeed, antigen-specific receptors constitute a minority of the repertoire, with experimental estimates of the fraction of the TCR repertoire reactive to a given epitope ranging between 10⁻⁶ and 10⁻⁴ [Yates, 2014]. It is at this task that diffRBM performs as well as the state-of-the-art methods. The performance by diffRBM decreases when we consider as negatives receptors with a different antigen specificity [Meysman et al., 2022], because the differential units can learn features shared between antigen-specific sequences (e.g., being expressed by CD8+ T cells). More fine-tuned choices of the background (e.g., with only CD8+ T-cell receptors) would be needed to improve performance for this particular task.

The diffRBM model of specific TCR response can be formulated in an equivalent way for the CDR3 region on the other chain of TCRs (α chain). In this case, one uses CDR3α sequences from healthy TCR repertoires and samples of antigen-specific CDR3α as respectively the background and the selected dataset. DiffRBM models for the α chain reach a discrimination performance comparable to the one for the β chain, as shown in [Meysman et al., 2022]. On the same footing, with single-cell TCR sequencing data becoming increasingly available, our approach could also be extended to modeling the pairs of TCR α and β chains, which have been suggested to play a synergistic role in determining antigen specificity [Carter et al., 2019; Montemurro et al., 2021; Milighetti et al., 2021].

Our primary goal was to learn, through diffRBM, generative models that can describe the probability distribution of data with ‘selected’ features compared to ‘background’ data. As such, they can be used to calculate a probabilistic measure of the diversity of the space sampled by the data, the entropy of the inferred probability distribution (Methods), to then estimate the change in diversity due to selection. The entropy values of background and selected sequences obtained through our models of immunogenicity and TCR specificity highlight a reduction of diversity in selected sequences stemming from the enrichment in distinctive sequence patterns (Figure S20). Such a reduction has substantial variability across epitopes in the case of TCR specificity models (Figure S20B), reflecting the existence of more or less well-defined sequence motifs underlying epitope specificity.

In this paper, the application of our approach has been clearly limited by the scarcity of experimentally confirmed immunogenic peptides for the majority of HLA alleles, as well as TCR specific to many other epitopes. Even for the HLA alleles and epitopes we considered, available datasets cover only a small fraction of HLA-specific peptides that have the potential to be immunogenic and TCR that can recognize given peptides, and pool together data obtained through a heterogeneous set of T-cell response assays, a source of additional noise in the data (Methods). Such limitations in the datasets might affect the performance of all the methods we tested. However, we emphasize here the methodological novelty of learning differential models of immunogenicity and antigen specificity. Our study provides robust proof-of-concept results that a differential modeling approach is an effective tool to extract the sequence features underlying the specificity of binding between antigens and T-cell receptors, and to generate insight into the structural basis of this process by identifying putative contact sites from the sequence alone. Furthermore, relying on unsupervised learning has the advantage that negative samples, whose construction is somewhat arbitrary, are not needed. In addition our model is generative, i.e., it could a priori be used to produce new, putative sequences of antigens or TCR with desired properties, testable in experimental setups.

The impact of biochemical properties, such as hydrophobicity, on immunogenicity is mediated by the structural arrangement of the peptide within the HLA binding groove and of TCR chains, hence
evidenced [Buckley et al., 2022]. We did an additional analysis as a preliminary check on the diffRBM performance in the cancer setting: we assigned scores to the HLA-A*02:01-presented neoantigens from the TESLA dataset [Wells et al., 2020] (11 immunogenic and 227 non-immunogenic peptides) and we obtained an AUC of immunogenic vs non-immunogenic discrimination differing significantly from the random expectation (0.6 for the diffRBM units' scores and 0.69 for the difference of diffRBM units' scores).

We also learnt differential models from datasets of epitope-specific T-cell receptors to capture the sequence features underlying their ability to bind to a given epitope, and we have successfully tested the models' power to discriminate epitope-specific from generic receptors and to identify CDR3β residues in contact with the antigen. We looked at the task of distinguishing a small number of antigen-specific receptors from the bulk repertoire (forming the background), a task which reflects many typical scenarios of application. Indeed, antigen-specific receptors constitute a minority of the repertoire, with experimental estimates of the fraction of the TCR repertoire reactive to a given epitope ranging between 10⁻⁶ and 10⁻⁴ [Yates, 2014]. It is at this task that diffRBM performs as well as the state-of-the-art methods. The performance by diffRBM decreases when we consider as negatives receptors with a different antigen specificity [Meysman et al., 2022], because the differential units can learn features shared between antigen-specific sequences (e.g., being expressed by CD8+ T cells). More fine-tuned choices of the background (e.g., with only CD8+ T-cell receptors) would be needed to improve performance for this particular task.

The diffRBM model of specific TCR response can be formulated in an equivalent way for the CDR3 region on the other chain of TCRs (α chain). In this case, one uses CDR3α sequences from healthy TCR repertoires and samples of antigen-specific CDR3α as respectively the background and the selected dataset. DiffRBM models for the α chain reach a discrimination performance comparable to the one for the β chain, as shown in [Meysman et al., 2022]. On the same footing, with single-cell TCR sequencing data becoming increasingly available, our approach could also be extended to modeling the pairs of TCR α and β chains, which have been suggested to play a synergistic role in determining antigen specificity [Carter et al., 2019; Montemurro et al., 2021; Milighetti et al., 2021].

Our primary goal was to learn, through diffRBM, generative models that can describe the probability distribution of data with ‘selected’ features compared to ‘background’ data. As such, they can be used to calculate a probabilistic measure of the diversity of the space sampled by the data, the entropy of the inferred probability distribution (Methods), to then estimate the change in diversity due to selection. The entropy values of background and selected sequences obtained through our models of immunogenicity and TCR specificity highlight a reduction of diversity in selected sequences stemming from the enrichment in distinctive sequence patterns (Figure S20). Such a reduction has substantial variability across epitopes in the case of TCR specificity models (Figure S20B), reflecting the existence of more or less well-defined sequence motifs underlying epitope specificity.

In this paper, the application of our approach has been clearly limited by the scarcity of experimentally confirmed immunogenic peptides for the majority of HLA alleles, as well as TCR specific to many other epitopes. Even for the HLA alleles and epitopes we considered, available datasets cover only a small fraction of HLA-specific peptides that have the potential to be immunogenic and TCR that can recognize given peptides, and pool together data obtained through a heterogeneous set of T-cell response assays, a source of additional noise in the data (Methods). Such limitations in the datasets might affect the performance of all the methods we tested. However, we emphasize here the methodological novelty of learning differential models of immunogenicity and antigen specificity. Our study provides robust proof-of-concept results that a differential modeling approach is an effective tool to extract the sequence features underlying the specificity of binding between antigens and T-cell receptors, and to generate insight into the structural basis of this process by identifying putative contact sites from the sequence alone. Furthermore, relying on unsupervised learning has the advantage that negative samples, whose construction is somewhat arbitrary, are not needed. In addition our model is generative, i.e., it could a priori be used to produce new, putative sequences of antigens or TCR with desired properties, testable in experimental setups.

The impact of biochemical properties, such as hydrophobicity, on immunogenicity is mediated by the structural arrangement of the peptide within the HLA binding groove and of TCR chains, hence
structural features associated to immunogenicity have been shown to be important for predictive power, both at the level of peptide binding to the MHC [Riley et al., 2019] and TCR binding to the pMHC complex [Lin et al., 2021]. Recent work [Lin et al., 2021] offers an approach for predicting immunogenic peptides by incorporating structural information into the estimation of energies of TCR-antigen binding using coarse-grained protein folding models. Its use is essentially limited to TCR-peptide pairs of known crystal structures, narrowing down the scope of its application. The work [Riley et al., 2019] proposes a neural network-based predictor of immunogenicity that uses the peptide sequence as well as information on the peptide/MHC complex (such as van der Waals interactions, Coulomb potentials, hydrogen bond energies, solvent accessible surface areas) estimated by rapid structural modeling with Rosetta [Chaudhury et al., 2010]. There the need for peptide/MHC structure is circumvented by threading generic peptide sequences into pre-defined template structures, a strategy which is successful due to the homogeneity of the backbone conformation across nonamer peptides presented by the same HLA-I, but might become problematic with peptides of different length or for HLA types admitting different binding modes with the peptide [Gfeller et al., 2018]. Given the large availability of peptide sequence data, many of them produced in clinical settings for personalized medicine purposes, it is extremely valuable to develop approaches that can achieve the more broadly applicable and computationally faster task of scoring peptides by their potential to be immunogenic given the peptide sequence alone. Nevertheless, in terms of future work, leveraging structural information on the TCR-pMHC complex and its estimated binding energy along with sequences, as investigated by [Riley et al., 2019, Lin et al., 2021, Milighetti et al., 2021, Karnaukhov et al., 2022], will contribute to improving the modeling strategy presented, especially given the expected increase in the number of available crystallographic structures. These tasks will require *ad hoc* adjustment of the model architecture, to integrate structural and energetic features into its input, and hence goes beyond the scope of the more general and flexible scheme we aimed to proposed here, which relies on a 'background' dataset of large size, and a small subset of selected data with distinctive property that justify their selection. Such flexibility suggests its broad applicability in several contexts, such as SELEX experiments [Ellington and Szostak, 1990, Tuerk and Gold, 1990, Sola et al., 2020], where each round performs a selection of a subset of molecules from the previous, more populated round (the 'background').
structural features associated to immunogenicity have been shown to be important for predictive power, both at the level of peptide binding to the MHC [Riley et al., 2019] and TCR binding to the pMHC complex [Lin et al., 2021]. Recent work [Lin et al., 2021] offers an approach for predicting immunogenic peptides by incorporating structural information into the estimation of energies of TCR-antigen binding using coarse-grained protein folding models. Its use is essentially limited to TCR-peptide pairs of known crystal structures, narrowing down the scope of its application. The work [Riley et al., 2019] proposes a neural network-based predictor of immunogenicity that uses the peptide sequence as well as information on the peptide/MHC complex (such as van der Waals interactions, Coulomb potentials, hydrogen bond energies, solvent accessible surface areas) estimated by rapid structural modeling with Rosetta [Chaudhury et al., 2010]. There the need for peptide/MHC structure is circumvented by threading generic peptide sequences into pre-defined template structures, a strategy which is successful due to the homogeneity of the backbone conformation across nonamer peptides presented by the same HLA-I, but might become problematic with peptides of different length or for HLA types admitting different binding modes with the peptide [Gfeller et al., 2018]. Given the large availability of peptide sequence data, many of them produced in clinical settings for personalized medicine purposes, it is extremely valuable to develop approaches that can achieve the more broadly applicable and computationally faster task of scoring peptides by their potential to be immunogenic given the peptide sequence alone. Nevertheless, in terms of future work, leveraging structural information on the TCR-pMHC complex and its estimated binding energy along with sequences, as investigated by [Riley et al., 2019, Lin et al., 2021, Milighetti et al., 2021, Karnaukhov et al., 2022], will contribute to improving the modeling strategy presented, especially given the expected increase in the number of available crystallographic structures. These tasks will require *ad hoc* adjustment of the model architecture, to integrate structural and energetic features into its input, and hence goes beyond the scope of the more general and flexible scheme we aimed to proposed here, which relies on a 'background' dataset of large size, and a small subset of selected data with distinctive property that justify their selection. Such flexibility suggests its broad applicability in several contexts, such as SELEX experiments [Ellington and Szostak, 1990, Tuerk and Gold, 1990, Sola et al., 2020], where each round performs a selection of a subset of molecules from the previous, more populated round (the 'background').
# 5 Materials and Methods

## 5.1 DiffRBM architecture

The core idea of what we refer to as a 'differential' probabilistic model is to learn a distribution for sequence data $\sigma$ with the parametric form:

$$P(\boldsymbol{\sigma}) = \frac{1}{Z} e^{-\mathcal{H}(\boldsymbol{\sigma})} \qquad \mathcal{H}(\boldsymbol{\sigma}) = \mathcal{H}^{\mathrm{b}}(\boldsymbol{\sigma}) + \mathcal{H}^{\mathrm{d}}(\boldsymbol{\sigma}) \qquad (1)$$

where $Z$ is simply a normalization factor, $Z = \sum_{\sigma} e^{-\mathcal{H}(\sigma)}$. $\mathcal{H}^{\mathrm{b}}$ specifies the background distribution $P^{\mathrm{b}}$ through $P^{\mathrm{b}}(\boldsymbol{\sigma}) = \frac{1}{Z^{\mathrm{b}}} e^{-\mathcal{H}^{\mathrm{b}}(\boldsymbol{\sigma})}$, where $Z^{\mathrm{b}} = \sum_{\sigma} e^{-\mathcal{H}^{\mathrm{b}}(\boldsymbol{\sigma})}$. The background distribution $P^{\mathrm{b}}(\boldsymbol{\sigma})$ is inferred from what we call a 'background' dataset $\mathcal{D}^{\mathrm{b}}$, a large-size dataset that enables $P^{\mathrm{b}}$ to embed the main statistical properties of the data of interest with high accuracy (thanks to its large size). On the other hand, $\mathcal{H}^{\mathrm{d}}$ contains the parameters that are learnt on top of the background distribution from a set of 'selected' data $\mathcal{D}^{\mathrm{s}}$, picking up their statistical differences with respect to the background distribution.

While this formulation allows for a rather generic choice of both $\mathcal{H}^{\mathrm{b}}$ and $\mathcal{H}^{\mathrm{d}}$, here we assume that they are parametrized in terms of Restricted Boltzmann Machines [Hinton, 2002, Hinton and Salakhutdinov, 2006]. Hence $P^{\mathrm{b}}(\boldsymbol{\sigma})$ can be written as:

$$P^{\mathrm{b}}(\boldsymbol{\sigma}) = \int \prod_{\mu=1}^{M^{\mathrm{b}}} dh_{\mu}^{\mathrm{b}} P^{\mathrm{b}}(\boldsymbol{\sigma}, \mathbf{h}^{\mathrm{b}}) \quad P^{\mathrm{b}}(\boldsymbol{\sigma}, \mathbf{h}^{\mathrm{b}}) \sim \exp \left( \sum_{i=1}^{N} g_{i}^{\mathrm{b}}(\sigma_{i}) - \sum_{\mu=1}^{M^{\mathrm{b}}} \mathcal{U}_{\mu}^{\mathrm{b}}(h_{\mu}^{\mathrm{b}}) + \sum_{i, \mu} h_{\mu}^{\mathrm{b}} w_{i\mu}^{\mathrm{b}}(\sigma_{i}) \right) \quad (2)$$

In a Restricted Boltzmann Machine, the probability of the data $P^{\mathrm{b}}(\boldsymbol{\sigma})$ is expressed as the marginal of a joint probability over the data $\boldsymbol{\sigma}$ (the 'observed' sequences of length $N$) and a set of $M^{\mathrm{b}}$ 'hidden' units $\mathbf{h}^{\mathrm{b}}$, playing the role of coordinates of low-dimensional representations of the data. The joint probability is specified by a set of parameters whose value is inferred from the data (here protein sequences), consisting of: a set of single-site biases $g_{i}^{\mathrm{b}}(\sigma_{i})$, capturing the amino acid usage at each sequence position, a potential $\mathcal{U}_{\mu}^{\mathrm{b}}$ for each hidden unit $h_{\mu}^{\mathrm{b}}$ and a set of parameters $w_{i\mu}^{\mathrm{b}}$, called weights, connecting the sites of observed sequences to each hidden unit. Thus for an RBM as background one obtains:

$$\mathcal{H}^{\mathrm{b}}(\boldsymbol{\sigma}) = - \left( \sum_{i=1}^{N} g_{i}^{\mathrm{b}}(\sigma_{i}) + \sum_{\mu=1}^{M^{\mathrm{b}}} \Gamma_{\mu}^{\mathrm{b}}(I_{\mu}^{\mathrm{b}}(\boldsymbol{\sigma})) \right) \qquad (3)$$

where we have set $I_{\mu}^{\mathrm{b}}(\boldsymbol{\sigma}) = \sum_{i} w_{i\mu}^{\mathrm{b}}(\sigma_{i})$ and $\Gamma_{\mu}^{\mathrm{b}}(I^{\mathrm{b}}(\boldsymbol{\sigma})) = \log \int dh^{\mathrm{b}} e^{-\mathcal{U}_{\mu}^{\mathrm{b}}(h^{\mathrm{b}})+h^{\mathrm{b}}I^{\mathrm{b}}(\boldsymbol{\sigma})}$. The parameters $g_{i}^{\mathrm{b}}$ and $w_{i\mu}^{\mathrm{b}}$, as well as the parameters defining the potentials $\mathcal{U}^{\mathrm{b}}$, are learnt from the background dataset $\mathcal{D}^{\mathrm{b}}$ by maximizing the log-likelihood:

$$\frac{1}{|\mathcal{D}^{\mathrm{b}}|} \sum_{\boldsymbol{\sigma} \in \mathcal{D}^{\mathrm{b}}} \log P^{\mathrm{b}}(\boldsymbol{\sigma}) = \langle \log P^{\mathrm{b}}(\boldsymbol{\sigma}) \rangle_{\mathcal{D}^{\mathrm{b}}} \qquad (4)$$

After the background distribution has been learnt, we learn the 'differential' units of the architecture that we have named differential Restricted Boltzmann Machine (diffRBM), corresponding to the probability distribution:

$$P(\boldsymbol{\sigma}, \mathbf{h}^{\mathrm{b}}, \mathbf{h}^{\mathrm{d}}) \sim \exp \left( \sum_{i=1}^{N} (g_{i}^{\mathrm{d}}(\sigma_{i}) + g_{i}^{\mathrm{b}}(\sigma_{i})) - \sum_{\mu'=1}^{M^{\mathrm{d}}} \mathcal{U}_{\mu'}^{\mathrm{d}}(h_{\mu'}^{\mathrm{d}}) - \sum_{\mu=1}^{M^{\mathrm{b}}} \mathcal{U}_{\mu}^{\mathrm{b}}(h_{\mu}^{\mathrm{b}}) + \sum_{i, \mu'} h_{\mu'}^{\mathrm{d}} w_{i\mu'}^{\mathrm{d}}(\sigma_{i}) + \sum_{i, \mu} h_{\mu}^{\mathrm{b}} w_{i\mu}^{\mathrm{b}}(\sigma_{i}) \right) \qquad (5)$$

see also Figure S1. Hence:

$$\mathcal{H}^{\mathrm{d}}(\boldsymbol{\sigma}) = - \left( \sum_{i=1}^{N} g_{i}^{\mathrm{d}}(\sigma_{i}) + \sum_{\mu'=1}^{M^{\mathrm{d}}} \Gamma_{\mu'}^{\mathrm{d}}(I_{\mu'}^{\mathrm{d}}(\boldsymbol{\sigma})) \right) \qquad (6)$$
# 5 Materials and Methods

## 5.1 DiffRBM architecture

The core idea of what we refer to as a 'differential' probabilistic model is to learn a distribution for sequence data $\sigma$ with the parametric form:

$$P(\boldsymbol{\sigma}) = \frac{1}{Z} e^{-\mathcal{H}(\boldsymbol{\sigma})} \qquad \mathcal{H}(\boldsymbol{\sigma}) = \mathcal{H}^{\mathrm{b}}(\boldsymbol{\sigma}) + \mathcal{H}^{\mathrm{d}}(\boldsymbol{\sigma}) \qquad (1)$$

where $Z$ is simply a normalization factor, $Z = \sum_{\sigma} e^{-\mathcal{H}(\sigma)}$. $\mathcal{H}^b$ specifies the background distribution $P^b$ through $P^b(\boldsymbol{\sigma}) = \frac{1}{Z^b} e^{-\mathcal{H}^b(\boldsymbol{\sigma})}$, where $Z^b = \sum_{\sigma} e^{-\mathcal{H}^b(\boldsymbol{\sigma})}$. The background distribution $P^b(\boldsymbol{\sigma})$ is inferred from what we call a 'background' dataset $\mathcal{D}^b$, a large-size dataset that enables $P^b$ to embed the main statistical properties of the data of interest with high accuracy (thanks to its large size). On the other hand, $\mathcal{H}^d$ contains the parameters that are learnt on top of the background distribution from a set of 'selected' data $\mathcal{D}^s$, picking up their statistical differences with respect to the background distribution.

While this formulation allows for a rather generic choice of both $\mathcal{H}^b$ and $\mathcal{H}^d$, here we assume that they are parametrized in terms of Restricted Boltzmann Machines [Hinton, 2002, Hinton and Salakhutdinov, 2006]. Hence $P^b(\boldsymbol{\sigma})$ can be written as:

$$P^{\mathrm{b}}(\boldsymbol{\sigma}) = \int \prod_{\mu=1}^{M^{\mathrm{b}}} dh_{\mu}^{\mathrm{b}} P^{\mathrm{b}}(\boldsymbol{\sigma}, \mathbf{h}^{\mathrm{b}}) \quad P^{\mathrm{b}}(\boldsymbol{\sigma}, \mathbf{h}^{\mathrm{b}}) \sim \exp \left( \sum_{i=1}^{N} g_{i}^{\mathrm{b}}(\sigma_{i}) - \sum_{\mu=1}^{M^{\mathrm{b}}} \mathcal{U}_{\mu}^{\mathrm{b}}(h_{\mu}^{\mathrm{b}}) + \sum_{i, \mu} h_{\mu}^{\mathrm{b}} w_{i\mu}^{\mathrm{b}}(\sigma_{i}) \right) \quad (2)$$

In a Restricted Boltzmann Machine, the probability of the data $P^b(\boldsymbol{\sigma})$ is expressed as the marginal of a joint probability over the data $\boldsymbol{\sigma}$ (the 'observed' sequences of length $N$) and a set of $M^b$ 'hidden' units $\mathbf{h}^b$, playing the role of coordinates of low-dimensional representations of the data. The joint probability is specified by a set of parameters whose value is inferred from the data (here protein sequences), consisting of: a set of single-site biases $g_i^b(\sigma_i)$, capturing the amino acid usage at each sequence position, a potential $\mathcal{U}_\mu^b$ for each hidden unit $h_\mu^b$ and a set of parameters $w_{i\mu}^b$, called weights, connecting the sites of observed sequences to each hidden unit. Thus for an RBM as background one obtains:

$$\mathcal{H}^{\mathrm{b}}(\boldsymbol{\sigma}) = - \left( \sum_{i=1}^{N} g_{i}^{\mathrm{b}}(\sigma_{i}) + \sum_{\mu=1}^{M^{\mathrm{b}}} \Gamma_{\mu}^{\mathrm{b}}(I_{\mu}^{\mathrm{b}}(\boldsymbol{\sigma})) \right) \qquad (3)$$

where we have set $I_\mu^b(\boldsymbol{\sigma}) = \sum_i w_{i\mu}^b(\sigma_i)$ and $\Gamma_\mu^b(I^b(\boldsymbol{\sigma})) = \log \int dh^b e^{-\mathcal{U}_\mu^b(h^b)+h^b I^b(\boldsymbol{\sigma})}$. The parameters $g_i^b$ and $w_{i\mu}^b$, as well as the parameters defining the potentials $\mathcal{U}^b$, are learnt from the background dataset $\mathcal{D}^b$ by maximizing the log-likelihood:

$$\frac{1}{|\mathcal{D}^b|} \sum_{\boldsymbol{\sigma} \in \mathcal{D}^b} \log P^b(\boldsymbol{\sigma}) = \langle \log P^b(\boldsymbol{\sigma}) \rangle_{\mathcal{D}^b} \qquad (4)$$

After the background distribution has been learnt, we learn the 'differential' units of the architecture that we have named differential Restricted Boltzmann Machine (diffRBM), corresponding to the probability distribution:

$$P(\boldsymbol{\sigma}, \mathbf{h}^{\mathrm{b}}, \mathbf{h}^{\mathrm{d}}) \sim \exp \left( \sum_{i=1}^{N} (g_i^{\mathrm{d}}(\sigma_i) + g_i^{\mathrm{b}}(\sigma_i)) - \sum_{\mu'=1}^{M^{\mathrm{d}}} \mathcal{U}_{\mu'}^{\mathrm{d}}(h_{\mu'}^{\mathrm{d}}) - \sum_{\mu=1}^{M^{\mathrm{b}}} \mathcal{U}_{\mu}^{\mathrm{b}}(h_{\mu}^{\mathrm{b}}) + \sum_{i, \mu'} h_{\mu'}^{\mathrm{d}} w_{i\mu'}^{\mathrm{d}}(\sigma_i) + \sum_{i, \mu} h_{\mu}^{\mathrm{b}} w_{i\mu}^{\mathrm{b}}(\sigma_i) \right) \qquad (5)$$

see also Figure S1. Hence:

$$\mathcal{H}^{\mathrm{d}}(\boldsymbol{\sigma}) = - \left( \sum_{i=1}^{N} g_{i}^{\mathrm{d}}(\sigma_{i}) + \sum_{\mu'=1}^{M^{\mathrm{d}}} \Gamma_{\mu'}^{\mathrm{d}}(I_{\mu'}^{\mathrm{d}}(\boldsymbol{\sigma})) \right) \qquad (6)$$
where we have defined $I_{\mu'}^d(\sigma) = \sum_i w_{i\mu'}^d(\sigma_i)$ and $\Gamma_{\mu'}^d(I^d(\sigma)) = \log \int dh^d e^{-U_{\mu'}^d(h^d)+h^d I^d(\sigma)}$. The parameters featuring in $\mathcal{H}^d(\sigma)$ and defining the 'differential' RBM units are learnt from the dataset $\mathcal{D}^s$ by maximizing:

$$ \frac{1}{|\mathcal{D}^s|} \sum_{\sigma \in \mathcal{D}^s} \log P(\sigma) = \langle \log P(\sigma) \rangle_{\mathcal{D}^s} \quad \text{with} \quad P(\sigma) = \int \prod_{\mu=1}^{M^b} dh_{\mu}^b \prod_{\mu'=1}^{M^d} dh_{\mu'}^d P(\sigma, \mathbf{h}^b, \mathbf{h}^d) \quad (7) $$

where $P(\sigma, \mathbf{h}^b, \mathbf{h}^d)$ is given by (5). Let us define a compact notation $\Psi^d$ for the set of parameters associated to the differential units, $\Psi^d = \{g_i^d(\sigma_i), w_{i\mu'}^d(\sigma_i), \xi_{\mu'}^d\}$, where $\xi_{\mu'}^d$ is a shorthand for the parameters specifying the shape of $U_{\mu'}^d$ (in all our implementations, both $U_{\mu'}^b$ and $U_{\mu'}^d$ are set to dReLU potentials, a choice that has been shown to confer high expressivity [Tubiana et al., 2019]). The rules to infer $\Psi^d$ are given by gradient-ascent equations to maximize the likelihood of the post-selection data $\mathcal{D}^s$ under the full RBM model (equation 1):

$$ \frac{\partial}{\partial \Psi^d} \langle \log P(\sigma) \rangle_{\mathcal{D}^s} = \left\langle \frac{\partial}{\partial \Psi^d} \mathcal{H}(\sigma) \right\rangle_m - \left\langle \frac{\partial}{\partial \Psi^d} \mathcal{H}(\sigma) \right\rangle_{\mathcal{D}^s} \quad (8) $$

where $\langle \rangle_m$ indicates the average under the full RBM model, $\langle f(v) \rangle_m = \sum_{\sigma} \frac{1}{Z} f(\sigma) e^{-\mathcal{H}^b(\sigma)-\mathcal{H}^d(\sigma)}$. Since the parameters $\Psi^d$ to learn appear only in $\mathcal{H}^d$, one has $\partial \mathcal{H}(\sigma)/\partial \Psi^d = \partial \mathcal{H}^d(\sigma)/\partial \Psi^d$; the contribution of the background $\mathcal{H}^b$ enters therefore only in the estimation of the average $\langle \rangle_m$, which requires, at each step of the training, to sample configurations with the probability (1).

In summary, the diffRBM architecture (equation 5) is essentially equivalent to an RBM with $M^b+M^d$ hidden units and where the overall observed fields are written as $g_i^b + g_i^d$: the weights $w_{i\mu}^b$ and potentials $U_{\mu}^b$ for the first $M^b$ units (as well as $g_i^b$) are learnt from the background data, then they are kept fixed and the weights $w_{i\mu'}^d$ and potentials $U_{\mu'}^d$ for additional $M^d$ units (as well as $g_i^d$) are learnt from the selected data (Figure S1).

Regularization terms can be added during the training to control the values of the inferred parameters. Here we used a $L_1$-type regularization over the weights, which enforces sparsity to prevent overfitting (see [Tubiana et al., 2019] for details); following the convention in [Tubiana et al., 2019] we denote by $\lambda_2^1$ the coefficient setting the magnitude of such regularization term. In addition, a $L_2$-type regularization on the fields was used for all models to control for large values stemming from undersampling, and following a standard choice [Cocco et al., 2018] its magnitude was set to the inverse of the size of the training dataset (hence to $1/|\mathcal{D}^b|$ for background models and $1/|\mathcal{D}^s|$ for differential models).

The predictions of peptide immunogenicity or epitope specificity rely on the assignment to sequences $\sigma$ of scores. Using (1), the score of the full RBM is given by the log-likelihood:

$$ \mathcal{L}(\sigma) = \log P(\sigma) = -\mathcal{H}^b(\sigma) - \mathcal{H}^d(\sigma) + \text{constant} \quad (9) $$

where $\mathcal{H}^b(\sigma)$ and $\mathcal{H}^d(\sigma)$ are given by (3) and (6) and the constant stands for a sequence-independent term coming from the partition function. Analogously, the background RBM score and the diffRBM units' score are respectively:

$$ \mathcal{L}^b(\sigma) = -\mathcal{H}^b(\sigma) + \text{constant} \quad (10) $$

$$ \mathcal{L}^d(\sigma) = -\mathcal{H}^d(\sigma) + \text{constant} \quad (11) $$

The entropies of background RBM and full RBM were estimated respectively as:

$$ S^{\text{Background}} = -\sum_{\sigma} P^b(\sigma) \log P^b(\sigma) \quad S^{\text{RBM}} = -\sum_{\sigma} P(\sigma) \log P(\sigma) \quad (12) $$

where log is meant as a natural logarithm, hence the entropy values in Figure S20 are expressed in nats.

In terms of software implementation, diffRBM is coded via a series of additional functions to execute a differential learning on top of the RBM Python implementation from [van der Plas et al., 2022] and is available at github.com/cossio/diffRBM. The codes used for its application to modeling antigen immunogenicity and TCR specificity are downloadable from github.com/bravib/diffRBM\_immunogenicity\_TCRspecificity.
where we have defined $I_{\mu'}^d(\sigma) = \sum_i w_{i\mu'}^d(\sigma_i)$ and $\Gamma_{\mu'}^d(I^d(\sigma)) = \log \int dh^d e^{-U_{\mu'}^d(h^d)+h^d I^d(\sigma)}$. The parameters featuring in $\mathcal{H}^d(\sigma)$ and defining the 'differential' RBM units are learnt from the dataset $\mathcal{D}^s$ by maximizing:

$$ \frac{1}{|\mathcal{D}^s|} \sum_{\sigma \in \mathcal{D}^s} \log P(\sigma) = \langle \log P(\sigma) \rangle_{\mathcal{D}^s} \quad \text{with} \quad P(\sigma) = \int \prod_{\mu=1}^{M^b} dh_{\mu}^b \prod_{\mu'=1}^{M^d} dh_{\mu'}^d P(\sigma, \mathbf{h}^b, \mathbf{h}^d) \quad (7) $$

where $P(\sigma, \mathbf{h}^b, \mathbf{h}^d)$ is given by (5). Let us define a compact notation $\Psi^d$ for the set of parameters associated to the differential units, $\Psi^d = \{g_i^d(\sigma_i), w_{i\mu'}^d(\sigma_i), \xi_{\mu'}^d\}$, where $\xi_{\mu'}^d$ is a shorthand for the parameters specifying the shape of $U_{\mu'}^d$ (in all our implementations, both $U_{\mu'}^b$ and $U_{\mu'}^d$ are set to dReLU potentials, a choice that has been shown to confer high expressivity [Tubiana et al., 2019]). The rules to infer $\Psi^d$ are given by gradient-ascent equations to maximize the likelihood of the post-selection data $\mathcal{D}^s$ under the full RBM model (equation 1):

$$ \frac{\partial}{\partial \Psi^d} \langle \log P(\sigma) \rangle_{\mathcal{D}^s} = \left\langle \frac{\partial}{\partial \Psi^d} \mathcal{H}(\sigma) \right\rangle_m - \left\langle \frac{\partial}{\partial \Psi^d} \mathcal{H}(\sigma) \right\rangle_{\mathcal{D}^s} \quad (8) $$

where $\langle \rangle_m$ indicates the average under the full RBM model, $\langle f(v) \rangle_m = \sum_{\sigma} \frac{1}{Z} f(\sigma) e^{-\mathcal{H}^b(\sigma)-\mathcal{H}^d(\sigma)}$. Since the parameters $\Psi^d$ to learn appear only in $\mathcal{H}^d$, one has $\partial \mathcal{H}(\sigma)/\partial \Psi^d = \partial \mathcal{H}^d(\sigma)/\partial \Psi^d$; the contribution of the background $\mathcal{H}^b$ enters therefore only in the estimation of the average $\langle \rangle_m$, which requires, at each step of the training, to sample configurations with the probability (1).

In summary, the diffRBM architecture (equation 5) is essentially equivalent to an RBM with $M^b+M^d$ hidden units and where the overall observed fields are written as $g_i^b + g_i^d$: the weights $w_{i\mu}^b$ and potentials $U_{\mu}^b$ for the first $M^b$ units (as well as $g_i^b$) are learnt from the background data, then they are kept fixed and the weights $w_{i\mu'}^d$ and potentials $U_{\mu'}^d$ for additional $M^d$ units (as well as $g_i^d$) are learnt from the selected data (Figure S1).

Regularization terms can be added during the training to control the values of the inferred parameters. Here we used a $L_1$-type regularization over the weights, which enforces sparsity to prevent overfitting (see [Tubiana et al., 2019] for details); following the convention in [Tubiana et al., 2019] we denote by $\lambda_2^1$ the coefficient setting the magnitude of such regularization term. In addition, a $L_2$-type regularization on the fields was used for all models to control for large values stemming from undersampling, and following a standard choice [Cocco et al., 2018] its magnitude was set to the inverse of the size of the training dataset (hence to $1/|\mathcal{D}^b|$ for background models and $1/|\mathcal{D}^s|$ for differential models).

The predictions of peptide immunogenicity or epitope specificity rely on the assignment to sequences $\sigma$ of scores. Using (1), the score of the full RBM is given by the log-likelihood:

$$ \mathcal{L}(\sigma) = \log P(\sigma) = -\mathcal{H}^b(\sigma) - \mathcal{H}^d(\sigma) + \text{constant} \quad (9) $$

where $\mathcal{H}^b(\sigma)$ and $\mathcal{H}^d(\sigma)$ are given by (3) and (6) and the constant stands for a sequence-independent term coming from the partition function. Analogously, the background RBM score and the diffRBM units' score are respectively:

$$ \mathcal{L}^b(\sigma) = -\mathcal{H}^b(\sigma) + \text{constant} \quad (10) $$

$$ \mathcal{L}^d(\sigma) = -\mathcal{H}^d(\sigma) + \text{constant} \quad (11) $$

The entropies of background RBM and full RBM were estimated respectively as:

$$ S^{\text{Background}} = -\sum_{\sigma} P^b(\sigma) \log P^b(\sigma) \quad S^{\text{RBM}} = -\sum_{\sigma} P(\sigma) \log P(\sigma) \quad (12) $$

where log is meant as a natural logarithm, hence the entropy values in Figure S20 are expressed in nats.

In terms of software implementation, diffRBM is coded via a series of additional functions to execute a differential learning on top of the RBM Python implementation from [van der Plas et al., 2022] and is available at github.com/cossio/diffRBM. The codes used for its application to modeling antigen immunogenicity and TCR specificity are downloadable from github.com/bravib/diffRBM\_immunogenicity\_TCRspecificity.
## 5.2 Data collection

### 5.2.1 Sequence datasets for the immunogenicity model

The differential models to describe immunogenicity were trained on sets of immunogenic peptides collected from the Immune Epitope Database (IEDB) [Vita et al., 2019], where the database entries were filtered through the following steps. Firstly, the curated set of HLA ligands tested in T cell assays was downloaded from IEDB (file *tcell_full_v3.csv* from http://www.iedb.org/database_export_v3.php, accessed in December 2021). We selected from this file linear, human peptides with a given HLA restriction (e.g. HLA-A*02:01), limiting the search to peptides of length 8-11 amino acids like in [Bravi et al., 2021b] and presented by HLA of class I (i.e., targeted epitopes of killer T cells). Following [Calis et al., 2013], we required the peptide (and not the full protein or the pathogen) to be the first immunogen (by setting the field *Antigen Epitope Relation* = ‘Epitope’) and we excluded T-cell response experiments with a restimulation step (by discarding ‘Restimulation in vitro’ from the field *In Vitro Process Type*). Immunogenic peptides were finally identified as the peptides for which positive responses by T cells were reported while negative ones were absent (field *Qualitative Measure* marked as ‘Positive’ or ‘Positive-High’ and never as ‘Negative’). Equivalently, non-immunogenic peptides were identified as the peptides for which negative responses by T cells were reported while positive ones were absent (field *Qualitative Measure* set to ‘Negative’ and never to ‘Positive’, ‘Positive-High’, ‘Positive-Intermediate’ or ‘Positive-Low’). To avoid oversampling, we removed duplicate entries. To check whether we needed an additional redundancy filtering, similarly to [Calis et al., 2013], we applied a reweighting scheme [Morcos et al., 2011] that reweighs each sequence by the inverse of the number of other sequences that have more than 80% of similarity, and we found that the models’ performance (Figure S12) is largely unchanged compared to the one without reweighting (Figure 4), indicating that there is no substantial need for additional sampling bias mitigation strategies. We choose only the HLA-I alleles for which the filtering steps just described allowed us to recover at least 200 immunogenic peptides and for which at least one TCR-pMHC structure was available in the Protein Data Bank (resulting in the choice of HLA-A*02:01, HLA-B*07:02 and HLA-B*35:01). The size of the final datasets of immunogenic peptides is: |D<sup>s</sup>| = 1682 for HLA-A*02:01, |D<sup>s</sup>| = 258 for HLA-B*07:02, |D<sup>s</sup>| = 215 for HLA-B*35:01. Sets of non-immunogenic peptides consist of 2301 sequences (HLA-A*02:01), 807 (HLA-B*07:02), 166 (HLA-B*35:01).

To train the antigen presentation model, playing here the role of the background model, we relied on the sets of 8-11 amino acid long peptides extracted from IEDB by the RBM-MHC algorithm as described in [Bravi et al., 2021b], choosing the option of peptides from HLA binding affinity assays rather than mass spectrometry, to avoid biases in the amino acid statistics of presented antigens that might be due to mass spectrometry techniques. The resulting training dataset sizes for the background models are |D<sup>b</sup>| = 4265 for HLA-A*02:01, |D<sup>b</sup>| = 1006 for HLA-B*07:02, |D<sup>b</sup>| = 1211 for HLA-B*35:01. For consistency with the type of datasets used to train the presentation model by RBM-MHC, scores of presentation from the algorithm NetMHCpan4.1 [Reynisson et al., 2020] are obtained with the option -BA (predictions from the training on binding assay data).

### 5.2.2 Sequence datasets for the T-cell specificity model

Each differential model of specific T-cell response to a given peptide was trained on sets of T-cell receptors experimentally validated to be specific to the peptide, which were collected from the VDJdb database [Shugay et al., 2018, Bagaev et al., 2020] (file *vdjdb.txt* downloaded from https://vdjdb.cdr3.net in July 2021). More precisely, in the database we selected all the human TCRβ chains fully annotated with their V and J segment and labelled to be specific to the given peptide (for example, for the Influenza M1<sub>58</sub>-specific model we set *antigen.epitope* = ‘GILGFVFTL’). We constructed the final training sets from the CDR3β sequence, the V and J segments identity of these entries, removing replicates with the same CDR3β and V/J annotation. We obtained datasets of epitope-specific receptors of size |D<sup>s</sup>| = 3464 (for the Influenza M1<sub>58</sub> model), |D<sup>s</sup>| = 4548 (for the CMV pp65<sub>495</sub> model), |D<sup>s</sup>| = 993 (for the EBV BMLF1<sub>280</sub> model), and |D<sup>s</sup>| = 315 (for the Sars-Cov-2 S<sub>269</sub> model).

For the background model (here meant as a model of the normal bulk TCRβ repertoire), we considered: the dataset assembled by [Isacchini et al., 2021] pooling together unique TCRβ clones from the
## 5.2 Data collection

### 5.2.1 Sequence datasets for the immunogenicity model

The differential models to describe immunogenicity were trained on sets of immunogenic peptides collected from the Immune Epitope Database (IEDB) [Vita et al., 2019], where the database entries were filtered through the following steps. Firstly, the curated set of HLA ligands tested in T cell assays was downloaded from IEDB (file *tcell_full_v3.csv* from http://www.iedb.org/database_export_v3.php, accessed in December 2021). We selected from this file linear, human peptides with a given HLA restriction (e.g. HLA-A*02:01), limiting the search to peptides of length 8-11 amino acids like in [Bravi et al., 2021b] and presented by HLA of class I (i.e., targeted epitopes of killer T cells). Following [Calis et al., 2013], we required the peptide (and not the full protein or the pathogen) to be the first immunogen (by setting the field *Antigen Epitope Relation* = ‘Epitope’) and we excluded T-cell response experiments with a restimulation step (by discarding ‘Restimulation in vitro’ from the field *In Vitro Process Type*). Immunogenic peptides were finally identified as the peptides for which positive responses by T cells were reported while negative ones were absent (field *Qualitative Measure* marked as ‘Positive’ or ‘Positive-High’ and never as ‘Negative’). Equivalently, non-immunogenic peptides were identified as the peptides for which negative responses by T cells were reported while positive ones were absent (field *Qualitative Measure* set to ‘Negative’ and never to ‘Positive’, ‘Positive-High’, ‘Positive-Intermediate’ or ‘Positive-Low’). To avoid oversampling, we removed duplicate entries. To check whether we needed an additional redundancy filtering, similarly to [Calis et al., 2013], we applied a reweighting scheme [Morcos et al., 2011] that reweighs each sequence by the inverse of the number of other sequences that have more than 80% of similarity, and we found that the models’ performance (Figure S12) is largely unchanged compared to the one without reweighting (Figure 4), indicating that there is no substantial need for additional sampling bias mitigation strategies. We choose only the HLA-I alleles for which the filtering steps just described allowed us to recover at least 200 immunogenic peptides and for which at least one TCR-pMHC structure was available in the Protein Data Bank (resulting in the choice of HLA-A*02:01, HLA-B*07:02 and HLA-B*35:01). The size of the final datasets of immunogenic peptides is: |D<sup>s</sup>| = 1682 for HLA-A*02:01, |D<sup>s</sup>| = 258 for HLA-B*07:02, |D<sup>s</sup>| = 215 for HLA-B*35:01. Sets of non-immunogenic peptides consist of 2301 sequences (HLA-A*02:01), 807 (HLA-B*07:02), 166 (HLA-B*35:01).

To train the antigen presentation model, playing here the role of the background model, we relied on the sets of 8-11 amino acid long peptides extracted from IEDB by the RBM-MHC algorithm as described in [Bravi et al., 2021b], choosing the option of peptides from HLA binding affinity assays rather than mass spectrometry, to avoid biases in the amino acid statistics of presented antigens that might be due to mass spectrometry techniques. The resulting training dataset sizes for the background models are |D<sup>b</sup>| = 4265 for HLA-A*02:01, |D<sup>b</sup>| = 1006 for HLA-B*07:02, |D<sup>b</sup>| = 1211 for HLA-B*35:01. For consistency with the type of datasets used to train the presentation model by RBM-MHC, scores of presentation from the algorithm NetMHCpan4.1 [Reynisson et al., 2020] are obtained with the option -BA (predictions from the training on binding assay data).

### 5.2.2 Sequence datasets for the T-cell specificity model

Each differential model of specific T-cell response to a given peptide was trained on sets of T-cell receptors experimentally validated to be specific to the peptide, which were collected from the VDJdb database [Shugay et al., 2018, Bagaev et al., 2020] (file *vdjdb.txt* downloaded from https://vdjdb.cdr3.net in July 2021). More precisely, in the database we selected all the human TCRβ chains fully annotated with their V and J segment and labelled to be specific to the given peptide (for example, for the Influenza M1<sub>58</sub>-specific model we set *antigen.epitope* = ‘GILGFVFTL’). We constructed the final training sets from the CDR3β sequence, the V and J segments identity of these entries, removing replicates with the same CDR3β and V/J annotation. We obtained datasets of epitope-specific receptors of size |D<sup>s</sup>| = 3464 (for the Influenza M1<sub>58</sub> model), |D<sup>s</sup>| = 4548 (for the CMV pp65<sub>495</sub> model), |D<sup>s</sup>| = 993 (for the EBV BMLF1<sub>280</sub> model), and |D<sup>s</sup>| = 315 (for the Sars-Cov-2 S<sub>269</sub> model).

For the background model (here meant as a model of the normal bulk TCRβ repertoire), we considered: the dataset assembled by [Isacchini et al., 2021] pooling together unique TCRβ clones from the
743 donors of the cohort in [Emerson et al., 2017], with a total of ~ 9 × 10⁷ sequences; the dataset of TCRβ sequences from healthy donors from [Dean et al., 2015], already preprocessed by [Isacchini et al., 2021] and composed of ~ 2 × 10⁶ sequences. For training the background RBM on these data, we used smaller datasets of 10⁶ sequences that could be more easily handled by randomly subsampling from the entire datasets from [Emerson et al., 2017] and [Dean et al., 2015].

### 5.2.3 Data preprocessing and formatting

RBM and PWM-based approaches require sequence inputs of fixed length, hence we performed an alignment. Background datasets are aligned to obtain same-length sequences, following the alignment procedures described in [Bravi et al., 2021b] for peptides and [Bravi et al., 2021a] for CDR3β amino acid sequences. The length of the alignment is set to 9 (9 being the typical length of HLA-I ligands) and to 20 in the case of CDR3β. These alignments serve as seeds to learn Hidden Markov Model profiles of length 9 and 20, in such a way that the selected datasets can be aligned against the profile built from the corresponding background dataset (see [Bravi et al., 2021b] for more details). In the models for TCR sequences, as a default, the input combines the aligned CDR3β amino acid sequence to the V segment type and the J segment type, all converted into numerical values varying within an interval of appropriate length (length = 21 for the CDR3β positions, standing for the 20 amino acids + 1 gap; length = 48 for the V type, length = 13 for the J type).

### 5.2.4 Crystallographic structures from PDB

We downloaded, from the Protein Data Bank www.rcsb.org [Berman et al., 2000] as of February 2022, the TCR-pMHC crystallographic structures with 9 amino acid-long peptides where the HLA complex is HLA-B*35:01, HLA-A*02:01, or HLA-B*07:02. We excluded the structures with modified/non-peptidic epitopes and with incomplete TCR chains. As a result, we obtained 5 structures for HLA-B*35:01, 56 for HLA-A*02:01, and 1 for HLA-B*07:02. Some of the 56 HLA-A*02:01 structures describe TCRs in contact with the peptides we considered for the differential models of specific TCR response (3 for the Sars-Cov-2 epitope YLQPRTFLL, 3 for the CMV epitope NLVPMVATV, 1 for the EBV epitope GLCTLVAML, and 8 for the Influenza epitope GILGFVFTL).

For each structure, we estimated the positions along the peptide in contact with the TCR, using a standard cutoff at 4 Å [Rossjohn et al., 2015, Schmidt et al., 2021, Lu et al., 2021] between heavy atoms. The availability of structures is highly skewed toward the limited set of epitopes that have been the focus of several studies, hence our final list of peptides has some redundancy, with same or similar peptides in complex with different TCRs. If more than 1 structure contain the same peptide and have same contact positions, we retain only one of such structures (resulting in 4 structures for HLA-B*35:01, 41 for HLA-A*02:01, and 1 for HLA-B*07:02). On the other hand, if the same ligand is annotated with different contact positions, we keep these as different entries but we re-weight their contribution to the average PPV and the frequency of contact positions (section 5.7). We followed the same steps to estimate CDR3β contacts with the peptide and peptide contacts with the HLA complex and to filter out redundant entries, opting for a slightly more restrictive cutoff distance (3.5 Å) in the case of peptide-HLA contacts. Since distance cutoffs can vary with the van der Waals' radii for single atoms [Sheriff et al., 1987], we also monitored the robustness of our results to changes in the choice of the cutoff (Figures S4D, S5D, S15D). The list of all structures and corresponding estimated contacts is provided in Table S1.

## 5.3 DiffRBM training and model selection

The first step of the diffRBM training consists of training the background model on the background dataset. For the model of immunogenicity, we trained allele-specific presentation models with an RBM architecture by running the RBM-MHC algorithm [Bravi et al., 2021b] with default parameters (10 hidden units, $\lambda_1^2 = 0.001$). The peptide data used by default as training set by RBM-MHC to build allele-specific presentation models are downloaded from IEDB, see [Bravi et al., 2021b]. The RBM-MHC algorithm internally aligns peptide sequences of its training dataset to the reference length of 9 amino acids.
743 donors of the cohort in [Emerson et al., 2017], with a total of ~ 9 × 10⁷ sequences; the dataset of TCRβ sequences from healthy donors from [Dean et al., 2015], already preprocessed by [Isacchini et al., 2021] and composed of ~ 2 × 10⁶ sequences. For training the background RBM on these data, we used smaller datasets of 10⁶ sequences that could be more easily handled by randomly subsampling from the entire datasets from [Emerson et al., 2017] and [Dean et al., 2015].

### 5.2.3 Data preprocessing and formatting

RBM and PWM-based approaches require sequence inputs of fixed length, hence we performed an alignment. Background datasets are aligned to obtain same-length sequences, following the alignment procedures described in [Bravi et al., 2021b] for peptides and [Bravi et al., 2021a] for CDR3β amino acid sequences. The length of the alignment is set to 9 (9 being the typical length of HLA-I ligands) and to 20 in the case of CDR3β. These alignments serve as seeds to learn Hidden Markov Model profiles of length 9 and 20, in such a way that the selected datasets can be aligned against the profile built from the corresponding background dataset (see [Bravi et al., 2021b] for more details). In the models for TCR sequences, as a default, the input combines the aligned CDR3β amino acid sequence to the V segment type and the J segment type, all converted into numerical values varying within an interval of appropriate length (length = 21 for the CDR3β positions, standing for the 20 amino acids + 1 gap; length = 48 for the V type, length = 13 for the J type).

### 5.2.4 Crystallographic structures from PDB

We downloaded, from the Protein Data Bank www.rcsb.org [Berman et al., 2000] as of February 2022, the TCR-pMHC crystallographic structures with 9 amino acid-long peptides where the HLA complex is HLA-B*35:01, HLA-A*02:01, or HLA-B*07:02. We excluded the structures with modified/non-peptidic epitopes and with incomplete TCR chains. As a result, we obtained 5 structures for HLA-B*35:01, 56 for HLA-A*02:01, and 1 for HLA-B*07:02. Some of the 56 HLA-A*02:01 structures describe TCRs in contact with the peptides we considered for the differential models of specific TCR response (3 for the Sars-Cov-2 epitope YLQPRTFLL, 3 for the CMV epitope NLVPMVATV, 1 for the EBV epitope GLCTLVAML, and 8 for the Influenza epitope GILGFVFTL).

For each structure, we estimated the positions along the peptide in contact with the TCR, using a standard cutoff at 4 Å [Rossjohn et al., 2015, Schmidt et al., 2021, Lu et al., 2021] between heavy atoms. The availability of structures is highly skewed toward the limited set of epitopes that have been the focus of several studies, hence our final list of peptides has some redundancy, with same or similar peptides in complex with different TCRs. If more than 1 structure contain the same peptide and have same contact positions, we retain only one of such structures (resulting in 4 structures for HLA-B*35:01, 41 for HLA-A*02:01, and 1 for HLA-B*07:02). On the other hand, if the same ligand is annotated with different contact positions, we keep these as different entries but we re-weight their contribution to the average PPV and the frequency of contact positions (section 5.7). We followed the same steps to estimate CDR3β contacts with the peptide and peptide contacts with the HLA complex and to filter out redundant entries, opting for a slightly more restrictive cutoff distance (3.5 Å) in the case of peptide-HLA contacts. Since distance cutoffs can vary with the van der Waals' radii for single atoms [Sheriff et al., 1987], we also monitored the robustness of our results to changes in the choice of the cutoff (Figures S4D, S5D, S15D). The list of all structures and corresponding estimated contacts is provided in Table S1.

## 5.3 DiffRBM training and model selection

The first step of the diffRBM training consists of training the background model on the background dataset. For the model of immunogenicity, we trained allele-specific presentation models with an RBM architecture by running the RBM-MHC algorithm [Bravi et al., 2021b] with default parameters (10 hidden units, $\lambda_1^2 = 0.001$). The peptide data used by default as training set by RBM-MHC to build allele-specific presentation models are downloaded from IEDB, see [Bravi et al., 2021b]. The RBM-MHC algorithm internally aligns peptide sequences of its training dataset to the reference length of 9 amino acids.
acids; we utilized the same alignment routine to align the immunogenic peptides against the seed given by the RBM-MHC training data.

The second step is training the diffRBM units on the 'selected' datasets, here immunogenic peptides and antigen-specific CDR3βs (section 5.2). Note that, for the applications discussed in this paper, the training on the background dataset and the one of the selected dataset are implemented consecutively. While the diffRBM architecture allows for a parallel training, the appeal to the sequential training here is justified by the size of the datasets at out disposal ($|D^b| \gg |D^s|$) and by the particular problem under consideration (modeling an enrichment of specific features, occurring on top of baseline features, which do not change as feedback from the former). The datasets $D^s$, collected as described in section 5.2, are divided into a training set with 80% of the data (used for training and for model selection) and a test set with the remaining 20% (used for model validation, see section 5.5), repeating this split 50 times.

Having fixed the background model, we used the largest available dataset of immunogenic sequences (the one for HLA-A*02:01) to perform model selection by cross-validation, as follows. We further divided randomly each of the 50 training sets into a set actually used for training and a validation set (with respectively 80% and 20% of the training set). We used this training/validation partitions to select optimal hyperparameters for the differential part (number of the number of hidden units and regularization penalty $\lambda_1^2$), by training diffRBM models on the training set at varying hyperparameters and monitoring the average diffRBM units' score (11) on the validation sets (Figure S3A-B). We also performed additional checks on the diffRBM units' AUC of immunogenic vs non-immunogenic discrimination with different hyperparameters (Figure S3C) and in a control case (Figure S3D).

For the model of T-cell specificity, the background is given by an RBM trained on a random subsample of $10^6$ CDR3β sequences either from the bulk repertoires in [Emerson et al., 2017] or the ones in [Dean et al., 2015] (section 5.2.2), choosing the optimal RBM architecture (100 hidden units, $\lambda_2^1 = 0.001$) by cross-validation (Figure S14A-B). For cross-validation, we used 50 random partitions of into $10^6$ data for training and $2 \times 10^5$ data for validation randomly chosen from the repertoire [Emerson et al., 2017] into training and validation sets (with respectively 80% and 20% of the data). The grid search for optimal hyperparameters for the differential part was then carried out at fixed background model, using the largest dataset of peptide-specific receptors (NLVPMVATV) and partitioning it 50 times at random into training and validation sets (with respectively 80% and 20% of the data, Figure S14C-F).

## 5.4 Alternative approaches

We compared diffRBM performance either to alternative architectures for differential models (diffRBM linear, PWM-based approach, SONIA) or to fully supervised, background-free models (classifiers).

### 5.4.1 DiffRBM linear.

DiffRBM linear is a diffRBM architecture where the differential part is specified by single-site fields $g_i^d(\sigma_i)$ (without the addition of weights $w_{i\mu'}^d(\sigma_i)$).

### 5.4.2 PWM-based approach.

We considered a type of differential model based entirely on the information in single-site amino frequency, where both the background model and the differential part are specified only by single-site fields ($g_i^b(\sigma_i)$ and $g_i^d(\sigma_i)$ respectively). The equations for learning $g_i^b(\sigma_i)$ and $g_i^d(\sigma_i)$ are:

$$ \frac{e^{g_i^b(\sigma_i)}}{\sum_{\sigma_i} e^{g_i^b(\sigma_i)}} + \lambda_2 g_i^b(\sigma_i) = f_i^b(\sigma_i) \qquad (13) $$

$$ \frac{e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))}}{\sum_{\sigma_i} e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))}} + \lambda_2(g_i^b(\sigma_i) + g_i^d(\sigma_i)) = f_i(\sigma_i) \qquad (14) $$

where $f_i^b(\sigma_i)$ is the frequency of amino acid $\sigma_i$ at position $i$ in the background dataset $D^b$, $f_i(\sigma_i)$ is the frequency of amino acid $\sigma_i$ at position $i$ in the dataset $D^s$ and $i = 1, ..., N$. We use the same $L_2$
acids; we utilized the same alignment routine to align the immunogenic peptides against the seed given by the RBM-MHC training data.

The second step is training the diffRBM units on the 'selected' datasets, here immunogenic peptides and antigen-specific CDR3βs (section 5.2). Note that, for the applications discussed in this paper, the training on the background dataset and the one of the selected dataset are implemented consecutively. While the diffRBM architecture allows for a parallel training, the appeal to the sequential training here is justified by the size of the datasets at out disposal (|D$^b$| ≫ |D$^s$|) and by the particular problem under consideration (modeling an enrichment of specific features, occurring on top of baseline features, which do not change as feedback from the former). The datasets D$^s$, collected as described in section 5.2, are divided into a training set with 80% of the data (used for training and for model selection) and a test set with the remaining 20% (used for model validation, see section 5.5), repeating this split 50 times.

Having fixed the background model, we used the largest available dataset of immunogenic sequences (the one for HLA-A*02:01) to perform model selection by cross-validation, as follows. We further divided randomly each of the 50 training sets into a set actually used for training and a validation set (with respectively 80% and 20% of the training set). We used this training/validation partitions to select optimal hyperparameters for the differential part (number of the number of hidden units and regularization penalty $\lambda_2^1$), by training diffRBM models on the training set at varying hyperparameters and monitoring the average diffRBM units' score (11) on the validation sets (Figure S3A-B). We also performed additional checks on the diffRBM units' AUC of immunogenic vs non-immunogenic discrimination with different hyperparameters (Figure S3C) and in a control case (Figure S3D).

For the model of T-cell specificity, the background is given by an RBM trained on a random subsample of 10⁶ CDR3β sequences either from the bulk repertoires in [Emerson et al., 2017] or the ones in [Dean et al., 2015] (section 5.2.2), choosing the optimal RBM architecture (100 hidden units, $\lambda_2^1 = 0.001$) by cross-validation (Figure S14A-B). For cross-validation, we used 50 random partitions of into 10⁶ data for training and 2 × 10⁵ data for validation randomly chosen from the repertoire [Emerson et al., 2017] into training and validation sets (with respectively 80% and 20% of the data). The grid search for optimal hyperparameters for the differential part was then carried out at fixed background model, using the largest dataset of peptide-specific receptors (NLVPMVATV) and partitioning it 50 times at random into training and validation sets (with respectively 80% and 20% of the data, Figure S14C-F).

## 5.4 Alternative approaches

We compared diffRBM performance either to alternative architectures for differential models (diffRBM linear, PWM-based approach, SONIA) or to fully supervised, background-free models (classifiers).

### 5.4.1 DiffRBM linear.

DiffRBM linear is a diffRBM architecture where the differential part is specified by single-site fields $g_i^d(\sigma_i)$ (without the addition of weights $w_{i\mu'}^d(\sigma_i)$).

### 5.4.2 PWM-based approach.

We considered a type of differential model based entirely on the information in single-site amino frequency, where both the background model and the differential part are specified only by single-site fields ($g_i^b(\sigma_i)$ and $g_i^d(\sigma_i)$ respectively). The equations for learning $g_i^b(\sigma_i)$ and $g_i^d(\sigma_i)$ are:

$$ \frac{e^{g_i^b(\sigma_i)}}{\sum_{\sigma_i} e^{g_i^b(\sigma_i)}} + \lambda_2 g_i^b(\sigma_i) = f_i^b(\sigma_i) \qquad (13) $$

$$ \frac{e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))}}{\sum_{\sigma_i} e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))}} + \lambda_2(g_i^b(\sigma_i) + g_i^d(\sigma_i)) = f_i(\sigma_i) \qquad (14) $$

where $f_i^b(\sigma_i)$ is the frequency of amino acid $\sigma_i$ at position $i$ in the background dataset D$^b$, $f_i(\sigma_i)$ is the frequency of amino acid $\sigma_i$ at position $i$ in the dataset D$^s$ and $i = 1, ..., N$. We use the same $L_2$
regularization penalty that we chosen for the RBM fields, and the parameter $\lambda_2$ sets its strength (like for the RBM-based approaches, $\lambda_2 = 1/|D^b|$ for background models and $\lambda_2 = 1/|D^s|$ for differential part). This approach is equivalent to learning from the background and selected datasets Position Weight Matrices (PWMs), probabilistic models of amino acid usage that treat all sequence positions as independent. By defining $P_i^b(\sigma_i) = e^{g_i^b(\sigma_i)} / \sum_{\sigma_i} e^{g_i^b(\sigma_i)}$ and $P_i(\sigma_i) = e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))/2} / \sum_{\sigma_i} e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))}$ and given the independent-site assumption of PWMs, the background data PWM probability and the selected data PWM probability are respectively recovered as $P^b(\sigma) = \prod_i P_i^b(\sigma_i)$ and $P(\sigma) = \prod_i P_i(\sigma_i)$. Since $P_i^b(\sigma_i)$ and $P_i(\sigma_i)$ are learnt to closely reproduce the frequencies $f_i^b(\sigma_i)$ and $f_i(\sigma_i)$, their predictions are the ones that are inferable from single-site amino acid frequency alone in the respective training datasets (hence in Figures 2, 5, S4, S5, S15 we refer to such predictions by ‘Amino acid frequency’). The prediction by the differential fields can be rewritten as:

$$g_i^d(\sigma_i) \sim \log \frac{P_i(\sigma_i)}{P_i^b(\sigma_i)} \qquad (15)$$

i.e., the differential parameters in this case simply corresponds to the log ratio of the PWM learnt on selected data and the PWM learnt on the background data, and measures the site-specific enrichment in amino usage in selected data compared to the background. For a sequence $\sigma$, the score based on the differential part is then simply given by the log-likelihood ratio:

$$\mathcal{L}^d(\sigma) = \sum_i \log \frac{P_i(\sigma_i)}{P_i^b(\sigma_i)} \qquad (16)$$

To estimate the entropy of PWM-based approaches (Figure S20C), we used expressions (12) with PWM probabilities for $P^b(\sigma)$ and $P(\sigma)$.

### 5.4.3 SONIA

We implemented a SONIA model using the software package soNNia from [Isacchini et al., 2021] and setting the option deep=False. For consistency with the background of diffRBM, we first learn a background model from $10^6$ randomly assembled from the universal donor repertoire from [Emerson et al., 2017], next we learn a model from epitope-specific repertoires and we take, as SONIA score, the difference between the scores assigned by these two models. The $L_2$-type regularization strength is set to the default value (zero); higher regularizations applied to one of our test sets did not provide visible improvement of the log-likelihood.

### 5.4.4 Classifier of immunogenic peptides

As a term of comparison for the immunogenicity model, we implemented first a linear (logistic) classifier. (see Model C1 in Table S2), which we trained by minimization of a binary cross-entropy loss. To find an optimal classifier, we also searched among different neural-network architectures trained to discriminate the immunogenic and non-immunogenic data retrieved from IEDB (section 5.2.1). We used the same 50 validation sets held out from the training set of immunogenic, HLA-A*02:01-presented peptides utilized for the model selection of diffRBM. In an analogous way, we randomly partitioned the set of non-immunogenic, HLA-A*02:01-presented peptides into training/test sets (with respectively 80%/20% of the data), repeating the partition 50 times, and we further partitioned the training sets to hold out at each repetition a validation set with 20% of sequences. We trained the models varying the number and width of hidden layers (Table S2), as well as the weight decay. The ‘optimal classifier’ performance for HLA-A*02:01 (Figure 4A) is obtained by selecting for each partition the best performing architecture on the validation set and evaluating its AUC on the test set. Next, we identified the architecture giving the maximal average AUC of discrimination between the immunogenic and non-immunogenic peptides in the validation sets (Model C8 in Table S2, see Figure S11) and we used it to estimate the ‘optimal classifier’ average AUC for the HLA-B*07:02 and HLA-B*35:01 models (Figure 4C-D). Our numerical implementation is based on the PyTorch library. Training is performed in mini-batches of 64 sequences,
regularization penalty that we chosen for the RBM fields, and the parameter $\lambda_2$ sets its strength (like for the RBM-based approaches, $\lambda_2 = 1/|D^b|$ for background models and $\lambda_2 = 1/|D^s|$ for differential part). This approach is equivalent to learning from the background and selected datasets Position Weight Matrices (PWMs), probabilistic models of amino acid usage that treat all sequence positions as independent. By defining $P_i^b(\sigma_i) = e^{g_i^b(\sigma_i)} / \sum_{\sigma_i} e^{g_i^b(\sigma_i)}$ and $P_i(\sigma_i) = e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))/2} / \sum_{\sigma_i} e^{(g_i^b(\sigma_i)+g_i^d(\sigma_i))}$ and given the independent-site assumption of PWMs, the background data PWM probability and the selected data PWM probability are respectively recovered as $P^b(\sigma) = \prod_i P_i^b(\sigma_i)$ and $P(\sigma) = \prod_i P_i(\sigma_i)$. Since $P_i^b(\sigma_i)$ and $P_i(\sigma_i)$ are learnt to closely reproduce the frequencies $f_i^b(\sigma_i)$ and $f_i(\sigma_i)$, their predictions are the ones that are inferable from single-site amino acid frequency alone in the respective training datasets (hence in Figures 2, 5, S4, S5, S15 we refer to such predictions by ‘Amino acid frequency’). The prediction by the differential fields can be rewritten as:

$$g_i^d(\sigma_i) \sim \log \frac{P_i(\sigma_i)}{P_i^b(\sigma_i)} \qquad (15)$$

i.e., the differential parameters in this case simply corresponds to the log ratio of the PWM learnt on selected data and the PWM learnt on the background data, and measures the site-specific enrichment in amino usage in selected data compared to the background. For a sequence $\sigma$, the score based on the differential part is then simply given by the log-likelihood ratio:

$$\mathcal{L}^d(\sigma) = \sum_i \log \frac{P_i(\sigma_i)}{P_i^b(\sigma_i)} \qquad (16)$$

To estimate the entropy of PWM-based approaches (Figure S20C), we used expressions (12) with PWM probabilities for $P^b(\sigma)$ and $P(\sigma)$.

### 5.4.3 SONIA

We implemented a SONIA model using the software package soNNia from [Isacchini et al., 2021] and setting the option deep=False. For consistency with the background of diffRBM, we first learn a background model from $10^6$ randomly assembled from the universal donor repertoire from [Emerson et al., 2017], next we learn a model from epitope-specific repertoires and we take, as SONIA score, the difference between the scores assigned by these two models. The $L_2$-type regularization strength is set to the default value (zero); higher regularizations applied to one of our test sets did not provide visible improvement of the log-likelihood.

### 5.4.4 Classifier of immunogenic peptides

As a term of comparison for the immunogenicity model, we implemented first a linear (logistic) classifier. (see Model C1 in Table S2), which we trained by minimization of a binary cross-entropy loss. To find an optimal classifier, we also searched among different neural-network architectures trained to discriminate the immunogenic and non-immunogenic data retrieved from IEDB (section 5.2.1). We used the same 50 validation sets held out from the training set of immunogenic, HLA-A*02:01-presented peptides utilized for the model selection of diffRBM. In an analogous way, we randomly partitioned the set of non-immunogenic, HLA-A*02:01-presented peptides into training/test sets (with respectively 80%/20% of the data), repeating the partition 50 times, and we further partitioned the training sets to hold out at each repetition a validation set with 20% of sequences. We trained the models varying the number and width of hidden layers (Table S2), as well as the weight decay. The ‘optimal classifier’ performance for HLA-A*02:01 (Figure 4A) is obtained by selecting for each partition the best performing architecture on the validation set and evaluating its AUC on the test set. Next, we identified the architecture giving the maximal average AUC of discrimination between the immunogenic and non-immunogenic peptides in the validation sets (Model C8 in Table S2, see Figure S11) and we used it to estimate the ‘optimal classifier’ average AUC for the HLA-B*07:02 and HLA-B*35:01 models (Figure 4C-D). Our numerical implementation is based on the PyTorch library. Training is performed in mini-batches of 64 sequences,