Of the 80 million unique TCR$\beta$ chains (defined by V-gene family and CDR3 sequence) in the 666 cohort repertoires, about 11 million chains are found in at least two individuals and referred to here as public chains (for a more nuanced examination of TCR chain sharing see [Elhanati et al., 2018]). The occurrence patterns of these public TCR$\beta$s—the subset of subjects in which each distinct chain occurs—can be thought of as forming a very large binary matrix $M$ with about 11 million rows and 666 columns. Entry $M_{i,j}$ contains a one or a zero indicating presence or absence, respectively, of TCR $i$ in the repertoire of subject $j$ (ignoring for the moment the abundance of TCR $i$ in repertoire $j$; Figure 1 depicts two illustrative sub-matrices of $M$). (Emerson et al., 2017) demonstrated that this binary occurrence matrix $M$ encodes information on subject genotype and immune history: they were able to successfully predict HLA-A and HLA-B allele type and CMV serostatus by learning sets of public TCR$\beta$ chains with occurrence patterns that were predictive of these features. Specifically, each feature—such as the presence of a given HLA allele (e.g. HLA-A*02:01) or CMV seropositivity—defines a subset of the cohort members positive for that feature, and can be encoded as a vector of 666 binary digits. This phenotype occurrence pattern of zeros and ones can be compared to the occurrence patterns of all the public TCR$\beta$ chains to identify similar patterns, as quantified by a p-value for significance of co-occurrence across the 666 subjects; thresholding on this p-value produces a subset of significantly associated TCR$\beta$ chains whose collective occurrence in a repertoire was found by Emerson et al. to be predictive of the feature of interest (in cross-validation and, for CMV, on an independent cohort). Generalizing from these results, it is reasonable to expect that other common immune exposures may be encoded in the occurrence matrix $M$, and that these encodings could be discovered if we had additional phenotypic data to correlate with TCR occurrence patterns. In this study, we set out to discover these encoded exposures de novo, without additional phenotypic correlates, by learning directly from the structure of the occurrence matrix $M$ and using as well the sequences of the TCR$\beta$ chains (both their similarities to one another and to TCR sequences characterized in the literature). We hypothesized that patterns of TCR co-occurrence (correlations between rows in the matrix $M$) might indicate shared responses to unknown immune exposures, that co-occurrence between TCR chains and HLA alleles (correlations between rows in $M$ and rows in the HLA allele occurrence matrix) could be used to help identify functional TCR chains, and that clustering TCR$\beta$ chains by co-occurrence and sequence could highlight functional associations (Figure 1). To support this effort we assembled additional HLA typing data for the subjects, now at 4-digit resolution (e.g., A*02:01 rather than A*02) and including MHC class II alleles, and we compiled a dataset of annotated TCR$\beta$ chains by combining online TCR sequence databases, structurally characterized TCRs, and published studies (see Materials and methods; [Shugay et al., 2018; Tickotsky et al., 2017; Berman et al., 2000; Dash et al., 2017; Glanville et al., 2017; Song et al., 2017; Kasprowicz et al., 2006]). Here we describe the outcome of this discovery process, and we report a number of intriguing general observations about the role of HLA in shaping the T cell repertoire.

Figure 1

Download asset Open asset

The results of our analysis are organized in the remaining five sections as follows. We begin with an examination of TCR co-occurrence patterns across the full cohort (first section, Figures 2–3). In the next section we examine patterns of TCR-HLA association (Table 1 and Figures 4–5). In the third section we analyze TCR co-occurrence within subsets of the cohort positive for specific HLA alleles, and we identify TCR clusters that may be reflective of shared immune exposures (Figures 6–8). In the fourth section we use our dataset of HLA-associated TCR$\beta$ chains to identify covariation between HLA and the TCR$\beta$ CDR3 sequence (Table 2 and Figure 9). In the final section we focus on CMV-responsive TCR$\beta$ chains, examining their degree of HLA-restriction and the extent to which they may be responding to other antigens (Figure 10). Figure 1 provides a graphical overview of the co-occurrence analysis.

We hypothesized that we could identify unknown immune exposures encoded in the public repertoire by comparing the occurrence patterns of individual TCR$\beta$ chains to one another. A subset of TCR$\beta$ chains that strongly co-occur across the 666 cohort subjects might correspond to an unmeasured immune exposure that is common to a subset of subjects. Since shared HLA restriction could represent an alternative explanation for significant TCR co-occurrence, we also compared the TCR occurrence patterns to the occurrence patterns for class I and class II HLA alleles. We began by analyzing TCR occurrence patterns over the full set of cohort members. For each pair of public TCR$\beta$ chains $t_1$ and $t_2$ we computed a co-occurrence p-value $P_{\text{CO}}(t_1,t_2)$ that reflects the probability of seeing an equal or greater overlap of shared subjects (i.e., subjects in whose repertoires both $t_1$ and $t_2$ are found) if the occurrence patterns of the two TCRs had been chosen randomly (for details, see Materials and methods and Figure 12). In a similar manner we computed, for each HLA allele $a$ and TCR $t$, an association p-value $P_{\text{HLA}}(a,t)$ that measures the degree to which TCR $t$ tends to occur in subjects positive for allele $a$. Finally, for each pair of strongly co-occurring ($P_{\text{CO}}<1\times {10}^{-8}$) TCR$\beta$ chains $t_1$ and $t_2$, we looked for a mutual HLA association that might explain their co-occurrence, by finding the allele having the strongest association with both $t_1$ and $t_2$, and noting its association p-value:

where $A$ denotes the set of all HLA alleles. In words, we take the p-value of the strongest HLA allele association with the TCR pair, where the association of an HLA allele with a TCR pair is defined by the weakest association of the allele among the individual TCRs.

Based on this analysis, we identified two broad classes of strongly co-occurring TCR pairs (Figure 2): those with a highly significant shared HLA association, where the co-occurrence of the two TCRs can be explained by a shared HLA allele association (i.e. a common HLA restriction), and those with only modest shared HLA-association p-value, for which another explanation of co-occurrence must be sought. Points above the dashed $y=x$ line correspond to pairs of TCRs for which there exists an HLA allele whose co-occurrence with each of the TCRs is stronger than their mutual co-occurrence, while for points below the line no such HLA allele was present in the dataset.

Figure 2

Download asset Open asset

Figure 2—source data 1

TCR pairs and corresponding $P_{\text{CO}}$ and $P_{\text{HLA}}$ values.

https://doi.org/10.7554/eLife.38358.005

Download elife-38358-fig2-data1-v2.tds

We used a neighbor-based clustering algorithm, DBSCAN (Ester et al., 1996), to link strongly co-occurring TCR pairs together to form larger correlated clusters (see Materials and methods), and then investigated phenotype associations with these clusters. At an approximate family-wise error rate of $0.05$ (see Materials and methods), we identified 28 clusters of co-occurring TCRs, with sizes ranging from 7 to 386 TCRs (Figure 3). Given one of these clusters of co-occurring TCRs, we can count the number of cluster member TCRs found in each subject’s repertoire. The aggregate occurrence pattern of the cluster can be visualized as a rank plot of this cluster TCR count over the subjects (the black curves in Figure 3C–D). This ranking can also be compared with other phenotypic or genotypic features of the same subjects. In particular, by comparing this aggregate occurrence pattern to a control pattern generated by repeatedly choosing equal numbers of subjects independently at random (dotted green lines in Figure 3C–D), we can identify a subset of the cohort with an apparent enrichment of cluster member TCRs and look for overlap between this subset and other defined cohort features. Performing this comparison against the occurrence patterns of class I and class II HLA alleles revealed that the majority of the TCR clusters were strongly associated with at least one HLA allele (as depicted for a DRB1*15:01-associated cluster in Figure 3C and summarized in Figure 3B).

