Method for identification of condition-associated public antigen receptor sequences

Diverse repertoires of hypervariable immunoglobulin receptors (TCR and BCR) recognize antigens in the adaptive immune system. The development of immunoglobulin receptor repertoire sequencing methods makes it possible to perform repertoire-wide disease association studies of antigen receptor sequences. We developed a statistical framework for associating receptors to disease from only a small cohort of patients, with no need for a control cohort. Our method successfully identifies previously validated Cytomegalovirus and type 1 diabetes responsive receptors.

Diverse repertoires of hypervariable immunoglobulin receptors (TCR and BCR) recognize antigens in the adaptive immune system. The development of immunoglobulin receptor repertoire sequencing methods makes it possible to perform repertoire-wide disease association studies of antigen receptor sequences. We developed a statistical framework for associating receptors to disease from only a small cohort of patients, with no need for a control cohort. Our method successfully identifies previously validated Cytomegalovirus and type 1 diabetes responsive receptors.
T-cell receptors (TCR) and B-cell receptors (BCR) are hypervariable immunoglobulins that play a key role in recognizing antigens in the vertebrate immune system. TCR and BCR are formed in the stochastic process of V(D)J recombination, creating a diverse sequence repertoire. Progress in high throughput sequencing now allows for deep profiling of T-cell repertoires, by establishing a near-complete list of unique receptor sequences, or "clonotypes," present in a sample.
Comparison of sequenced repertoires has revealed that in any pair of individuals, large numbers of TCR sequences have the same amino acid sequence [1]. Several mechanisms leading to the repertoire overlap have been identified so far. The first mechanism is convergent recombination. Due to biases in V(D)J recombination process, the probability of generation of some receptors is very high, making them appear in almost every individual multiple times and repeatedly sampled in repertoire profiling experiments [2]. This sharing does not result from a common specificity or function of the shared clonotypes, and may in fact correspond to cells from the naive compartment in both donors [3], or from functionally distinct subsets such as CD4 and CD8 T-cells. The second possible reason for TCR sequence sharing is specific to identical twins, who may share T cell clones as a consequence of cord blood exchange in utero via a shared placenta [4]. The third and most interesting mechanism for sharing receptor sequences is convergent selection, in response to a common antigen. From functional studies, such as sequencing of MHC-multimer specific T-cells, it is known that the antigen-specific repertoire is often biased, and the same antigen-specific TCR sequences can be found in different individuals [5][6][7].
Reproducibility of a portion of the antigen-specific Tcell repertoire in different patients creates an opportunity for disease association studies using T-cell repertoire datasets [8,9]. These studies analyse the TCR sequence overlap in large cohorts of healthy controls and patients to identify shared sequences overrepresented in the patient cohort. Here we propose a novel computational method to identify clonotypes which are likely to be shared because of selection for their response to a common antigen, instead of convergent recombination. Our approach is based on a mechanistic model of TCR recombination and is applicable to small cohorts of patients, without the need for a healthy control cohort.
As a proof of concept, we applied our method to two large publicly available TCR beta datasets from Cytomegalovirus (CMV)-positive [9] and type 1 diabetes (T1D) [10] patients. In both studies the authors found shared public TCR clonotypes that are specific to CMVpeptides or self-peptides, respectively. Specificity of these clonotypes was defined using MHC-multimers. We show that clonotypes functionally associated with CMV and T1D in these studies are identified as outliers by our method.
The main ingredient of our approach is to estimate the probability of generation of shared clonotypes, and to use this probability to determine the source of sharing (see Fig. 1 We estimate the theoretical Pgen(σ) for each sequence σ and compare it to P data (σ), which is empirically derived from sharing pattern of that sequence in the cohort. Comparison of these two values allows us to calculate the analog of a p-value, namely the posterior probability that the sharing pattern is explained by convergent recombination alone, with no selection for a common antigen.
any TCR sequencing experiment, chances to sample the same clonotype twice are low, unless this clonotype is easy to generate convergently, with many independent generation events with the same TCR amino acid sequence in each individual (convergent recombination), or if this clonotype underwent clonal expansion, making its concentration in blood high (convergent selection). Thus, we reasoned that convergently selected clonotypes should have a lower generative probability than typical convergently recombined clonotypes. To test this, we estimated the generative probability of the TCR's beta chain Complementarity Determining Region 3 (CDR3) amino-acid sequences that were shared between several patients. Since no algorithm exists that can compute this generative probability directly, our method relies on the random generation and translation of massive numbers of TCR nucleotide sequences using a mechanistic statistical model of V(D)J recombination [11], as can be easily performed e.g. using the IGoR software [12]. In Fig.2A we plot for each clonotype the number of donors sharing that clonotype against its generation probability. Disease-specific TCR variants validated by functional tests in source studies are circled in red. Note that validated disease-specific TCRs have a much lower generation probability than the typical sequences shared by the same number of donors. We developed a method of axis transformation (see SI) to compare the model prediction with data values on the same scale (Fig.2B), so that outliers can be easily identified by their distance to identity line. Our method can be used to narrow down the potential candidates for further experimental validation of responsive receptors.
Our method also identifies other significant outliers than reported in the source studies (shown in red, and obtained after multiple-test correction -see Methods), which may have three possible origins. First, they may be associated with the condition, but were missed by the source studies. Second, they may be due to other factors shared by the patients, such as features involved in thymic or peripheral selection, or reactivity to other common conditions than CMV (e.g. influenza infection). Third, they can be the result of intersample contamination. Our approach is able to diagnose the last explanation by estimating the likelihood of sharing at the level of nucleotide sequences (i.e. synonymously), as detailed in the Methods section.
Our approach can be used on other hypervariable receptor chains (TCR alpha, BCR heavy and light chains), as well as other species (mice, fish, etc.). Recent advances in computational methods allow us to extract TCR repertoires from existing RNA-seq data [13,14]. Huge numbers of available RNAseq datasets from patients with various conditions can be used for analysis and identification of novel virus, cancer, and self reactive TCR variants using our method. The more immunoglobulin receptors with known specificity are found using this type of association mapping, the more clinically relevant information can be extracted from immunoglobulin repertoire data.

ACKNOWLEDGMENTS
This work was supported by Russian Science Foundation grant 15-15-00178, and and partially supported by European Research Council Starting Grant 306312.

Problem formulation
Our framework is applicable to analyze the outcome of a next generation sequencing experiment probing the immune receptor repertoires of n individuals with a given condition, e.g. CMV or Type 1 diabetes. We denote by M i the number of unique amino acid TCR sequences in patient i, i = 1, . . . , N . For a given TCR amino acid sequence σ, we set x i = 1 to indicate that σ is present in patient i's repertoire, and x i = 0 otherwise. For a given shared sequence σ, we want to know how likely its sharing pattern is under the null hypothesis of convergent recombination, correcting for the donors' different sampling depths. In other words, is σ overrepresented in the population of interest? If σ is significantly overrepresented, we also want to quantify the size of this effect.

Overview
Under the null hypothesis, the presence of σ in a certain number of donors is explained by independent convergent V(D)J recombination events in each donor. Given the total number of recombination events that led to the sequenced sample of donor i, N i , the presence of given amino acid sequence σ in donor is Bernoulli distributed with probability where P post (σ) is the model probability that a recombined product found in a blood sample has sequence σ under the null hypothesis. It is formed by the product of P gen (σ), the probability to generate the sequence σ, estimated using a V(D)J recombination model (see I A 3), and Q, a constant correction factor accounting for thymic selection (see I A 4). The number of independent recombination events N i leading to the observed unique sequences in a sample i is unknown, because of convergent recombination events within the sample, but it can be estimated from the number of unique sequences M i , using the model distribution P post (see I A 6). We also calculate the posterior distribution of P data (σ), corresponding to the empirical counterpart of P post (σ) in the cohort, inferred from the sharing pattern of σ across donors. We use information about the presence of σ in our donors, x 1 , . . . , x n and the sequencing depth for each donor, N 1 , . . . , N n (see I A 5), yielding the posterior density: ρ(P data |x 1 , . . . , x N ).
Finally, we estimate the probability, given the observations, that the true value of P data is smaller than the theoretical value P post predicted using V(D)J recombination model, analogous to a p-value and used to identify significant effects: P(P post > P data ) = Ppost 0 ρ(P data |x 1 , . . . , x n )dP data .
(3) To estimate the effect size q(σ) we compare P data to P post , To procedure outlined above requires to calculate P gen (σ), the probability to generate a given CDR3 amino acid sequence. Methods exist to calculate the probability of TCR and BCR nucleotide sequences from a given recombination model [11,12], but are impractical to calculate the probability of amino acid sequences, because of the large number of codon combinations that can lead to the same amino acid sequence, L a=1 n codons (σ(a)), where L is the sequence length, and n codons (τ ) the number of codons coding for amino acid τ . The number is about 1.4 × 10 7 for a typical CDR3 length of 15 amino acid.
Instead, we estimated P gen (σ) using a simple Monte-Carlo approach. We randomly generated a massive number (N sim = 2×10 9 ) of recombination scenarios according to the validated recombination model [11]: The resulting sequences were translated, truncated to only keep the CDR3, and counted. P gen (σ) was approximated by the fraction of events thus generated that led to sequence σ. This approximation becomes more accurate as N sim increases, with an error on P gen (σ) scaling as (P gen (σ)/N sim ) 1/2 .

Estimation of the correction factor Q
Not all generated sequences pass selection in the thymus. P gen systematically underestimates the frequency of recombination event that eventually make it into the observed repertoire. To correct for this effect, we estimate a correction factor Q, as was suggested in [15]: Contrary to [15], which learned a sequence-specific factor for each individual, here we assume that all observed sequences passed thymic selection. Q is a normalization factor accounting for the fact that just a fraction Q −1 of sequences pass thymic selection. This factor is determined for each VJ-combination as an offset when plotting log P gen against log P * data (see I A 5 for definition of P * data ), using least squares fitting.

Estimation of P data (σ), the probability of sequence occurrence in data
The variable x i indicates the presence or absence of a given TCR amino acid sequence σ in the ith dataset with N i recombination events per donor. We want to estimate P data (σ), which is a fraction of recombination events leading to σ in the population of interest. According to Bayes' theorem, for a given σ, the probability density function of P data reads: ρ(P data |x 1 , . . . , x n ) = P(x 1 , . . . , x n |P data )ρ prior (P data ) 1 0 P(x 1 , . . . , x n |P data )ρ prior (P data ) dP data .
The likelihood is given by a product of Bernouilli probabilities: and a flat prior ρ prior (P data ) = const is used.
We estimate P * data (shown in Fig. 2B) as the maximum of the posterior distribution: 6. Estimation of Ni, the number of recombination events The total number N i of recombination events in ith dataset is unknown, but we can count the number of unique CD3 acid sequences M i observed in the sequencing experiment. For a typical TRB experiment, convergent recombination is relatively rare and one could use N i ≈ M i as an approximation. However, for less diverse loci (e.g TRA), or for much higher sequencing depths, one should correct for convergent recombination, as the the observed number of unique aminoacid sequences could be much lower than the actual number of corresponding recombination events.
The average number of unique sequences resulting from N i recombination events is, in theory: where T is the set of sequences that can pass thymic selection. To estimate that number, we generate a very large number N sim of recombinations, leading to N uni unique CDR3 amino acid sequences for which P gen is estimated as explained above. We take T to be a random subset of unique sequences, T ⊂ {σ 1 , . . . , σ Nuni }, of size |T | = N uni /Q, and we apply Eq. 10.
Using this equation we plot the calibration curve for the TRBV5-1 TRBJ2-6 VJ datasets in Fig. 3. For comparison the case of no thymic selection (Q = 1) is shown in red. The inversion of this curve yields N i as a function of M i .

B. Pipeline description
In this section we describe how to apply our algorithm to real data. All the code and data necessary to reproduce our analysis is available online on github (https://github.com/pogorely/vdjRec/). We start with annotated TCR datasets (CDR3 amino acid sequence, V-segment, J-segment), one per donor. Such datasets are produced by MiXCR [16], immunoseq (http://www.adaptivebiotech.com/immunoseq) and most other software for NGS repertoire data preprocessing.
Data we used was in immunoseq format, publicly available from https: //clients.adaptivebiotech.com/immuneaccess database.
We proceed as follows: 1. Split datasets by VJ combinations. The resulting datasets correspond to lists of unique CDR3 amino acid sequences for each donor and VJ combination. All following steps should be done independently for each VJ combination.
2. (Optional). Filter out sequences present in only one donor to speed up the downstream analysis.
3. Generate a large amount of simulated nucleotide TCR sequences for a given VJ combination. Extract and translate their CDR3, and count how many times each sequence appears in the simulated set (restricting to sequences actually observed in donors for better efficiency). The resulting number divided by the total number of simulated sequences is an estimate of P gen .
4. Estimate P * data for each sequence in the dataset, see I A 5.
5. Using P * data and P gen , estimate for each VJ combination the normalization Q by minimizing n j=1 (log P * data (σ j ) − log P gen (σ j ) − log Q) 2 , see I A 4, where σ j , j = 1, . . . , n are the shared sequences.
C. Usage example

Data sources
Data from [9] and [10] is publicly available from the immuneaccess database: https://clients. adaptivebiotech.com/immuneaccess. For our analysis, we only considered VJ combinations for which the authors identified condition-associated clonotypes with MHC-multimer proved specificity. CDR3 aminoacid sequences and V and J segment of these TCR clonotypes are given in Table I.

Analysis results
We applied our pipeline to identify CMV-specific and self-specific TCR sequences listed in Table I CASSLVGGPSSEAFF TRBV05-01 TRBJ1-1 self [10,17]  each dataset we followed our pipeline described in I B.
We found that sequences reported in the source studies as being both significantly enriched in the patient cohort, and antigen-specific according to MHC-multimers, were the most significant in 3 out of 4 datasets. In the remaining TRBV12 dataset, the sequence of interest was the top 40 most significant out of 27, 699 sequences present in at least two CMV-positive donors.

D. Identifying contaminations
Intersample contamination may complicate highthroughput sequencing data analysis in many ways. It could occur both during library preparation or the sequencing process itself [18]. Contaminations have the same nucleotide and amino acid sequence in all datasets, and so our method identifies them as outliers, because their sharing cannot be explained by a high recombina-tion probability.
Our method provides a tool to diagnose contamination. Given an amino-acid sequence present in many donors, we measure its theoretical nucleotide diversity using the same simulation approach we used to calculate the generative probability P gen of the amino acid sequence (see I A 3). If the diversity of the simulated nucleotide sequences is much larger than observed in the data, it is a sign of contamination.
We applied this approach to the CDR3 sequence CASSLVGGPSSEAFF associated to Type 1 diabetes, and found 19 recombination events consistent with that amino acid sequence out of our simulated dataset. We found 18 different nucleotide variants out of the 19 total possible. In contrast, in the data this clononotype had the same nucleotide variant in all of the 8 donors in which it was present. That variant was absent from the simulated set. A one-sided Fisher exact test gives a p < 10 −6 probability of this happening by chance, indicating contamination as a likely source of sharing.