Our data consist of the sequenced T-cell repertoires of 4,144 subjects with HLAs genotyped via next generation sequencing (NGS) (47) (see Supplementary Figure S1 for demographic distributions). The median number of unique T cells sequenced from each individual is ~ 227,000 and 90% of subjects have counts between ~ 74,000 and ~ 610,000. The majority of our samples are taken from healthy adults residing in the United States; ~ 5% and ~ 20% are Lyme and Covid positive, respectively. For a given HLA, we separate subjects into cases and controls defined as those with and without the HLA, respectively. Subjects expressing an HLA that is in the same p-group ¹ as the HLA of interest are excluded from the control group, as such HLAs have identical amino acid sequences in the peptide binding region and thus may share TCR specificity. HLA-DP and -DQ are treated as heterodimers, with cases and controls defined by α-β pairs. The α chain of HLA-DR is invariant and thus we treat these HLAs as monomers, similar to class I HLAs. We randomly select a fixed 80% and 20% of the samples for training and validation, respectively, and evaluate model generalization on an entirely independent cohort from Kanazawa, Japan.

We identify sets of public HLA-associated TCRβs using a statistical approach that enforces the assumption that any TCRβ is associated with at most one HLA allotype (we test this assumption below). This assumption enables us to disentangle the effects of linkage disequilibrium (LD) among HLA loci (48) that would otherwise result in a large number of spurious HLA-TCR associations ( Supplementary Figure S2A ). We use exact matching of the TCRβ V-gene, J-gene and CDR3 to identify sequences thatare over-represented in subjects with a given HLA allotype. Thus, our association of TCRβs with HLAs is agnostic to the specific amino acid sequence, it solely relies on it being observed in multiple repertoires.

Our assumption that TCRs typically associate with a single HLA allotype is important for resolving LD but is also supported by current biological evidence. While TCRs may be theoretically capable of broad cross-reactivity (42), there is limited experimental evidence that public TCRs functionally recognize multiple unrelated pHLAs. Rather, evidence suggests that functional cross-reactivity may be rare, with TCRs typically exhibiting specificity for particular epitopes and HLAs (49–51). Nonetheless, the true extent of functional cross-reactivity in physiologically relevant settings remains an open question. Our sequence identification procedure is primarily powered to detect public TCRβs associated with a single HLA and has limited sensitivity to detect sequences with broad cross-reactivity, if such sequences exist.

The identification procedure works as follows (see Methods for precise details): for each HLA allotype, we first create a set of candidate HLA-associated TCRβs using a one-sided Fisher’s Exact Test (FET) to identify TCRβs over-represented in cases. We adopt a pre-specified fiducial p-value threshold, p ^*. For each unique candidate TCRβ, we fit an L1-regularized Logistic Regression (L1LR) model which predicts the presence of that TCRβ in subjects given their HLAs (represented as abinary vector of indicator variables). We tune the L1 hyperparameter λ to be the smallest value for which exactly one HLA parameter is non-zero. In other words, we determine which single HLA allotype best predicts the observed distribution of a given TCRβ in the repertoires of our training sample. We test all TCRβs with p-values < p ^* and retain only TCRβs which associate most strongly with the HLA being modeled. As our interest here is primarily in the characteristics of titsets of HLA-associated TCRβs, we set p ^* to a permissive value of p ^* = 10^-4 and use the hold-out repertoires for validation. Note that due to exclusion of p-group matched HLAs a small number of TCRβs are assigned to multiple HLAs due to variations in the training data ( Supplementary Figure S2B ).

We associate ~ 10⁶ TCRβs to specific HLAs, for a median of 2,400 TCRβs per HLA with ~ 70% of HLAs having a total number of associated TCRβs in the range of 1,600-5,000. To maintain consistency with other work associating TCRβs with disease (37, 38, 52, 53), we refer to these HLA-associated TCRβs as enhanced sequences (ES).

HLA homozygosity has been epidemiologically linked to poor clinical prognosis in the context of both chronic HIV infection (56) and cancer checkpoint-inhibitor immunotherapy (57), possibly due to the reduced size of the HLA-restricted peptide repertoire available for T-cell recognition. This reduced size of the antigen repertoire, coupled with higher relative surface concentration of pHLAs associated with the homozygous protein, may directly impact the TCR repertoire by increasing the probability of clonal expansion of T cells expressing cognate TCRs. Consistent with this hypothesis, we find that the distribution of ES counts is elevated among homozygous individuals ( Figures 3A-F ). Across all HLAs, the distribution of ESs is(on average) about one standard deviation higher for homozygous as compared to heterozygous individuals ( Figures 3G, H ). Notably, the ES distribution is an additional standard deviation higher among individuals homozygous at both the α and β locus for HLA-DP or -DQ (“double homozygous”, Figure 3H ).