Figure 3 with 3 supplements see all

Download asset Open asset

Figure 3—source data 1

Cluster sizes and HLA-allele association p-values.

https://doi.org/10.7554/eLife.38358.010

Download elife-38358-fig3-data1-v2.tds

In addition, there were two large clusters of TCRs which were not strongly associated with any of the typed HLA alleles (clusters $6$ and $7$ in Figure 3). Visual inspection of the CDR3 regions of TCRs in one of these clusters revealed a distinctive ‘YV’ C-terminal motif that is characteristic of the TRBJ2-7*02 allele (Figure 3—figure supplement 1), and indeed the 41 subjects whose repertoires indicated the presence of this genetic variant were exactly the 41 subjects enriched for members of this TCR cluster (Figure 3D). This demonstrated that population diversity in germline allele sets manifests as occurrence pattern clustering. The other large, non-HLA associated TCR cluster had a number of distinctive features as well: strong preference for the TRBV06 family, followed by TRBV20 and TRBV04 (Figure 3—figure supplement 2); low numbers of inserted ‘N’ nucleotides; and a skewed age distribution biased toward younger subjects (Figure 3—figure supplement 3). These features, together with the lack of apparent HLA restriction, suggested that this cluster represented an invariant T cell subset, specifically MAIT (mucosal-associated invariant T) cells (Kjer-Nielsen et al., 2012; Venturi et al., 2013; Pogorelyy et al., 2017). Since MAIT cells are defined primarily by their alpha chain sequences, we searched in a recently published paired dataset (Howie et al., 2015) for partner chains of the clustered TCR$\beta$ chain sequences, and found a striking number that matched the MAIT consensus (TRAV1-2 paired with TRAJ20/TRAJ33 and a 12 residue CDR3, Figure 3—figure supplement 3D). We also looked for these clustered TCRs in a recently published MAIT cell sequence dataset (Howson et al., 2018) and found that 93 of the 138 cluster member TCRs occurred among the 31,654 unique TCRs from this dataset; of these 93 TCR$\beta$ chains, 27 were found among the 78 most commonly occurring TCRs in the dataset (the TCRs occurring in at least 7 of the 24 sequenced repertoires), a highly significant overlap ($P<2\times {10}^{-52}$ in a one-sided hypergeometric test). These concordances indicate that our untargeted approach has detected a well-studied T cell subset de novo through analysis of occurrence patterns.

These analyses suggested to us that TCR co-occurrence patterns across the full cohort of subjects are strongly influenced by the distribution of the HLA alleles, in accordance with the expectation that the majority of $\alpha \beta$ TCRs are HLA-restricted. Covariation between TCRs responding to the same HLA-restricted epitopes would only be expected in subjects positive for the restricting alleles, with TCR presence and absence outside these subjects likely introducing noise into the co-occurrence analysis. We therefore decided to analyze patterns of TCR co-occurrence within subsets of the cohort positive for specific HLA alleles, and to restrict our co-occurrence analysis to TCRs having a statistically significant association with the specific allele defining the cohort subset. To begin, we performed a comprehensive analysis of TCR-HLA association.

At a false discovery rate of 0.05 (estimated from shuffling experiments; see Materials and methods), we were able to assign 16,951 TCR$\beta$ sequences to an HLA allele (or alleles: DQ and DP alleles were analyzed as $\alpha \beta$ pairs, and there were 5 DR/DQ haplotypes whose component alleles were so highly correlated across our cohort that we could not assign TCR associations to individual DR or DQ components; see Materials and methods). Table 1 lists the top 50 HLA-associated TCR sequences by association p-value and top 10 associated TCRs for the well-studied A*02:01 allele.

We find that 8 of the top 10 A*02:01-associated TCRs have been previously reported and annotated as being responsive to viral epitopes, specifically influenza M1${}_{58}$ and Epstein-Barr virus (EBV) BMLF1${}_{280}$ (Shugay et al., 2018; Tickotsky et al., 2017). Moreover, each of these 8 TCR$\beta$ chains is present in a recent experimental dataset (Dash et al., 2017) that included tetramer-sorted TCRs positive for these two epitopes; each TCR has a clear similarity to one of the consensus epitope-specific repertoire clusters identified in that work, with the EBV TRBV20, TRBV29, and TRBV14 TCRs, respectively, matching the three largest branches of the BMLF1${}_{280}$ TCR tree, and the three influenza M1${}_{58}$ TCRs all matching the dominant TRBV19 ‘RS’ motif consensus (Figure 4—figure supplement 2). TCRs with annotation matches are sparser in the top 50 across all other alleles, which is likely due in part to a paucity of experimentally characterized non-A*02 TCRs, however we again see EBV-epitope responsive TCRs (with B*08:01 and B*35:01 restriction).

A global comparison of TCR feature distributions for HLA-associated versus non-HLA-associated TCRs provides further evidence of functional selection. As shown in Figure 4A, HLA-associated TCRs are on average more clonally expanded than a set of background, non-HLA associated TCRs with matching frequencies in the cohort. They also have lower generation probabilities—are harder to make under a simple random model of the VDJ rearrangement process—which suggests that their observed cohort frequencies may be elevated by selection (Figure 4B, see Materials and methods for further details on the calculation of clonal expansion indices and generation probabilities; also see (Pogorelyy et al., 2018)). Examination of two-dimensional feature distributions suggests that these shifts are correlated, with HLA-associated TCRs showing an excess of lower-probability, clonally expanded TCRs (Figure 4C); this trend appears stronger for class-I associated TCRs than for class II-associated TCRs (Figure 4—figure supplement 1).

Figure 4 with 2 supplements see all

Download asset Open asset

Figure 4—source data 1

Generation probabilities, clonal expansion indices, and allele associations for the TCRs analyzed here.

https://doi.org/10.7554/eLife.38358.015

Download elife-38358-fig4-data1-v2.tds

To give a global picture of TCR-HLA association, we counted the number of significant TCR associations found for each HLA allele in the dataset, and plotted this number against the number of subjects in the cohort with that allele (Figure 5). As expected, the more common HLA alleles have on average greater numbers of associated TCRs (since greater numbers of subjects permit the identification of more public TCRs, and the statistical significance assigned to an observed association of fixed strength grows as the number of subjects increases). What was somewhat more surprising is that the slope of the correlation between cohort frequency and number of associated TCRs varied dramatically among the HLA loci, with HLA-DRB1 alleles having the largest number of associated TCRs for a given allele frequency and HLA-C alleles having the smallest. The best-fit slope for the five DR/DQ haplotypes (12.2) was roughly the sum of the DR (7.99) and DQ (3.39) slopes, suggesting as expected that these haplotypes were capturing TCRs associated with both the DR and DQ component alleles. The smaller rate of TCR association observed at the HLA-C locus could be explained by a relatively lower level of cell surface expression of HLA-C alleles as well as their greater tendency to interact with killer cell immunoglobulin-like receptors (KIR) on natural killer (NK) cells (Kaur et al., 2017).

Figure 5 with 1 supplement see all

Download asset Open asset

Figure 5—source data 1

Allele frequencies and numbers of associated TCRs.

https://doi.org/10.7554/eLife.38358.018

Download elife-38358-fig5-data1-v2.tds

We assessed the accuracy of our TCR:HLA associations in two ways. First, we compared our HLA allele assignments to those given in the VDJdb database (which provides the peptide:MHC target and hence a putative HLA restriction for all entries; [Shugay et al., 2018]) and found that 90% of the VDJdb assignments for TCR$\beta$ chains present in both sets matched our associations. This agreement increases to 96% after filtering for the highest level of supporting evidence (VDJdb score of 3). Interestingly, two of the mismatches with VDJdb score three were from the protein structural database: the allo-complex between the B*08-restricted LC13 TCR and HLA-B*44:05 (Macdonald et al., 2009), and the structure of the A*02-restricted JM22 TCR bridged to a class II allele by a staphylococcal superantigen (Saline et al., 2010). In both of these cases, our data predict the canonical association: B*08 for the LC13 TCR$\beta$ chain and A*02 for the JM22 TCR$\beta$ chain. Second, we looked for HLA-associated public TCR$\beta$ chains in sequenced repertoires from T cell populations that were sorted for the presence of CD4/CD8 surface markers. One would expect that TCR$\beta$ chains associated with class I MHC molecules should be preferentially found in CD8+ populations, while class II-associated TCRs should be found in CD4+ populations. We selected four repertoire datasets (Emerson et al., 2013; Rubelt et al., 2016; Li et al., 2016; Oakes et al., 2017) with matched CD4+ and CD8+ repertoires from a total of $63$ individuals, and we analyzed the occurrence patterns of our HLA-associated TCR$\beta$ chains in these sequence datasets, producing for each TCR$\beta$ counts of the number of CD4+ and CD8+ repertoires it was observed in. Figure 5—figure supplement 1 shows that if we assign each TCR$\beta$ to the class (CD4+ or CD8+) with the higher count, these assignments are largely concordant with the MHC class of its associated HLA allele, and moreover this agreement increases as we increase either the stringency of HLA association or the stringency of the CD4/CD8 assignment (i.e., the minimum absolute difference between the CD4 and CD8 repertoire counts; see Materials and methods).

Having identified a set of HLA-associated TCR$\beta$ chains, we next sought to identify TCR clusters that might represent HLA-restricted responses to shared immune exposures. We performed this analysis for each HLA allele individually, restricting our clustering to the set of TCR chains significantly-associated with that allele and comparing occurrence patterns only over the subset of subjects positive for that allele. To reduce spurious co-occurrence signals driven by the presence/absence of other HLA alleles, we excluded TCR chains that were more strongly associated with a different HLA allele (i.e., not the one defining the cohort subset). The smaller size of many of these allele-positive cohort subsets reduces our statistical power to detect significant clusters using co-occurrence information. To counter this effect, we used the TCRdist similarity measure (Dash et al., 2017) to leverage the TCR sequence similarity which is often present within epitope-specific responses (Dash et al., 2017; Glanville et al., 2017) (see for example the A*02:01 TCRs in Table 1 and Figure 4—figure supplement 2). We augmented the probabilistic similarity measure used to define neighbors for DBSCAN clustering to incorporate information about TCR sequence similarity (as measured by TCRdist), in addition to cohort co-occurrence (see Materials and methods). We independently clustered each allele’s associated TCRs and merged the clustering results from all alleles; using the Holm multiple testing criterion (Holm, 1979) to limit the approximate family-wise error rate to 0.05, we found a total of $78$ significant TCR clusters.

We analyzed the sequences and occurrence patterns of the TCRs belonging to these 78 clusters in order to assess their potential biological significance and prioritize them for further study (Table 3). Each cluster was assigned two scores (Figure 6): a size score ($S_{\text{size}}$, $x$-axis), reflecting the significance of seeing a cluster of that size given the total number of TCRs clustered for its associated allele, and a co-occurrence score ($Z_{\text{CO}}$, $y$-axis), reflecting the degree to which the TCRs in that cluster co-occur within its allele-positive cohort subset (see Materials and methods). In computing the co-occurrence score, we defined a subset of individuals with an apparent enrichment for the member TCRs in each cluster; the size of this enriched subset of subjects is given in the ‘Subjects’ column in Table 3. We rank ordered the 78 clusters based on the sum of their size and co-occurrence scores (weighted to equalize dynamic range); the top five clusters are presented in greater detail in Figure 7 and Figure 8. HLA associations, member TCR and enriched subject counts, cluster center TCR sequences, scores, and annotations for all 78 clusters are given in Table 3.

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 6—source data 1

Paired TCR$\alpha$ chain sequences from the pairSEQ dataset of (Howie et al., 2015) for all clusters with at least $2$ matched TCR$\beta$ chains, along with a score for each cluster that assesses the degree of sequence similarity among the partner chains.

https://doi.org/10.7554/eLife.38358.021

Download elife-38358-fig6-data1-v2.tds

Figure 7

Download asset Open asset

Figure 8

Download asset Open asset

We found that a surprising number of the most significant HLA-restricted clusters had links to common viral pathogens. For example, the top cluster by both size and co-occurrence (Figure 7, upper panels) is an A*24:02-associated group of highly similar TCR$\beta$ chains, five of which can be found in a set of 12 TCR$\beta$ sequences reported to respond to the parvovirus B19 epitope FYTPLADQF as part of a highly focused CD8+ response to acute B19 infection (Kasprowicz et al., 2006). The subject TCR-counts curve for this cluster (Figure 7, top right panel) shows a strong enrichment of member TCRs in roughly 30% of the A*24:02 repertoires, which is on the low end of prevalence estimates for this pathogen (Heegaard and Brown, 2002) and may suggest that, if cluster enrichment does correlate with B19 exposure, there are likely to be other genetic or epidemiologic factors that determine which B19-exposed individuals show enrichment. The second most significant cluster by both measures is an A*02:01-associated group of TRBV19 TCRs with a high frequency of matches to the influenza M1${}_{58}$ response (41/43 TCRs, labeled ‘INF-pGIL’ for the first three letters of the GILGFVFTL epitope). Notably, the cluster member sequences recapitulate many of the core features of the tree of experimentally identified M1${}_{58}$ TCRs (Figure 4—figure supplement 2): a dominant group of length 13 CDR3 sequences with an ‘RS’ sequence motif together with a smaller group of length 12 CDR3s with the consensus CASSIG.YGYTF.

Rounding out the top five, the third and fifth most significant clusters also appear to be pathogen-associated. Cluster #3 brings together a diverse set of DRB1*07:01-associated TCR$\beta$ chains (Figure 8, top dendrogram), none of which matched our annotation database. However, it was strongly associated with CMV serostatus: As is evident in the subject TCR-counts panel for this cluster (Figure 8, top right), there is a highly significant ($P<3\times {10}^{-19}$) association between CMV seropositivity (blue dots at the bottom of the panel) and cluster enrichment (here defined as a subject TCR count $\ge 3$). Finally, the B*08:01-associated cluster #5 (bottom panels in Figure 8) appears to be EBV-associated: four of the TCR$\beta$ chains in this cluster match TCRs annotated as binding to EBV epitopes (two matches for the B*08:01-restricted FLRGRAYGL epitope and two for the B*08:01-restricted RAKFKQLL epitope). The fact that this cluster brings together sequence-dissimilar TCRs that recognize different epitopes from the same pathogen supports the hypothesis that at least some of the observed co-occurrence may be driven by a shared exposure.