While these results are consistent with increased relative antigen abundance increasing clonal expansion of associated T cells, an alternative hypothesis is that a decrease in the relative surface expression or decreased antigenic diversity in homozygous subjects results in less crowding out by other HLAs or TCRs, respectivly. This crowding hypothesis implies an increased breadth across multiple loci within a given class for homozygous subjects. We test this hypothesis by examining whether a homozygous subject at one locus has higher breadth at another class-matched loci. For example, if HLA and/or TCRs crowd each other, we expect that a subject homozygous for HLA-A will also have a higher average breadth in their HLA-B and/or HLA-C response due to lower diversity at the class I loci. We do not find evidence of such an effect (see Supplementary Figure S10 ) and thus conclude that crowding out by HLAs and/or TCRs may not be the primary driver of increased breadth unless it is restricted to a particular locus.

Taken together, these results suggest homozygosity at a particular locus increases the breadth of the T-cell response against peptides presented by that HLA, possibly through increased surface expression and antigen presentation.

The clear separation of ES counts in HLA cases versus controls implies that HLA allotypes can be easily imputed from HLA-associated TCR β s alone. To this end, we fit a simple logistic regression model for each HLA allotype observed in at least 30 training samples (representing ∼ 1 % expression frequency), predicting whether an individual expresses that HLA as a function of observed ES and total-unique-rearrangement (log) counts (see Figures 4A, B for representative examples, Supplementary Figures S3 - S8 for all HLAs tested). The number of subjects with each HLA allotype based on HLA genotyping and the model performance for each imputed allotype are provided in Table 2 (class I) and Table 3 (class II); model performance is also shown in Figure 4C .

Over all the imputation accuracy is extremely high with area under the receiver operating characteristic curve (AUC-ROC) scores of ≥ 0.9 for all but 3 of the 120 HLA-A,-B, and -DR allotypes modeled. This accuracy highlights the specificity of HLA ESs,indicating HLAs at these loci can be accurately imputed from immunosequencing alone. Furthermore, this accuracy confirms that HLA β alone is typically sufficient to identify shared antigen specificity of public T cells. Among class I HLAs, HLA-C allotypes are a notable outlier with relatively lower classification performance even among models with a significant number of positive training examples ( Figure 4D ). We hypothesize that this reduced performance is due in part to the ∼ 10 × -lower surface expression of HLA-C compared to allotypes expressed by other class I genes (58).

Model performance among HLA-DP and -DQ heterodimers is substantially more variable than performance at the other loci. We find that 8 of the 30 HLA-DP and 37 of 81 of HLA-DQ fail to achieve AUC-ROC scores of > 0.9 even among heterodimers expressed in a high number of individuals in our training population ( Figure 4E ). Notably, HLA-DQ and -DP are the only loci we study with polymorphic α chains, suggesting that heterodimer incompatibility may explain the lack of associated TCRs and the corresponding inability to impute expression of these heterodimers. We explore this issue in the next section.

We assess the generalizability of our imputation models and rule out overfitting or data leakage by evaluating model performance on an entirely independent, out-of-distribution cohort of 136 individuals from Kanazawa, Japan. Subjects in this cohort have HLAs genotyped by direct sequencing and none are included in our training or internal holdout sets. This population differs both genetically and geographically from the predominantly U.S.-based training cohort, thus providing a stringent test of the model’s generalizability to out-of-distribution data. For this analysis, we use a set of 135 HLAs across all six loci that yield models with cross-validation precision > 0.9 , recall > 0.8 and ≥ 30 positive training cases.

We apply the 135 HLA models to the Japanese cohort, resulting in an average of 8.5 imputations per individual, compared to 11.0 per individual in our internal validation cohort. This difference is expected and reflects the shift in HLA allele frequencies between the Japanese and U.S-based cohorts. Despite these differences, we observe comparable model performance for the subset of overlapping HLAs. The overall F1 score ³ for the Japanese cohort, aggregated across all imputations, is statistically consistent ( < 3 σ ) to that observed in the internal validation cohort ( 0.925 ± 0.006 and 0.940 ± 0.002 , respectively). These results confirm that model performance generalizes across genetic and geographic backgrounds and that the high accuracy we observe is not a consequence of data leakage, but rather reflects a strong HLA-specific signal in TCR repertoires.