As a preliminary validation of the clusters identified here, we examined the occurrence patterns of cluster member TCRs in two independent cohorts: a set of 120 individuals (‘Keck120’) that formed the validation cohort for the original Emerson et al. study, and a set of 86 individuals (‘Brit86’) taken from the aging study of (Britanova et al., 2016). Whereas the Keck120 repertoires were generated using the same platform as our 666-member discovery cohort, the Brit86 repertoires were sequenced from cDNA libraries using 5’-template switching and unique molecular identifiers. In the absence of HLA typing information for these subjects, we simply evaluated the degree to which each cluster’s member TCRs co-occurred over the entirety of each of these validation cohorts, using the co-occurrence score described above ($Z_{\text{CO}}^{\text{Keck}120}$ and $Z_{\text{CO}}^{\text{Brit}86}$ columns in Table 3). Although rare alleles and cluster-associated exposures may not occur with sufficient frequency in these smaller cohorts to generate co-occurrence signal, co-occurrence scores support the validity of the clusterings identified on the discovery cohort: 94% of the Keck120 scores and 92% of the Brit86 scores are greater than 0, indicating a tendency of the clustered TCRs to co-occur (smoothed score distributions are shown in Figure 6—figure supplement 1).

Given our large dataset of HLA-associated TCR$\beta$ sequences, we set out to look for correlations between CDR3 sequence and HLA allele sequence. Previous studies have identified correlations between TCR V-gene usage and HLA alleles (Sharon et al., 2016; Blevins et al., 2016); these correlations are consistent with a picture of TCR:peptide:MHC interactions in which the CDR1 and CDR2 loops (whose sequence is determined by the V gene) primarily contact the MHC while the CDR3 loops contact the peptide. To complement these studies and leverage our large set of HLA-associated sequences, we set out to look for correlations between the CDR3 sequence itself and the HLA allele. In our previous work on epitope-specific TCRs (Dash et al., 2017), we identified a significant negative correlation between CDR3 charge and peptide charge, suggesting a tendency toward preserving charge complementarity across the TCR:pMHC interface. Although the CDR3 loop primarily contacts the MHC-bound peptide, computational analysis of solved TCR:peptide:MHC structures in the Protein Data Bank (Berman et al., 2000) (see Materials and methods) identified a number of HLA sequence positions that are frequently contacted by CDR3 amino acids (Table 2). For each frequently-contacted HLA position with charge variability among alleles we computed the covariation between HLA allele charge at that position and average CDR3 charge for allele-associated TCRs. Since portions of the CDR3 sequence are contributed by the V- and J-gene germline sequences, and covariations are known to exist between HLA and V-gene usage, we also performed a covariation analysis restricting to ‘non-germline’ CDR3 sequence positions whose coding sequence is determined by at least one non-templated insertion base (based on the most parsimonious VDJ reconstruction; see Materials and methods). We found a significant negative correlation ($R=-0.47,P<4\times {10}^{-4}$ for the full CDR3 sequence; $R=-0.52,P<7\times {10}^{-5}$ for the non-germline CDR3 sequence) between CDR3 charge and the charge at position 70 of the class II beta chain (correcting these p-values for the fact that we considered $7$ positions yields $2.3\times {10}^{-3}$ and $4.3\times {10}^{-4}$). We did not see a significant correlation for the frequently contacted position on the class II alpha chain, perhaps due to the lack of sequence variation at the DR$\alpha$ locus and/or the more limited number of DQ$\alpha$ and DP$\alpha$ alleles. None of the five class I positions showed significant correlations, which could be due to their lower contact frequencies, a smaller average number of associated TCRs ($51$ for class I versus $309$ for class II), bias toward A*02 in the structural database, or noise introduced from multiple contacted positions varying simultaneously. Further analysis of the class II correlation suggested that it was driven largely by HLA-DRB1 alleles: position 70 correlations were $-0.56$ versus $-0.10$ for DR and DQ, respectively, over the full CDR3 and $-0.64$ vs $-0.38$ for the non-germline CDR3. Figure 9 provides further detail on this DRB1-TCR charge anti-correlation, including a structural superposition showing the proximity of position 70 to the TCR$\beta$ CDR3 loop.

Figure 9

Download asset Open asset

Figure 9—source data 1

Charge at position 70 and average CDR3 charge of allele-associated TCRs for $30$ HLA-DRB1 alleles.

https://doi.org/10.7554/eLife.38358.026

Download elife-38358-fig9-data1-v2.tds

We analyzed the HLA associations of strongly CMV-associated TCR$\beta$ chains to gain insight into their predictive power across genetically diverse individuals. Here we change perspective somewhat from earlier sections, in that we select TCRs based on their CMV association and then evaluate HLA association, rather than the other way around. In their original study, Emerson et al. identified a set of TCR$\beta$ chains that were enriched in CMV seropositive individuals and showed that by counting these CMV-associated TCR$\beta$ chains in a query repertoire they could successfully predict CMV serostatus both in cross-validation and on an independent test cohort. The success of this prediction strategy across a diverse cohort of individuals raises the intriguing question of whether these TCR$\beta$s are primarily HLA-restricted in their occurrence and in their association with CMV, or whether they span multiple HLA types. To shed light on this question we focused on a set of $68$ CMV-associated TCR$\beta$ chains whose co-occurrence with CMV seropositivity was significant at a p-value threshold of 1.5e-5 (corresponding to an FDR of $0.05$; see Materials and methods). For each CMV-associated TCR$\beta$ chain, we identified its most strongly associated HLA allele and compared the p-value of this association to the p-value of its association with CMV (Figure 10A). From this plot we can see that the majority of the CMV-associated chains do appear to be HLA-associated, having p-values that exceed the FDR $0.05$ threshold for HLA association. The excess of highly significant HLA-association p-values for these CMV-associated TCR$\beta$s can be seen in Figure 10B, which compares the observed p-value distribution to a background distribution of HLA association p-values for randomly selected frequency-matched public TCR$\beta$s.

Figure 10

Download asset Open asset

Figure 10—source data 1

Full and subsetted CMV- and HLA- association p-values for $68$ CMV-associated TCRs.

https://doi.org/10.7554/eLife.38358.028

Download elife-38358-fig10-data1-v2.tds

As a next step we looked to see whether these HLA associations fully explained the CMV association, in the sense that the CMV association was only present in subjects positive for the associated allele. For each of the $68$ CMV-associated TCRs, we divided the cohort into subjects positive for its most strongly associated HLA allele and subjects negative for that allele. Here we considered both 2- and 4-digit resolution alleles when defining the most strongly associated allele, to allow for TCRs whose association extends beyond a single 4-digit allele. We computed association p-values between TCR occurrence and CMV seropositivity over these two cohort subsets independently and compared them (Figure 10C). We see that the majority of the points lie below the $y=x$ line—indicating a stronger CMV-association on the subset of the cohort positive for the associated allele—and also below the line corresponding to the expected minimum of $68$ uniform random variables (i.e. the expected upper significance limit in the absence of CMV association on the allele-negative cohort subsets). There are however a few TCR$\beta$s which do not appear strongly HLA-associated and for which the CMV-association remains strong in the absence of their associated allele (the points above the line $y=x$ in Figure 10C). For example, the public TCR$\beta$ chain defined by TRBV07 and the CDR3 sequence CASSSDSGGTDTQYF (which corresponds to the highest point in Figure 10C) is strongly CMV-associated ($22/23$ subjects with this chain are CMV positive; $P<3\times {10}^{-7}$) but does not show evidence of HLA association in our dataset. TCRs with HLA promiscuity may be especially interesting from a diagnostic perspective, since their phenotype associations may be more robust to differences in genetic background.

Finally, we looked to see whether CMV association completely explained the observed HLA associations, in the sense that a response to one or more CMV epitopes was likely the only driver of HLA association, or whether there might be evidence for other epitope-specific responses by these TCR$\beta$ chains or a more general affinity for the associated allele, perhaps driven by common self antigens. Put another way, do we see evidence for pre-existing enrichment of any of these TCR$\beta$ chains when their associated allele is present, even in the absence of CMV, which might suggest that the CMV response recruits from a pre-selected pool enriched for TCRs with intrinsic affinity for the restricting allele? To approach this question we split the cohort into CMV seropositive and seronegative subjects and computed, for each of the $68$ CMV-associated TCRs, the strength of its association with its preferred allele over these two subsets separately. Figure 10D compares these HLA-association p-values computed over the subsets of the cohort positive (289 individuals, $x$-axis) and negative (352 individuals, $y$-axis) for CMV. We can see in this case that all of the associations on the CMV-positive subset are stronger than those on the CMV-negative subset, and indeed the CMV-negative p-values do not appear to exceed random expectation given the number of comparisons performed. Thus, the apparent lack of any significant HLA-association on the CMV-negative cohort subset suggests that the HLA associations of these CMV-predictive chains are largely driven by CMV exposure. A limitation of this analysis is that, although the CMV-negative subset of the cohort is larger than the CMV-positive subset, the number of TCR occurrences in the CMV-negative subset is likely lower than in the CMV-positive subset for these CMV-associated chains, which will limit the strength of the HLA associations that can be detected.