HLAs are inherited as a haplotype, such that one full set is inherited on a single chromosome from each parent (59). The α and β chains of HLA-DP and -DQ ⁴ pair after synthesis, yielding a phenotype of up to four unique heterodimers in each individual: two formed in cis, where both subunits are encoded on the same chromosome, and two in trans, where the subunits are encoded on opposite chromosomes. Given the high degree of polymorphism observed in HLA-DP and -DQ subunits, it is perhaps unsurprising that some pairs of α and ta chains do not form stable heterodimers (60–62). Based on structural and sequence analysis of HLA-DQ, Tollefsen et al. (62) propose specific group pairings that likely form stable heterodimers (though they note a small number of exceptions to the pairing rules).

In the context of the present study, the proposed existence of incompatible (and thus non-functional) α and β chains implies two testable hypotheses: (1) that co-inheritance of incompatible pairs on the same chromosome will be under strong negative selection (62); and (2) that incompatible pairs are unable to elicit a T-cell response and thus will not be associated with any public TCRβs.

The first hypothesis implies that co-expression of incompatible pairs almost always results from trans-complementation; as such, incompatible pairs will be in linkage equilibrium since they are not co-inherited. Conversely, cis-complementation should necessarily result in functional pairs thus all pairs forming from subunits in linkage disequilibrium should be functional. For two subunits α and β, with respective expression frequencies f α and f β and co-expression frequency f α + β , linkage equilibrium results in an expected co-expression frequency of approximately E L E [ f α + β ] = 2 f α f β (63). Overall, we find that a majority (75 of 119) HLA-DP and -DQ alpha + β pairs appear to be in linkage equilibrium, i.e., not co-inherited ( Figure 5A ; only subunits comprising heterodimers observed in at least 30 individuals are considered). Notably, all of the pairs that Tollefsen et al. (62) propose to be incompatible are in equilibrium ( Figure 5B ).

Given the tight linkage between genes encoding subunits, any pair expressed in cis in at least one individual is likely to be in LD in a large cohort. Consistent with this hypothesis, every individual heterozygous at both the α and a locus has at least two of four α/β pairs in LD; conversely, no individual should have more than two α/β pairs in equilibrium. We confirm that no subject in our cohort has more than two heterodimers formed from subunits in linkage equilibrium, as expected. Taken together, we conclude that LD is a strong proxy for cis-complementarity, and that, given strong selection pressure, incompatible pairs are only expressed in trans.

If incompatible pairs are truly non-functional (with respect to antigen presentation and T-cell recognition), then there should not be any TCRβs that are specific to such pairs. As an example, we are unable to identify distinguishing TCRβs forDQA1*01:02+DQB1*03:01, which both violates the Tollefsen et al. (62) pairing rules and is in linkage equilibrium ( Figure 5C ). Moreover, all of our high-frequency, poor-performing HLA-DP and -DQ imputation models are for heterodimers formingfrom subunits that are in linkage equilibrium, and thus are almost certainly expressed in trans ( Figure 5D ; see also Figure 3C ). Moreover, almost all pairs that violate the Tollefsen et al. (62) pairing rules have lower-than-expectedimputation performance ( Figure 5D ).

Notably, while heterodimer incompatibility implies both linkage equilibrium and poor model performance, not all heterodimers forming from subunits in linkage equilibrium result in poor performing models. Thus, trans-complementation may yield heterodimers which can drive a T-cell response resulting in identifiable TCRβs and a high-accuracy imputation model ( Figure 5D ).

Using these results on model performance and gene linkage, we can extend the Tollefsen et al. (62) pairing rules to the HLA-DP locus: DPA1*02 appears to form unstable heterodimers when paired with DPB1*02 and DPB1*04; DPA1*01 is apparentlyunrestricted, forming stable heterodimers with subunits from all DPB1 groups in our sample.

The T-cell repertoire consists of a mixture of naive and memory T cells. In principle, HLA-restricted antigen presentation will bias both thymic selection and clonal expansion, and thus HLA-specific signatures may exist in both compartments. However, ourstatistical approach to identifying HLA-specific public TCRβs likely favors identification of TCRs from memory T-cells.