TCR$\beta$ repertoire sequence data for the 666 members of the discovery cohort was downloaded from the Adaptive biotechnologies website using the link provided in the original (Emerson et al., 2017) publication (https://clients.adaptivebiotech.com/pub/Emerson-2017-NatGen). The repertoire sequence data for the 120 individuals in the ‘Keck120’ validation set was included in the same download. Repertoire sequence data for the 86 individuals in the ‘Brit86’ validation set was downloaded from the NCBI SRA archive using the Bioproject accession PRJNA316572 (Britanova et al., 2016) and processed using scripts and data supplied by the authors (https://github.com/mikessh/aging-study) in order to demultiplex the samples and remove technical replicates. Repertoire sequence data for TCR$\beta$ chains from MAIT cells was downloaded from the NCBI SRA archive using the Bioproject accession PRJNA412739 (Howson et al., 2018). Repertoire sequence data for TCR$\beta$ chains from T cells sorted for CD4/CD8 surface markers were taken from the following studies: (Emerson et al., 2013), available for download at https://clients.adaptivebiotech.com/pub/emerson-2013-jim; (Rubelt et al., 2016), downloaded from the NCBI SRA archive using the Bioproject accession PRJNA300878; (Li et al., 2016), downloaded from the NCBI SRA archive using the Bioproject accession PRJNA348095; and (Oakes et al., 2017), downloaded from the NCBI SRA archive using the Bioproject accession PRJNA390125.

V and J genes were assigned by comparing the TCR nucleotide sequences to the IMGT/GENE-DB (Giudicelli et al., 2005) nucleotide sequences of the human TR genes (sequence data downloaded on 9/6/2017 from http://www.imgt.org/genedb/). CDR3 nucleotide and amino acid sequences and most-parsimonious VDJ recombination scenarios were assigned by the TCRdist pipeline (Dash et al., 2017) (the most parsimonious recombination scenario, used for identifying non-germline CDR3 amino acids, is the one requiring the fewest non-templated nucleotide insertions). To define the occurrence matrix of public TCRs and assess TCR-TCR, TCR-HLA and TCR-CMV association, a TCR$\beta$ chain was identified by its CDR3 amino acid sequence and its V-gene family (e.g., TRBV6-4*01 was reduced to TRBV06). TCR sequence reads for which a unique V-gene family could not be determined (due to equally well-matched V genes from different families, a rare occurrence in this dataset) were excluded from the analysis. The matrix $M$ of public TCR$\beta$ occurrences across the discovery cohort, HLA allele occurrence patterns, and other associated data needed to reproduce the findings of this study have been deposited in the Zenodo database (doi:10.5281/zenodo.1248193).

A preliminary analysis of TCR sharing at the nucleotide level was conducted to identify potential cross-contamination in the discovery cohort repertoires. Each TCR$\beta$ nucleotide sequence that was found in multiple repertoires was assigned a generation probability ($P_{\text{gen}}$, see below) in order to identify nucleotide sequences with suspiciously high sharing rates among repertoires. Visual comparison of the sharing rate (the number of repertoires in which each TCR$\beta$ nucleotide sequence was found) to the generation probability (Figure 11) showed that the majority of highly-shared TCRs had correspondingly high generation probabilities; it also revealed a cluster of TCR chains with unexpectedly high sharing rates. Examination of the sequences of these highly-shared TCRs revealed them to be variants of the consensus sequence CFFKQKTAYEQYF (coding sequence: tgttttttcaagcagaagacggcatacgagcagtacttc). Consultation with scientists at Adaptive Biotechnologies confirmed that these sequences were likely to represent a technical artifact of the sequencing pipeline. We elected to remove all TCR$\beta$ nucleotide sequences whose sharing rates put them outside the decision boundary indicated by the black line in Figure 11, which eliminated the vast majority of the artifactual variants as well as a handful of other highly shared, low-probability sequences ($592$ nucleotide sequences in total were removed).

Figure 11

Download asset Open asset

Figure 11—source data 1

TCR$\beta$ nucleotide sequences excluded from our analysis.

https://doi.org/10.7554/eLife.38358.030

Download elife-38358-fig11-data1-v2.tds

Each public TCR$\beta$ chain was assigned a clonal expansion index ($I_{\text{exp}}$) determined by its frequencies in the repertoires in which it was found. First, the unique TCR$\beta$ chains present in each repertoire were ordered based on their inferred nucleic acid template count (Carlson et al., 2013), and assigned a rank ranging from $0$ (lowest template count) to $S-1$ (highest template count), where $S$ is the total number of chains present in the repertoire. TCRs with the same template count were assigned the same tied rank equal to the midpoint of the tied group. In order to compare across repertoires, the ranks for each repertoire were then normalized by dividing by the number of unique sequences in the repertoire. The clonal expansion index for a given public TCR $t$ was taken to be its average normalized rank for the repertoires in which it occurred:

where the sum is taken over the $m$ repertoires in which $t$ is found, $r_i$ is the template-count rank of TCR $t$ in repertoire $i$, and $S_i$ is the total size of repertoire $i$.

HLA genotyping was performed and confirmed by molecular means, including sequence specific oligonucleotide probe typing (SSOP), Sanger sequencing (SBT) or next generation sequencing (NGS) (Smith et al., 2014). Independently, HLA alleles were imputed using data generated by high density single-nucleotide polymorphism arrays as previously described (Martin et al., 2017). Imputed alleles were compared with HLA typing data from SBT and NGS, and used to resolve ambiguous HLA codes generated by SSOP and provide a uniform set of four digit allele assignments. HLA typing data availability varied across loci as follows: HLA-A ($629$ subjects), HLA-B ($630$ subjects), HLA-C ($629$ subjects), HLA-DRB1 ($630$ subjects), HLA-DQA1 ($522$ subjects), HLA-DQB1 ($630$ subjects), HLA-DPA1 ($606$ subjects), and HLA-DPB1 ($472$ subjects). When calculating the association p-values between TCR$\beta$ chains and HLA alleles reported in Table 1, the cohort was restricted to the subset of subjects with available HLA typing at the relevant locus. For comparing TCR association rates across loci in Figure 5, associations were calculated over the cohort subset ($522$ subjects) with typing data at all compared loci (A, B, C, DRB1, DQA1, and DQB1) in order to avoid spurious differences in association strengths arising from differential data availability among the loci. Due to their very strong linkage on our cohort, five DR-DQ haplotypes were treated as single allele units for association calculations and clustering: DRB1*03:01-DQA1*05:01-DQB1*02:01, DRB1*15:01-DQA1*01:02-DQB1*06:02, DRB1*13:01-DQA1*01:03-DQB1*06:03, DRB1*10:01-DQA1*01:05-DQB1*05:01, and DRB1*09:01-DQA1*03:02-DQB1*03:03.

We implemented a version of the probabilistic model proposed by Walczak and co-workers (Murugan et al., 2012) in order to assign to each public TCR$\beta$ chain (defined by a V-gene family and a CDR3 amino acid sequence) a generation probability, $P_{\text{gen}}$, which captures the probability of seeing that TCR$\beta$ in the preselection repertoire. $P_{\text{gen}}$ is calculated by summing the probabilities of the possible VDJ rearrangements that could have produced the observed TCR:

where $S$ represents the set of possible VDJ recombination scenarios capable of producing the observed TCR V family and CDR3 amino acid sequence. To compute the probability of a given recombination scenario $s$, we use the factorization proposed by Marcou et al. (2018), which captures observed dependencies of V-, D-, and J-gene trimming on the identity of the trimmed gene and of inserted nucleotide identity on the identity of the preceding nucleotide:

Here the recombination scenario $s$ consists of a choice of V gene ($V_s$), D gene ($D_s$), J gene ($J_s$), number of nucleotides trimmed back from the end of the V gene (${{\displaystyle \text{del}}}_{s}V$) or J gene (${{\displaystyle \text{del}}}_{s}J$) or D gene (${{\displaystyle \text{del}}}_{s}D{5}^{\prime }$ and ${{\displaystyle \text{del}}}_{s}D{3}^{\prime }$), number of nucleotides inserted between the V and D genes (${{\displaystyle \text{Ins}}}_{s}VD$) and between the D and J genes (${{\displaystyle \text{Ins}}}_{s}DJ$) and the identities of the inserted nucleotides ($\left\{n_i\right\}$ and $\left\{m_i\right\}$ respectively). At the start of the calculation, the CDR3 amino acid sequence is converted to a list of potential degenerate coding nucleotide sequences (here degenerate means that nucleotide class symbols such as W (for A and T) and R (for A and G) are allowed). Since each amino acid other than Leucine, Serine, and Arginine has a single degenerate codon (P=CCN, N = AAY, K = AAR, etc.) and these three amino acids have two such codons (S={TCN,AGY}, R={CGN,AGR}, L={CTN,TTR}), this list of nucleotide coding sequences is generally not too long. The generation probability is then taken to be the sum of the probabilities of these degenerate nucleotide sequences. Since the total number of possible recombination scenarios is in principle quite large, we make a number of approximations to speed the calculation: we limit excess trimming of genes to at most three nucleotides, where excess trimming is defined to be trimming back a germline gene nucleotide which matches the target CDR3 nucleotide (therefore requiring non-templated reinsertion of the same nucleotide); at most two palindromic nucleotides are allowed; sub-optimal D gene alignments are only considered up to a score gap of 2 matched nucleotides relative to the best match. The parameters of the probability model are fit by a simple iterative procedure in which we generate rearranged sequences using an initial model, compare the statistics of those sequences to statistics derived from observed out-of-frame rearrangements in the dataset, and adjust the probability model parameters to iteratively improve agreement. We compared the nucleotide sequence generation probabilities computed using our software with those computed using the published tool IGoR (Marcou et al., 2018) and found good overall agreement: a linear regression analysis of the ${\displaystyle {\log }_{10}\left(P_{\text{gen}}\right)}$ values gives a correlation coefficient $R=0.97$ with slope of $0.98$ and an intercept of $0.22$ for a set of $800$ randomly selected TCR$\beta$ chains.

We performed an analysis of covariation across the cohort for pairs of TCR chains and for TCR chains and HLA alleles (Figure 12). We used the hypergeometric distribution to assess the significance of an observed overlap between two subsets of the cohort (for example, the subset of subjects positive for a given HLA allele and the subset of subjects with a given TCR$\beta$ chain in their repertoires), taking our significance p-value to be the probability of seeing an equal or greater overlap if the two subsets had been chosen at random:

where $k$ is the size of the overlap, $N_1$ and $N_2$ are the sizes of the two subsets, and $N$ is the total cohort size (i.e., the number of individuals in the cohort). We use $P_{\text{overlap}}$ to assess the significance of an overlap $C_a\cap C_t$ between an HLA allele $a$ found in the cohort subset $C_a$ and a TCR$\beta$ chain $t$ found in the cohort subset $C_t$ as follows:

where $\left|C\right|$ denotes the cardinality of the set $C$. A complication arises when assessing TCR-TCR co-occurrence in the presence of variable-sized repertoires: TCRs are more likely to come from the larger repertoires than the smaller ones, which violates the assumptions of the hypergeometric distribution and leads to inflated significance scores. In particular, when we use the hypergeometric distribution to model the overlap between the sets of subjects in which two TCR chains are found, we implicitly assume that all subjects are equally likely to belong to a TCR chain’s subject set. If the subject repertoires vary in size, this assumption will not hold. For example, in the limit of a subject with an empty repertoire, no TCR subject sets will contain that subject, which will inflate all the overlap p-values since we are effectively overstating the size $N$ of the cohort by $1$. On the other hand, if one of the subject repertoires contains all the public TCR chains, then each TCR-TCR overlap will automatically contain that subject, again inflating the p-values since we are artificially adding $1$ to each of $k$, $N_1$, $N_2$, and $N$. We developed a simple heuristic to correct for this effect using a per-subject bias factor by defining

where $S_i$ is the size of repertoire $i$ and $N$ is the cohort size. To score an overlap between the occurrence patterns of two TCR$\beta$ chains $t$ and $t^{\prime }$, where $t$ is found in the subset $C_t$ of the cohort, $t^{\prime }$ is found in the subset $C_{t^{\prime }}$, and their overlap $C_t\cap C_{t^{\prime }}$ contains the $k$ subjects $s_1,\mathrm{...},s_k$, we adjust the overlap p-value ($P_{\text{overlap}}$) by the product of the bias factors of the subjects in the overlap:

Here we are multiplying the hypergeometric p-value ($P_{\text{overlap}}$) by a term that corrects for the fact that not all overlaps of size $k$ are equally likely (the product of the $k$ bias factors captures the relative bias toward the observed overlap). This has the effect of decreasing the significance assigned to overlaps involving larger repertoires, yet remains fast to evaluate, an important consideration given that the all-vs-all TCR co-occurrence calculation involves about ${10}^{14}$ pairwise comparisons (and this calculation is repeated multiple times with shuffled occurrence patterns to estimate false-discovery rates). When clustering by co-occurrence, we augmented this heuristic p-value correction by also eliminating repertoires with very low (fewer than 30,000) or very high (more than 120,000) numbers of public TCR$\beta$ chains (nonzero entries in the occurrence matrix $M$), as well as five additional repertoires which showed anomalously high levels of TCR nucleotide sharing with another repertoire—all with the goal of reducing potential sources of spurious TCR-TCR co-occurrence signal.

Figure 12

Download asset Open asset

We used the approach of (Storey and Tibshirani, 2003) to estimate false-discovery rates for detecting associations between TCRs and HLA alleles and between TCRs and CMV seropositivity. Briefly, for a fixed significance threshold $P$ we estimate the false-discovery rate (FDR) by randomly permuting the HLA allele or CMV seropositivity assignments $20$ times and computing the average number of significant associations discovered at the threshold $P$ in these shuffled datasets. The estimated FDR is then the ratio of this average shuffled association number to the number of significant associations discovered in the true dataset at the same threshold. In order to estimate a false-discovery rate for TCR-TCR co-occurrence over the full cohort, we performed $20$ co-occurrence calculations on shuffled occurrence matrices, preserving the per-subject bias factors during shuffling by resampling each TCR’s occurrence pattern with the bias distribution $\left\{b_i\right\}$ determined by the subject repertoire sizes.

We assessed the accuracy of our TCR:HLA associations by looking for HLA-associated public TCR$\beta$ chains in sequenced repertoires from T cell populations that were sorted for the presence of CD4/CD8 surface markers. We selected four repertoire datasets with matched CD4+ and CD8+ repertoires from a total of $63$ individuals (see the section Datasets for access details; [Emerson et al., 2013; Rubelt et al., 2016; Li et al., 2016; Oakes et al., 2017]). We analyzed the occurrence patterns of HLA-associated TCR$\beta$ chains in these sequence datasets, producing for each TCR$\beta$ counts of the number of CD4+ and CD8+ repertoires it was observed in ($N_{\text{CD}4}$ and $N_{\text{CD}8}$). TCR$\beta$ abundance levels within the individual repertoires were ignored; each occurrence in a repertoire contributed a single count to the respective CD4 or CD8 total (which therefore range between $0$ and $63$). Given a threshold $\delta$ on the CD4/CD8 counts difference, we assign to the CD4 compartment all TCRs for which $N_{\text{CD}4}-N_{\text{CD}8}\ge \delta$, and we assign to the CD8 compartment all TCRs for which $N_{\text{CD}8}-N_{\text{CD}4}\ge \delta$. Figure 5—figure supplement 1 shows the concordance between these assignments and inferences based on the HLA class of the most strongly associated HLA allele, for all significantly associated TCR$\beta$ chains and for various threholds $\delta$.

We used the DBSCAN (Ester et al., 1996) algorithm to cluster public TCR$\beta$ chains by their occurrence patterns. DBSCAN is a simple and robust clustering procedure that requires two input parameters: a similarity/distance threshold ($T_{\text{sim}}$) at which two points in the dataset are considered to be neighbors, and a minimum number of neighbors ($N_{\text{core}}$) for a point to be considered a core, as opposed to a border, point. DBSCAN clusters consist of the connected components of the neighbor-graph over the core points, together with any border point neighbors the core cluster members have. To prevent the discovery of fictitious clusters, $T_{\text{sim}}$ and $N_{\text{core}}$ can be selected so that core points (points with at least $N_{\text{core}}$ neighbors) are unlikely to occur by chance. There is a trade-off between the two parameter settings: as $T_{\text{sim}}$ is relaxed, points will tend to have more neighbors on average and thus $N_{\text{core}}$ should be increased, which biases toward discovery of larger clusters; conversely, more stringent settings of $T_{\text{sim}}$ are compatible with smaller values for $N_{\text{core}}$ which permits the discovery of smaller, more tightly linked clusters.

For clustering TCRs by co-occurrence over the full cohort, we used a threshold of $T_{\text{sim}}$$={10}^{-8}$ and chose a value for $N_{\text{core}}$ ($6$) such that no core points were found in any of the $20$ shuffled datasets. In other words, two TCRs $t_1$ and $t_2$ were considered to be neighbors for DBSCAN clustering if $P_{\text{CO}}(t_1,t_2)<{10}^{-8}$; a TCR was considered a core point if it had at least $6$ neighbors. Choosing parameters for HLA-restricted TCR clustering was slightly more involved due to the variable number of clustered TCRs for different alleles, and the more complex nature of the similarity metric, whose dependence on TCR sequence makes shuffling-based approaches more challenging. To begin, we transformed the TCRdist sequence-similarity measure into a significance score $P_{\text{TCRdist}}$ which captures the probability of seeing an observed or smaller TCRdist score for two randomly selected TCR$\beta$ chains. Since public TCR$\beta$ chains are on average shorter and closer to germline than private TCRs, we derived the $P_{\text{TCRdist}}$ CDF by performing TCRdist calculations on randomly selected public TCRs seen in at least $5$ repertoires. We identified neighbors for DBSCAN clustering using a similarity score $P_{\text{sim}}$ that combines co-occurrence and TCR sequence similarity:

where the transformation by ${\displaystyle f\left(x\right)=x-x\log \left(x\right)}$ corrects for taking the product of two p-values because $f\left(x\right)$ is the cumulative distribution function of the product of two uniform random variables. Thus, if $P_{\text{TCRdist}}$ and $P_{\text{CO}}$ are independent and uniformly distributed, the same will be true of $P_{\text{sim}}$.

For HLA-restricted clustering using this combined similarity measure we set a fixed value of $T_{\text{sim}}$$={10}^{-4}$ and adjusted the $N_{\text{core}}$ parameter as a function of the total number of TCRs clustered for each allele. As in global clustering, our goal was to choose $N_{\text{core}}$ such that core points were unlikely to occur by chance (more precisely, had a per-allele probability less than $0.05$). We estimated the probability of seeing core points by modeling neighbor number using the binomial distribution, assuming that the observed neighbor number of a given TCR during clustering is determined by $M-1$ independent Bernoulli-distributed neighborness tests with rate $r$, where $M$ is the number of clustered TCRs. Rather than assuming a fixed neighbor-rate $r$ across TCRs, we captured the observed variability in neighbor-rate (due, for example, to unequal V-gene frequencies and variable CDR3 lengths) by using a mixture of 20 rates $\left\{r_j\right\}$ estimated from similarity comparisons on randomly chosen public TCRs. More precisely, we choose the smallest value of $N_{\text{core}}$ for which the following inequality holds (where $M$ is the number of clustered TCRs for the allele in question):

We also used this neighbor-number model to assign a p-value ($P_{\text{size}}$) to each cluster reflecting the likelihood of seeing the observed degree of clustering by chance. Since DBSCAN clusters are effectively single-linkage-style partitionings of the core points (together with any neighboring border points), they can have a variety of shapes, ranging from densely interconnected graphs, to extended clusters held together by local neighbor relationships (Ester et al., 1996). Modeling the total size of these arbitrary groupings is challenging, so we took the simpler and more conservative approach of assigning p-values based on the size of the largest TCR neighborhood (set of neighbors for a single TCR) contained within each cluster. We identified the member TCR with the greatest number of neighbors in each cluster (the cluster center) and computed the likelihood of seeing an equal or greater neighbor-number under the mixture model described above. This significance estimate is conservative in that it neglects clustering contributions from TCRs outside the neighborhood of the cluster center, however in practice we observed that the majority of TCR clusters were dominated by a single dense region of repertoire space and therefore reasonably well-captured by a single neighborhood. To control false discovery when combining DBSCAN clusters from independent clustering runs for different HLA alleles, we used the Holm method (Holm, 1979) applied to the sorted list of cluster $P_{\text{size}}$ values, with a target family-wise error rate (FWER) of $0.05$ (i.e., we attempted to limit the overall probability of seeing a false cluster to $0.05$). In the Holm FWER calculation we set the total number of hypotheses equal to the total number of TCRs clustered across all alleles minus the cumulative neighbor-count of the cluster centers (we exclude cluster center neighbors since their neighbor counts are not independent of the neighbor count of the cluster center). When performing HLA-restricted clustering, each TCR$\beta$ chain was assigned to its most strongly associated HLA allele. Where two alleles had identical or nearly identical (within a factor of $1.25$) association p-values, the TCR chain was included in the clustering analysis for both alleles.

For each (global or HLA-restricted) TCR cluster, we analyzed the occurrence patterns of the member TCRs in order to identify a subset of the (full or allele-positive) cohort enriched for those TCRs. We counted the number of cluster member TCRs found in each subject’s repertoire and sorted the subjects by this TCR count (rank plots in Figure 3B–C and in the right panels of Figure 7). For comparison, we generated control TCR count plots by independently resampling the subjects for each member TCR, preserving the frequency of each TCR and biasing by subject repertoire size. Each complete resampling of the cluster member TCR occurrence patterns produced a subject TCR rank plot; we repeated this resampling process $1000$ times and averaged the rank plots to yield the green (‘randomized’) curves in Figure 3B–C and Figure 7. To compare the observed and randomized curves, we took a signed difference

between the observed counts $C_j$ and the randomized counts $R_j$, where the value of the subject index $i=i_{\text{max}}$ that maximizes the right-hand side in the equation above represents a switchpoint below which the observed counts generally exceed the randomized counts and above which the reverse is true (both sets of counts are sorted in decreasing order). We take this switchpoint $i_{\text{max}}$ as an estimate of the number of enriched subjects for the given cluster (this is the value given in the ‘Subjects’ column in Table 3).

Since the raw $D_{\text{CO}}$ values are not comparable between clusters of different sizes and for different alleles, we transformed these values to a Z-score ($Z_{\text{CO}}$) by generating, for each cluster, $1000$ additional random TCR count curves and computing the mean (${\mu }_D$) and standard deviation (${\sigma }_D$) of their $D_{\text{CO}}^{\text{rand}}$ score distribution:

We used this co-occurrence score $Z_{\text{CO}}$ together with a log-transformed version of the cluster size p-value,

for visualizing clustering results in Figure 6 ($S_{\text{size}}$ on the $x$-axis and $Z_{\text{CO}}$ on the $y$-axis) and prioritizing individual clusters for detailed follow-up.

We annotated public TCRs in our dataset by matching their sequences against two publicly available datasets: VDJdb (Shugay et al., 2018), a curated database of TCR sequences with known antigen specificities (downloaded on 3/29/18; about $17,000$ human TCR$\beta$ entries) and McPAS-TCR (Tickotsky et al., 2017), a curated database of pathogen-associated TCR sequences (downloaded on 3/29/18; about $9,000$ human TCR$\beta$ entries). VDJdb entries are associated with a specific MHC-presented epitope, whereas McPAS-TCR also includes sequences of TCRs isolated from diseased tissues whose epitope specificity is not defined. We added to this merged annotation database the sequences of structurally characterized TCRs of known specificity (see below), as well as literature-derived TCRs from a handful of primary studies (Dash et al., 2017; Glanville et al., 2017; Song et al., 2017; Kasprowicz et al., 2006). For matches between HLA-associated TCRs and database TCRs of known specificity, we filtered for agreement (at 2-digit resolution) between the associated HLA allele in our dataset and the presenting allele from the database. In other words, TCRs belonging to B*08:01-restricted clusters were not annotated with matches to database TCRs that bind to A*02:01-presented peptides.

We analyzed a set of experimentally determined TCR:peptide-MHC structures to find MHC positions frequently contacted by the CDR3$\beta$ loop. Crystal structures of complexes involving human TCRs and human class I or class II HLA alleles (Table 4) were identified using BLAST (Altschul et al., 1997) searches against the RCSB PDB (Berman et al., 2000) sequence database (ftp://ftp.wwpdb.org/pub/pdb/derived_data/pdb_seqres.txt). Structural coverage of HLA loci and alleles is sparse and highly biased toward well studied alleles such as HLA-A*02. Given the high degree of structural similarity among class I and among class II MHC structures solved to date, we elected to share contact information across loci using trans-locus sequence alignments. For class I we used the merged alignment (ClassI_prot.txt) available from the IPD-IMGT/HLA (Robinson et al., 2015) database. Starting with multiple sequence alignments for individual class II loci from the IPD-IMGT/HLA database, we inserted gaps as needed in order to created merged alignments for the class II $\alpha$ and $\beta$ chains. These alignments provided a common reference frame in which to combine residue-residue contacts from the TCR:peptide-MHC structures. We considered two amino acid residues to be in contact if they had a side chain heavyatom contact distance less than or equal to $4.5$Å. The CDR3$\beta$ contact frequency for an alignment position (class I, class II-$\alpha$, or class II-$\beta$) was defined to be the total number of contacted CDR3$\beta$ amino acids observed for that position, divided by the total number of structures analyzed. Redundancy in the structural database was assessed at the level of TCR and HLA sequence, ignoring the sequence of the peptide. Contacts from a set of $n$ structures all containing the same TCR and HLA were given a weight of $1/n$ when computing the residue contact frequencies. The statistical significance of correlations between HLA allele charge and average HLA-associated TCR CDR3 charge were computed using a 2-sided test as implemented in the function scipy.stats.linregress.

C++ source code implementing the clustering, generation probability, and correlation algorithms described here is available at https://github.com/phbradley/pubtcrs (copy archived at https://github.com/elifesciences-publications/pubtcrs [Bradley, 2018]).

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This work was supported in part through the NIH/NCI Cancer Center Support Grant P30 CA015704 and by NIH NHLBI grant R01-HL105914 to JH, as well as R01 GM113246 and U19 AI117891. The research of Frederick Matsen was supported in part by a Faculty Scholar grant from the Howard Hughes Medical Institute and the Simons Foundation. We gratefully acknowledge superlative computing support from Fred Hutch scientific computing and thank Paul Thomas and Jeremy Crawford for helpful comments on a preliminary version of this manuscript.

Human subjects: All samples were collected and analyzed, and informed consent and consent to publish were obtained, according to research protocols approved by the Fred Hutchinson Cancer Research Center (FHCRC) Institutional Review Board.

1. Received: May 14, 2018
2. Accepted: August 21, 2018
3. Accepted Manuscript published: August 28, 2018 (version 1)
4. Version of Record published: September 28, 2018 (version 2)

© 2018, DeWitt et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.