We investigate characteristics of HLA associated sequences by comparing the distribution of ESs observed in the memory and naive compartments sequenced from 45 individuals who were not included in the original study. For each subject, we sequence five separate repertoires: the memory and naive compartments of CD8⁺ and CD4⁺ T cells, respectively, and the unsorted repertoire. As sequence-based typing was unavailable for these individuals, we treat the imputed HLAs from the the unsorted repertoire as ground truth (limiting to 90 models with >0.95 precision and recall). For each repertoire and each possible HLA, we compute the weighted breadth of the HLA-specific ESs observed in the repertoire. Across all 45 unsorted repertoires, the weighted breadth of both class I and class II ESs is substantially higher for the HLAs an individual expresses compared to those they do not express ( Figure 6A ; compare HLA-positive to -negative).

Within the sorted compartments, a striking pattern emerges: in the naive compartments, there is little difference in the breadth of ES specific for an individual’s expressed HLAs compared to background ( Figures 6C, E ), while in the memory compartments, ESs specific to the individuals’ expressed HLAs have far higher breadth ( Figures 6B, D ). Moreover, within the memory compartment, CD8⁺ cells are closely linked to high relative breadth of class I HLA ESs, while CD4⁺ cells are closely linked to high relative breadth of class II HLA ESs. This result provides further confirmation that ESs are correctly mapped to HLAs despite the challenges of HLA LD.

The centrality of the memory compartment in driving our HLA imputation signal raises several additional hypotheses. The first is that clonal expansion increases the likelihood of detection for ESs. This increased likelihood implies that, while ESs are frequently observed in individuals without the associated HLA, they will tend to be at higher repertoire frequency among individuals who do express the HLA. Indeed, we observe a notable increase in the distribution of clonal frequency among cases compared controls ( Figure 7A ).

The second hypothesis is that clonal expansion results from antigen exposure, which will also result in polyclonal expansion of T cells that express distinct TCRs that all respond to the same antigen. An extreme form of polyclonal expansion is convergent recombination, in which multiple, distinct TCR DNA sequences encode identical amino acid sequences. Consistent with this hypothesis, the per-ES convergent recombination rates are higher when ESs are observed in cases as compared to controls ( Figure 7B ; conditional on the ES being present in the repertoire).

The publicity of HLA ESs combined with their apparent antigen-specific clonal expansion suggests that these sequences are responding to antigens which are common in the human population (64). We test this hypothesis in the context of SARS-CoV-2 exposure. Because of its novel nature, SARS-CoV-2 is especially well-suited for this task as we are able to confidently assign Covid-19 negativity to samples collected before 2020. In our training sample for deriving HLA restricted sequences, 694 ofthe 4,144 subjects (~ 20%) are Covid-19 positive based on PCR labels; the remaining samples are collected before 2020 and thus Covid-19 negative. Because SARS-CoV-2 exposure is relatively high in our training sample, we expect that some of the HLA-specific ESs we identity are SARS-CoV-2 specific.

We identify 6866 SARS-CoV-2-specific ESs using a set of 1523 SARS-CoV-2 PCR-positive samples that have no overlap with the typed HLA training samples and 4386 controls (1008 controls overlap with the typed HLA samples). A high proportion (80% of the most confident SARS-CoV-2 specific ESs) are also HLA-specific ESs ( Figure 7C , blue line). This overlap in ESs results from the fact that ~ 20% of our typed HLA training samples are SARS-CoV-2 positive, i.e., SARS-CoV-2 specific ESs are predictiveof HLAs as well. In contrast, if we define an alternative set of HLA-ESs using only 3,450 repertoires which were sampled prior to 2020, we find very little overlap ( Figure 7C , orange line). The limited overlap we do observe may be due to cross-reactivityto homologous epitopes from other coronaviruses and/or false positive sequences in the SARS-CoV-2/HLA ES set.

The intersection of SARS-CoV-2 specific ES set with the and HLA specific ES sets yields 1,880 TCRβs (using a threshold of p * = p < 10 − 4 ) that are associated with a particular HLA in the context of SARS-CoV-2 infection. Thus, these ESs are almost certainly specific to commonly targeted SARS-CoV-2 T-cell epitopes. To confirm this specificity, we impute HLAs for samples in the SARS-CoV-2 training data ⁵ and then compute the fraction of HLA-associated SARS-CoV-2 ES observed in their repertoire. For “HLA +” and “HLA -” subjects we only count sequences if the subject has or does not have the HLA which the SARS-CoV-2 ES is associated with, respectively. Notably, both class-I and class-II associated TCRβs are far more commonly observed in individuals with the restricting HLA and known SARS-CoV-2 infection ( Figure 7D ), thus confirming the HLA- and pathogen-specificity of these TCRβs and further demonstrating that while TCRβs may be cross-reactive to many distinct antigens, the likelihood of cross-reactivity in vivo is low (65).

We conclude that, taken together, these results demonstrate that the majority of HLA-specific ESs are the result of clonal expansion of T cells responding to common antigens, thus establishing an immunological foundation for HLA imputation from TCRβrepertoires.

TCRβ and HLA sequence data from human samples used for these studies were aggregated from several independent study collections described below. All necessary patient/participant consent has been obtained for each study and the appropriate institutional forms have been archived.

We use the ImmunoSEQ assay developed by Adaptive Biotechnologies to measure the TCRβ repertoire. ImmunoSEQ is a multiplexed PCR-based method targeting rearranged TCRβ sequences (45). The assay uses a panel of primers specific to all functional TRBV and TRBJ gene segments to amplify the full V(D)J region, including part of the TRBV segment, the entire diversity (D) and a portion of the J segment. The sequenced portions of the V and J segments are sufficient to determine gene usage. This region also includes the complete hypervariable complementarity-determining region 3 (CDR3), which accounts for most of the TCR diversity and is the key determinant of TCR antigen specificity. The reverse primer binds near the 3′ end of the J region, adjacent to the constant region and may extend into the 5′ end of the TCR. PCR primers contain universal priming sites on the 3′-end and Illumina sequencing adaptors on the 5′-end. This design enables robust amplification across a wide range of rearrangements while controlling for bias. High-throughput sequencing of these amplicons provides accurate quantification of unique TCRβ clonotypes.

Flow cytometry and cell sorting 3 − 5 × 10 7 PBMC were stained using a cocktail of antibodies that included CD3 (clone UCHT1, Biolegend), CD4 (clone OKT4, Biolegend), CD8 (clone RPA-T8, Biolegend), CD45RA(clone HI100, BD Bioscience), CCR7 (clone G043H7, Biolegend), CD95 (clone DX2, Biolegend), and CD28 (clone CD28.2, Biolegend) for 10 minutes at 4 degrees C. PBMC were washed with MACS buffer (Miltenyi Biotec) and then PE+ CD3+ T cells were enriched using anti-PE Microbeads with LS columns (Miltenyi Biotec) following the manufacturer’s protocol. The enriched CD3 sample was sorted using a FACSAria Fusion Flow Cytometer (BD Biosciences) to isolate naive CD4 (CD4⁺ CD3⁺ CD45RA⁺ CD95^– CD28⁺ CCR7⁺), memory CD4 (CD4⁺ CD3⁺ non-naive CD4), naive CD8 (CD8⁺ CD3⁺ CD45RA⁺ CD95^– CD28⁺ CCR7⁺) and memory CD8 (CD8⁺CD3⁺ non-naive CD8). Between 150,000 and 3 × 10 6 T cells were sorted and sent in for sequencing for each subset. TCRβ sequencing was carried out using the ImmunoSEQ assay at Adaptive Biotechnologies.

We independently train a model for each HLA using the following procedure:

Below we describe training steps in detail using the following definitions: Let H be the set of N_H typed HLA allotypes in our training set and h i ∈ H , i = 1 , … , N H refer to any single HLA in the set. The training set of h_i is denoted as t i ∈ T where T is the set of N_T training samples. We similarly define a test (or holdout) set for each HLA t i ′ ∈ T ′ where T′ is the set of N T ′ holdout samples. Below we drop the subscript i unless needed for clarity as the procedure is applied to each HLA independently.

Step 1: T and T′ are derived from a random 80/20 split of all typed data with N T + N T ′ = 4 , 144 . Supplementary Figure S1 shows the age, sex and ethnicity distribution of sample subjects.

Step 2: When building a model for a given HLA h i , samples are labeled as cases (label = 1) or controls (label = 0) if a subject expresses or does not express h i , respectively. A subset of controls are unlabeled if a p-group matched HLA to h i is expressed by the subject. P-group matched HLAs share the same antigen binding domain and thus may share T-cell receptor (TCR) specificities.

Step 3: For any HLA h, an initial set of ESs, E, is defined using Fisher’s Exact Test (FET) [74] applied to training set t of HLA h. The a range of p-value thresholds are used to identify ESs and the threshold is treated as ahyperparameter in the modeling; it is independently derived for each HLA we model.

Step 4: If HLA h is in linkage disequilibrium (LD) with other HLAs, the enhanced sequences specific to HLAs in LD will be present in the ES set of HLA h because of their co-occurrence in subjects (this is how LD is defined). The goal of our procedure is to identify and remove these sequences from the final set of ESs associated to h. Thus, a fundamental assumption of the model is that for a given train set, an ES may be associated with only a single HLA.

We enforce the assumption that each TCRβ is specific to one HLA by fitting a logistic regression model with L1 regularization and we refer to this as the L1LR method. We independently associate each ES in E with a single HLA. Thus, the assumption that a given TCRβ is associated to a single HLA is not globally enforced, i.e., it is possible that due to variations in training data that the same TCRβ may be associated to multiple HLAs though this is rare (see Supplementary Figure S2 ).

When associating individual ESs to HLAs, we model the presence/absence of any given ES e j ∈ E in the training set t as a logistic regression classification. Here j denotes the individual TCRβs in the ES set. For any sample t l ∈ t , where l denotes individual samples in the training set t, the label for sequence e j is either 0 or 1 if the sequences is present or absent in t l , respectively. The feature vector for associating a given ES e j includes the indicator vector for the typed set of HLAs, i.e. I k l = 1 if sample l expresses HLA h k and I k l = 0 if sample l does not express HLA h k . Here, k indicates the subset of HLAs in H is in LD with HLA h (defined by FET p-value < 10 − 3 ). Thus, the feature vector is a binary vector indicating whether a subject has or does not have a given HLA which is observed to co-occur with the HLA being modeled. When associating sequences we include the total number of unique rearrangements in each sample repertoire, d l as a covariate not subject to L1 regularization. We model each ES independently and thus drop the subscript j in the following for clarity. Taken together a single ES e is modeled such that.

In Equation 1, α and β are free parameters, b is the bias term and λ is the regularization strength applied to β. We find the best-fit parameters by minimizing the log-loss and tune λ to be the smallest value that yields precisely one non-zero coefficient in β k . If the non-zero coefficient is β k = i , then the sequence is associated with the HLA being modeled and is retained as part of the final ES set. If the non-zero coefficient is β k ≠ i . the sequences is excluded from the final ES set of HLA h which is denoted as E ˜ . We weight each sample by the square root of the convergent recombination count when finding the best-fit parameters and zero count samples are given a weight of unity. The effectiveness of our procedure to resolve LD is demonstrated in Supplementary Figure S2 .

Step 5: We derive a per-sequence weight, w j for each sequence e j ∈ E ˜ by fitting a two feature logistic regression classifier. Here j refers to individual sequences in E ˜ . Each sequence is fit independently. For each sequence the target of the classifier is whether the sample is a case or control of HLA h i and the features are the convergent recombination count of the sequence and total number of unique rearrangements in the given repertoire sample. Essentially, we build a model to discriminate cases from controls using the count of a single sequence and take the coefficient as a weight. In detail, we standardize the features by subtracting the mean and normalizing to the standard deviation of the features across all samples. We fit the default logistic regression model in the Scikit-learn library which includes a fixed amount of L2 regularization ( λ = 1 ). This procedure yields a weight for each sequence in e j ∈ E ˜ .

Step 6: The final classifier model is a two feature model where the features F for subject l are.

In Equation 2 c l j is the convergent recombination count of e j in sample l. These features, F l , are used in a the standard Scikit-learn logistic regression to predict the presence/absence of an HLA for a given sample.

Step 7: To model any given HLA, we fit the standard Scikit-Learn logistic regression classifier using features generated in Step 6). We treat the p-value cutoff for identifying ESs via the FET as described in Step 3) as a hyperparameter of each HLA model. We train and evaluate a model using a five-fold cross validation strategy for p-value cutoffs = [ 10 − 3 , 10 − 4 ,   10 − 5 , 10 − 6 , 10 − 7 , 10 − 8 ] . We adopt p-value cutoff of the highest performing cross-validation model.

Step 8: We generate our final model by adopting the p-value threshold determined in step 7) and train on the full set of training data. We evaluate the final model on the holdout set ( Supplementary Figure S3 – S8 ).