Tools and Resources

Method for identification of condition-associated public antigen receptor sequences

Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Russia
Skolkovo Institute of Science and Technology, Russia
Central European Institute of Technology, Czech republic
Moscow State University, Russia
CNRS, Sorbonne University, Paris-Diderot University, École Normale Supérieure, France
CNRS, Sorbonne University, École Normale Supérieure, France

Mar 13, 2018

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

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: March 28, 2018 (This version)
Accepted Manuscript published: March 13, 2018 (Go to version)
Accepted: March 12, 2018
Received: October 24, 2017

1. Of interest
Generative power of a protein language model trained on multiple sequence alignments

Damiano Sgarbossa, Umberto Lupo, Anne-Florence Bitbol

Research Article Updated Mar 24, 2023
Further reading

Article
Figures and data
Abstract
Introduction
Results
Discussion
Materials and methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

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 one diabetes responsive TCR $β$ sequences .

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

Introduction

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. These receptors consist of two hypervariable chains, the $α$ and $β$ chains in the case of TCR. Progress in high throughput sequencing now allows for deep profiling of TCR $α$ and TCR $β$ chain repertoires, by establishing a near-complete list of unique receptor chain sequences, or ‘clonotypes’, present in a sample. Most sequencing data available correspond to TCR $β$ only, but the same principles discussed below apply to TCR $α$ repertoires, or to paired $α β$ repertoires.

Comparison of sequenced repertoires has revealed that in any pair of individuals, large numbers of TCR $β$ sequences have the same amino acid sequence (Venturi et al., 2011). 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 TCR $β$ sequences is very high, making them appear in almost every individual multiple times and repeatedly sampled in repertoire profiling experiments (Britanova et al., 2014). This sharing does not result from a common specificity or function of T-cells corresponding to the shared TCR $β$ clonotypes, and may in fact correspond to cells from the naive compartment in both donors (Quigley et al., 2010), 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 (Pogorelyy et al., 2017). Note that in that scenario both the $β$ and $α$ chains are shared together. The third and most interesting mechanism for sharing the sequence of either the $β$ or $α$ or both chains 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 $β$ or $α$ chain sequences can be found in different individuals (Miles et al., 2011; Dash et al., 2017; Glanville et al., 2017).

Reproducibility of a portion of the antigen-specific T-cell repertoire in different patients creates an opportunity for disease association studies using TCR $β$ repertoire datasets (Faham et al., 2017; Emerson et al., 2017). 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.

Results

As a proof of concept, we applied our method to two large publicly available TCR $β$ datasets from Cytomegalovirus (CMV)-positive (Emerson et al., 2017) and type one diabetes (T1D) (Seay et al., 2016) patients. In both studies the authors found shared public TCR $β$ clonotypes that are specific to CMV-peptides or self-peptides, respectively. Specificity of these clonotypes was defined using MHC-multimers. We show that TCR $β$ chain sequences 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 Figure 1). Due to the limited sampling depth of any TCR sequencing experiment, chances to sample the same TCR $β$ 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 corresponding T-cell clone 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 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 (Murugan et al., 2012), as can be easily performed for example using the IGoR software (Marcou et al., 2017).

Figure 1

Download asset Open asset

Method principle and pipeline.

(Top left) Sequence overlap between two TCR or BCR repertoires. (Bottom left) There are two major mechanisms for sequence sharing between two repertoires: convergent recombination and convergent selection. Because convergent recombination favors sequences with high generation probabilities, these two classes of sequences have different distributions of the generative probability, $P_{gen} (σ)$ . (Right) We estimate the theoretical $P_{gen} (σ)$ for each sequence $σ$ and compare it to $P_{data} (σ)$ , which is empirically derived from the 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 the convergent recombination alone, with no selection for a common antigen.

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

In Figure 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 TCR $β$ sequences 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 Materials and methods) to compare the model prediction with data values on the same scale (Figure 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. Additional information, like the expansion of the identified TCR $β$ clonotype in the inflammation site, the presence of the same clonotype in the repertoire of activated or memory T-cells, or absence in a cohort of healthy controls, could provide additional evidence for functional association of identified candidates with a given condition.

Figure 2

Download asset Open asset

Identification of condition-associated clonotypes using generative probability

(A) CDR3aa of antigen specific clonotypes (red circles) have less generative probability than other clonotypes shared among the same number of donors. The number of in silico rearrangements obtained for each TCR $β$ sequence in our simulation (which is proportional to generation probability for each clonotype in a given VJ combination $P_{post} (σ)$ ), plotted against the number of patients with that TCR $β$ clonotype. (B) Model prediction of generative probabilities agrees well with data. To directly compare $P_{post} (σ)$ to data, we estimate the empirical probability of occurrence of sequences, $P_{data} (σ)$ , from its sharing pattern across donors (see Materials and methods). In A. and B. red dots indicate significant results (adjusted $P < 0.01$ , Holm’s multiple testing correction), while red circles point to the responsive clonotypes identified in the source studies.

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

Our method also identifies other significant outliers than reported in the source studies (shown in red, and obtained after multiple-test correction – see Materials and 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 Materials and methods section.

Discussion

Antigen receptor sequencing currently has little clinical applications. One of the most important ones is diagnostics and tracking of malignant T-cell and B-cell clones in lymphomas, where it allows for directly measuring the abundances of certain clones at different timepoints. Our method allows for a sequence-based theoretical prediction of T-cell abundances at the population level, and for the identification of T-cell clones associated with infectious and autoimmune conditions. Extensive databases of condition-associated clones can provide a means of disease diagnostics and extend the clinical utility of antigen receptor repertoire sequencing technologies.

This method may also be useful in the analysis of known antigen-specific TCR clonotypes. The typical source of such TCR sequences are MHC-multimer positive cells isolated from one or a few donors (Shugay et al., 2018; Tickotsky et al., 2017). Some of these antigen-specific clonotypes are private, and are hard to find in other patients, providing limited diagnostic value. Our method is able to distinguish these clones from publicly responding clonotypes that are likely to be shared by many patients using only their CDR3 amino acid sequences.

The cohort size necessary for the identification of antigen-specific clonotypes with our method varies (see ‘Designing the experiment’ subsection in Materials and methods). It depends on the strength and diversity of the response to the given antigen. CMV and other Herpesviridae (EBV, HSV), are able to cause a persistent infection, and a large fraction of the TCR repertoire of CMV-positive donors are believed to be specific to them—on average, up to 10% of CD8 +cells are specific to a single CMV epitope in elderly individuals (Khan et al., 2004). However, it was shown that in a human acute infection model of yellow fever vaccination, virus-specific T-cell clones are one of the most abundant in the TCR repertoire and occupy up to 12% of the CD8 +T cell repertoire. This response is short-lived and contracts significantly a month after immunization (Miller et al., 2008). So the peak of an immune response is the best timepoint to search for antigen-specific TCRs in acute infections using this method. T-cell response to herpesviruses is also not unique in terms of public clonotype involvement—in ankylosing spondylitis (Faham et al., 2017), 30–40% of patients share a certain TCR $β$ aminoacid sequence, which is more than the fraction of patients sharing CMV-specific clonotypes that we analysed in this study.

Our approach can be used on other hypervariable receptor chains (TCR $α$ , BCR heavy and light chains), as well as other species (mice, fish, etc.). Both $α$ and $β$ chains contribute to T-cell receptor specificity. Single-cell or paired sequencing technologies (Zemmour et al., 2018) could identify partner receptor chains for condition-associated TCR $α$ or $β$ chain sequences identified with our approach. Antigenic peptides recognized by complete T-cell receptors could then be recovered in vitro using yeast-display libraries of peptide-MHC (Gee et al., 2018). As paired sequencing becomes more widespread, our method can be extended to the analysis of full paired TCR by applying the exact same analysis using the joint recombination probability of $α β$ clonotypes.

Recent advances in computational methods allow us to extract TCR repertoires from existing RNA-Seq data (Bolotin et al., 2017; Brown et al., 2015). Huge numbers of available RNA-Seq 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.

Materials and methods

Statistical analysis

Problem formulation

Request a detailed protocol

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, for example CMV or Type one diabetes. We denote by $M_{i}$ the number of unique amino acid TCR sequences in patient $i$ , $i = 1, \dots, 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

p_{i} = ⟨ x_{i} ⟩ = {(1 - P_{p o s t} (σ))}^{N_{i}},

P_{p o s t} (σ) = P_{g e n} (σ) \times Q,

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 the following subsection), and $Q$ , a constant correction factor accounting for thymic selection (see Estimation of the correction factor $Q$ subsection). 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 Estimation of $N_{i}$ subsection).

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}, \dots, x_{n}$ and the sequencing depth for each donor, $N_{1}, \dots, N_{n}$ (see Estimation of $P_{data} (σ)$ subsubsection), yielding the posterior density: $ρ (P_{data} | x_{1}, \dots, 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_{p o s t} > P_{d a t a}) = \int_{0}^{P_{p o s t}} ρ (P_{d a t a} | x_{1}, \dots, x_{n}) d P_{d a t a} .

To estimate the effect size $q (σ)$ we compare $P_{data}$ to $P_{post}$ ,

q (σ) = \frac{P_{data} (σ)}{P_{post} (σ)} .

Estimation of $P_{gen}$ , the probability of generation of a TCR CDR3 amino acid sequence

Request a detailed protocol

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 (Murugan et al., 2012; Marcou et al., 2017), 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, $\prod_{a = 1}^{L} 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 \times 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 \times 10^{9}$ ) of recombination scenarios according to the validated recombination model (Murugan et al., 2012):

\begin{array}{ll} P_{r e a r r}^{β} & = P (V) P (D, J) P (d e l V | V) P (i n s V D) \\ \times P (d e l D l, d e l D r | D) P (i n s D J) P (d e l J | J) . \end{array}

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$

Request a detailed protocol

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 (Elhanati et al., 2014):

P_{post} (σ) = P_{gen} (σ) \times Q .

Contrary to (Elhanati et al., 2014), 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 the following subsection for definition of $P_{data}^{*}$ ), using least squares fitting.

Estimation of $P_{data} (σ)$ , the probability of sequence occurrence in data

Request a detailed protocol

The variable $x_{i}$ indicates the presence or absence of a given TCR amino acid sequence $σ$ in the $i$ th 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_{d a t a} | x_{1}, \dots, x_{n}) = \frac{P (x_{1}, \dots, x_{n} | P_{d a t a}) ρ_{p r i o r} (P_{d a t a})}{\int_{0}^{1} P (x_{1}, \dots, x_{n} | P_{d a t a}) ρ_{p r i o r} (P_{d a t a}) d P_{d a t a}} .

The likelihood is given by a product of Bernouilli probabilities:

P (x_{1}, \dots, x_{n} | P_{data}) = \prod_{i = 1}^{N} {[1 - {(1 - P_{data})}^{N_{i}}]}^{x_{i}} {[{(1 - P_{data})}^{N_{i}}]}^{1 - x_{i}},

and a flat prior $ρ_{prior} (P_{data}) = const$ is used.

We estimate $P_{data}^{*}$ (shown in Figure 2B) as the maximum of the posterior distribution:

P_{d a t a}^{*} = \arg max_{P_{d a t a}} ρ (P_{d a t a} | x_{1}, \dots, x_{n}) .

Estimation of $N_{i}$ , the number of recombination events

Request a detailed protocol

The total number $N_{i}$ of recombination events in $i$ th 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} \approx 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:

⟨ M_{i} ⟩ = \sum_{σ \in T} {(1 - P_{post} (σ))}^{N_{i}} .

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 \subset {σ_{1}, \dots, σ_{N_{uni}}}$ , of size $| T | = N_{uni} / Q$ , and we apply Equation 8.

Using this equation we plot the calibration curve for the TRBV5-1 TRBJ2-6 VJ datasets in Figure 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}$ .

Figure 3

Download asset Open asset

Calibration curve for TRBV5-1 TRBJ2-6 combination.

Here we plot the fraction of unique amino acid sequences to recombination events against the logarithm of the number of recombination events. The blue line corresponds to the theoretical solution with selection, the red line corresponds to the theoretical solution without selection.

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

Pipeline description

Request a detailed protocol

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/; copy archived at https://github.com/elifesciences-publications/vdjRec/).

We start with annotated TCR datasets (CDR3 amino acid sequence, V-segment, J-segment), one per donor. Such datasets are produced by MiXCR (Bolotin et al., 2015), 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:

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.
(Optional). Filter out sequences present in only one donor to speed up the downstream analysis.
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}$ .
Estimate $P_{data}^{*}$ for each sequence in the dataset, see Estimation of $P_{data} (σ)$ subsection.
Using $P_{data}^{*}$ and $P_{gen}$ , estimate for each VJ combination the normalization $Q$ by minimizing $\sum_{j = 1}^{n} {(\log P_{data}^{*} (σ_{j}) - \log P_{gen} (σ_{j}) - \log Q)}^{2}$ , see Estimation of the correction factor $Q$ subsection, where $σ_{j}$ , $j = 1, \dots, n$ are the shared sequences.
Calculate $P_{post} = Q \times P_{gen}$ . Calculate the p-value (Equation 1) and effect size (Equation 2).

Usage example

Data sources

Request a detailed protocol

Data from (Emerson et al., 2017) and (Seay et al., 2016) 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 1.

Table 1

Published antigen-specific clonotypes used to test the algorithm.

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

CDR3aa	V-segment	J-segment	Antigen source	Ref.
CASSLAPGATNEKLFF	TRBV07-06	TRBJ1-4	CMV	(Emerson et al., 2017)
CASSPGQEAGANVLTF	TRBV05-01	TRBJ2-6	CMV	(Emerson et al., 2017)
CASASANYGYTF	TRBV12-3,−4	TRBJ1-2	CMV	(Emerson et al., 2017)
CASSLVGGPSSEAFF	TRBV05-01	TRBJ1-1	self	(Seay et al., 2016; Gebe et al., 2009)

Analysis results

Request a detailed protocol

We applied our pipeline to identify CMV-specific and self-specific TCR sequences listed in Table 1. For our analysis we used only case cohorts, without controls. For each dataset we followed our pipeline described in Pipeline description subsection. 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 (See Table 2). 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.

Table 2

Output of the algorithm for sequences from Table 1.

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

CDR3aa	V	J	Ag.source .	p-value rank	p-value	Effect size
CASSLAPGATNEKLFF	07–06	1–4	CMV	1/1637	$1.2 \times 10^{- 17}$	8.8
CASSPGQEAGANVLTF	5–01	2–6	CMV	1/5549	$1.8 \times 10^{- 17}$	42.3
CASASANYGYTF	12–3,−4	1–2	CMV	40/27669	$2.5 \times 10^{- 14}$	28.8
CASSLVGGPSSEAFF	5–01	1–1	self	1/2646	$9.5 \times 10^{- 19}$	524

Identifying contaminations

Request a detailed protocol

Intersample contamination may complicate high-throughput sequencing data analysis in many ways. It could occur both during library preparation or the sequencing process itself (Sinha et al., 2017). 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 recombination 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 Estimation of $P_{gen}$ subsection). 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 one 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 eight 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.

Designing the experiment

Request a detailed protocol

Our approach also allows us to obtain important estimates for experiment design. A number of variables affect detection of an antigen-specific clone using our approach: the abundance of the clone in the general population (represented by $P_{data}$ in our approach), the cohort size, the sequencing depth $N_{i}$ in each donor in the cohort, and also the effect size. Fixing any two of these variables results in a constraint between the other two and the affects the probability to detect an antigen-specific clonotype, which translates into the statistical power of the method. As an example of such an analysis, we fix the cohort size at 10, 30 or 100 donors (see Figure 4A–C respectively) and the sequencing depth at $N_{i} = 1000$ unique clones sequenced per repertoire for a given VJ-combination in each donor in the cohort. We ask how frequently a disease specific clone with ${\tilde{P}}_{data}$ abundance in the population and effect size $q = {\tilde{P}}_{data} / P_{post}$ is detected with our method. To address this question for each value ${\tilde{P}}_{data}$ we perform a simulation: we simulate $x_{1}, x_{2}, \dots, x_{n}$ Bernoulli variables, each with a $p_{i} = 1 - e^{- N_{i} {\tilde{P}}_{data}}$ success probability. For a given value of ${\tilde{P}}_{data}$ and $q$ there is a single value of $P_{post} = {\tilde{P}}_{data} / q$ . Then we calculate

Figure 4

Download asset Open asset

Simulation of the method performance with different cohort sizes, sequencing depths, effect sizes and target clone abundances in population.

In panels (**A. B. C**) we plot the number of simulations (out of 100) where a clone with a given effect size $q$ (line color, see legend) and ${\tilde{P}}_{data}$ (x-axis) is found to be significant using our approach, for cohort sizes of 10, 30 and 100 donors respectively. Larger cohort sizes and effect sizes make it possible to resolve clonotypes with lower abundance in the population. In panel (D) we show the effect of sequencing depth for fixed $q = 10$ : larger numbers of clonotypes sequenced per donor allow us to resolve less frequent clones, since a clone of a given ${\tilde{P}}_{data}$ is detected in a larger fraction of donors (panel E).

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

P (P_{p o s t} > P_{d a t a}) = \int_{0}^{P_{p o s t}} ρ (P_{d a t a} | x_{1}, \dots, x_{n}) d P_{d a t a},

where $ρ (P_{data} | x_{1}, \dots, x_{n})$ is the posterior density, and check if $P (P_{p o s t} > P_{d a t a})$ is below a significance threshold of 0.0001. Such a low significant threshold in this example is chosen to take into account the multiple testing correction: we assume that about 1000 shared clones would be tested in a such analysis and p¡0.01 after multiple testing is chosen as the significance threshold in this study, which gives p¡0.0001 before the Bonferroni multiple testing correction. Then we plot the number of simulations in which a significant result was obtained for given effect size $q$ and ${\tilde{P}}_{data}$ for the clone of interest and the fraction of donors with this sequence in the simulated cohort (see Figure 4E, blue curve). Unsurprisingly, the effect size plays a role in the probability to detect an antigen specific clone, and the detection is not possible at all if the clone is not shared between several donors in the cohort (in our example this happens for ${\tilde{P}}_{data} < 10^{- 5}$ ) irrespective to the effect size. Larger cohort sizes can help to resolve clones with lower abundances, but sequencing depth also has a strong effect on the power of the approach. In Figure 4D and E we show simulation results for a fixed $q = 10$ and different sequencing depths $N_{i}$ of 100, 1000 or 10000 clones per donor in a given VJ combination. Interestingly, a large sequencing depth (black curve) can lead to a situation when an abundant and frequently generated clone will not be detected by the algorithm, because it will be found in all donors in the cohort. An additional test that checks the predictions by lowering the sequencing depth in silico by downsampling can solve this problem.

Another complicated question is how $P_{data}$ is related to the number of clones and the fraction of the repertoire involved in the response to the infection in a given donor. If the same antigen-specific clone is present in every donor, $P_{data}$ is close to the average abundance of this clone in the repertoire. However one can imagine an opposite situation where the response is so diverse and private that different clones respond to a given antigen in each donor. It was previously shown that the diversity and publicness of responding T-cell clonotypes varies a lot across antigens (Dash et al., 2017). Our approach is restricted to the identification of public antigen-specific clonotypes, which may not exist for all antigens.

Data availability

The following previously published data sets were used

1. Emerson RO
2. DeWitt WS
3. Vignali M
4. Gravley J
5. Hu JK
6. Osborne EJ
7. Desmarais C
8. Klinger M
9. Carlson CS
10. Hansen JA
11. Rieder M
12. Robins HS
(2017) Immunosequencing identifies signatures of cytomegalovirus exposure history and HLA-mediated effects on the T cell repertoire
Publicly available at ImmuneAccess.

http://dx.doi.org/10.21417/B7001Z
1. Seay HR
2. Yusko E
3. Rothweiler SJ
4. Zhang L
5. Posgai AL
6. Campbell-Thompson M
7. Vignali M
8. Emerson RO
9. Kaddis JS
10. Ko D
11. Nakayama M
12. Smith MJ
13. Cambier JC
14. Pugliese A
15. Atkinson MA
16. Robins HS
17. Brusko TM
(2016) Tissue distribution and clonal diversity of the T and B cell repertoire in type 1 diabetes
Publicly available at ImmuneAccess.

http://doi.org/10.21417/B73S3K

References

(2017) Antigen receptor repertoire profiling from RNA-seq data
Nature Biotechnology 35:908–911.
https://doi.org/10.1038/nbt.3979
- PubMed
- Google Scholar
(2015) MiXCR: software for comprehensive adaptive immunity profiling
Nature Methods 12:380–381.
https://doi.org/10.1038/nmeth.3364
- PubMed
- Google Scholar
(2014) Age-related decrease in TCR repertoire diversity measured with deep and normalized sequence profiling
The Journal of Immunology 192:2689–2698.
https://doi.org/10.4049/jimmunol.1302064
- PubMed
- Google Scholar
(2015) Profiling tissue-resident T cell repertoires by RNA sequencing
Genome Medicine 7:125.
https://doi.org/10.1186/s13073-015-0248-x
- PubMed
- Google Scholar
1. Dash P
2. Fiore-Gartland AJ
3. Hertz T
4. Wang GC
5. Sharma S
6. Souquette A
7. Crawford JC
8. Clemens EB
9. Nguyen THO
10. Kedzierska K
11. La Gruta NL
12. Bradley P
13. Thomas PG
(2017) Quantifiable predictive features define epitope-specific T cell receptor repertoires
Nature 547:89–93.
https://doi.org/10.1038/nature22383
- PubMed
- Google Scholar
1. Elhanati Y
2. Murugan A
3. Callan CG
4. Mora T
5. Walczak AM
(2014) Quantifying selection in immune receptor repertoires
PNAS 111:9875–9880.
https://doi.org/10.1073/pnas.1409572111
- PubMed
- Google Scholar
1. Emerson RO
2. DeWitt WS
3. Vignali M
4. Gravley J
5. Hu JK
6. Osborne EJ
7. Desmarais C
8. Klinger M
9. Carlson CS
10. Hansen JA
11. Rieder M
12. Robins HS
(2017) Immunosequencing identifies signatures of cytomegalovirus exposure history and HLA-mediated effects on the T cell repertoire
Nature Genetics 49:659–665.
https://doi.org/10.1038/ng.3822
- PubMed
- Google Scholar
1. Faham M
2. Carlton V
3. Moorhead M
4. Zheng J
5. Klinger M
6. Pepin F
7. Asbury T
8. Vignali M
9. Emerson RO
10. Robins HS
11. Ireland J
12. Baechler-Gillespie E
13. Inman RD
(2017) Discovery of T cell receptor β motifs specific to HLA-B27-positive ankylosing spondylitis by deep repertoire sequence analysis
Arthritis & Rheumatology 69:774–784.
https://doi.org/10.1002/art.40028
- PubMed
- Google Scholar
1. Gebe JA
2. Yue BB
3. Unrath KA
4. Falk BA
5. Nepom GT
(2009) Restricted autoantigen recognition associated with deletional and adaptive regulatory mechanisms
The Journal of Immunology 183:59–65.
https://doi.org/10.4049/jimmunol.0804046
- PubMed
- Google Scholar
1. Gee MH
2. Han A
3. Lofgren SM
4. Beausang JF
5. Mendoza JL
6. Birnbaum ME
7. Bethune MT
8. Fischer S
9. Yang X
10. Gomez-Eerland R
11. Bingham DB
12. Sibener LV
13. Fernandes RA
14. Velasco A
15. Baltimore D
16. Schumacher TN
17. Khatri P
18. Quake SR
19. Davis MM
20. Garcia KC
(2018) Antigen identification for orphan T cell receptors expressed on tumor-infiltrating lymphocytes
Cell 172:1–15.
https://doi.org/10.1016/j.cell.2017.11.043
- PubMed
- Google Scholar
1. Glanville J
2. Huang H
3. Nau A
4. Hatton O
5. Wagar LE
6. Rubelt F
7. Ji X
8. Han A
9. Krams SM
10. Pettus C
11. Haas N
12. Arlehamn CSL
13. Sette A
14. Boyd SD
15. Scriba TJ
16. Martinez OM
17. Davis MM
(2017) Identifying specificity groups in the T cell receptor repertoire
Nature 547:94–98.
https://doi.org/10.1038/nature22976
- PubMed
- Google Scholar
1. Khan N
2. Hislop A
3. Gudgeon N
4. Cobbold M
5. Khanna R
6. Nayak L
7. Rickinson AB
8. Moss PA
(2004) Herpesvirus-specific CD8 T cell immunity in old age: cytomegalovirus impairs the response to a coresident EBV infection
The Journal of Immunology 173:7481–7489.
https://doi.org/10.4049/jimmunol.173.12.7481
- PubMed
- Google Scholar
Preprint
(2017) IGoR: a tool for high-throughput immune repertoire analysis
bioRxiv.

https://arxiv.org/abs/1705.08246
- Google Scholar
(2011) Bias in the αβ T-cell repertoire: implications for disease pathogenesis and vaccination
Immunology and Cell Biology 89:375–387.
https://doi.org/10.1038/icb.2010.139
- PubMed
- Google Scholar
1. Miller JD
2. van der Most RG
3. Akondy RS
4. Glidewell JT
5. Albott S
6. Masopust D
7. Murali-Krishna K
8. Mahar PL
9. Edupuganti S
10. Lalor S
11. Germon S
12. Del Rio C
13. Mulligan MJ
14. Staprans SI
15. Altman JD
16. Feinberg MB
17. Ahmed R
(2008) Human effector and memory CD8+ T cell responses to smallpox and yellow fever vaccines
Immunity 28:710–722.
https://doi.org/10.1016/j.immuni.2008.02.020
- PubMed
- Google Scholar
1. Murugan A
2. Mora T
3. Walczak AM
4. Callan CG
(2012) Statistical inference of the generation probability of T-cell receptors from sequence repertoires
PNAS 109:16161–16166.
https://doi.org/10.1073/pnas.1212755109
- PubMed
- Google Scholar
1. Pogorelyy MV
2. Elhanati Y
3. Marcou Q
4. Sycheva AL
5. Komech EA
6. Nazarov VI
7. Britanova OV
8. Chudakov DM
9. Mamedov IZ
10. Lebedev YB
11. Mora T
12. Walczak AM
(2017) Persisting fetal clonotypes influence the structure and overlap of adult human T cell receptor repertoires
PLoS Computational Biology 13:e1005572.
https://doi.org/10.1371/journal.pcbi.1005572
- PubMed
- Google Scholar
1. Quigley MF
2. Greenaway HY
3. Venturi V
4. Lindsay R
5. Quinn KM
6. Seder RA
7. Douek DC
8. Davenport MP
9. Price DA
(2010) Convergent recombination shapes the clonotypic landscape of the naive T-cell repertoire
PNAS 107:19414–19419.
https://doi.org/10.1073/pnas.1010586107
- PubMed
- Google Scholar
1. Seay HR
2. Yusko E
3. Rothweiler SJ
4. Zhang L
5. Posgai AL
6. Campbell-Thompson M
7. Vignali M
8. Emerson RO
9. Kaddis JS
10. Ko D
11. Nakayama M
12. Smith MJ
13. Cambier JC
14. Pugliese A
15. Atkinson MA
16. Robins HS
17. Brusko TM
(2016) Tissue distribution and clonal diversity of the T and B cell repertoire in type 1 diabetes
JCI Insight 1:1–19.
https://doi.org/10.1172/jci.insight.88242
- Google Scholar
1. Shugay M
2. Bagaev DV
3. Zvyagin IV
4. Vroomans RM
5. Crawford JC
6. Dolton G
7. Komech EA
8. Sycheva AL
9. Koneva AE
10. Egorov ES
11. Eliseev AV
12. Van Dyk E
13. Dash P
14. Attaf M
15. Rius C
16. Ladell K
17. McLaren JE
18. Matthews KK
19. Clemens EB
20. Douek DC
21. Luciani F
22. van Baarle D
23. Kedzierska K
24. Kesmir C
25. Thomas PG
26. Price DA
27. Sewell AK
28. Chudakov DM
(2018) VDJdb: a curated database of T-cell receptor sequences with known antigen specificity
Nucleic Acids Research 46:D419–D427.
https://doi.org/10.1093/nar/gkx760
- PubMed
- Google Scholar
Preprint
1. Sinha R
2. Stanley G
3. Gulati GS
4. Ezran C
5. Travaglini KJ
6. Wei E
7. Chan CKF
8. Nabhan AN
9. Su T
10. Morganti RM
11. Conley SD
12. Chaib H
13. Red-Horse K
14. Longaker MT
15. Snyder MP
16. Krasnow MA
17. Weissman IL
(2017) Index switching causes “spreading-of-signal” among multiplexed samples in illumina hiseq 4000 dna sequencing
bioRxiv.
https://doi.org/10.1101/125724
- Google Scholar
(2017) McPAS-TCR: a manually curated catalogue of pathology-associated T cell receptor sequences
Bioinformatics 33:2924–2929.
https://doi.org/10.1093/bioinformatics/btx286
- PubMed
- Google Scholar
1. Venturi V
2. Quigley MF
3. Greenaway HY
4. Ng PC
5. Ende ZS
6. McIntosh T
7. Asher TE
8. Almeida JR
9. Levy S
10. Price DA
11. Davenport MP
12. Douek DC
(2011) A mechanism for TCR sharing between T cell subsets and individuals revealed by pyrosequencing
The Journal of Immunology 186:4285–4294.
https://doi.org/10.4049/jimmunol.1003898
- PubMed
- Google Scholar
1. Zemmour D
2. Zilionis R
3. Kiner E
4. Klein AM
5. Mathis D
6. Benoist C
(2018) Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR
Nature Immunology 19:291–301.
https://doi.org/10.1038/s41590-018-0051-0
- PubMed
- Google Scholar

Decision letter

Arup K Chakraborty

Reviewing Editor; Massachusetts Institute of Technology, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Method for identification of condition-associated public antigen receptor sequences" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Arup Chakraborty as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor/Senior Editor has drafted this decision to help you prepare a revised submission.

Your manuscript describes a computational pipeline that identifies particular TCR sequences that are over represented in groups of individuals. In effect, it identifies the outlier TCR sequences populations from a computationally determined normal distribution. You go on to show that, based on TCR deep sequencing datasets, this method can retrospectively identify TCRβ chains that are over represented following CMV infections, and TCRβ sequences that have reactivity to defined autoantigens in type 1 diabetes. The manuscript is theoretically interesting, useful, timely, and clearly written. However, we have some significant concerns about the general applicability of this tool and its value for prospective, rather than retrospective, studies. These concerns are listed below and need to be addressed in a revision.

1) Methods to identify outlier populations are understood and statistical approaches to achieve this goal are available. Indeed, a related approach to this question is the basis for lymphoid tumor diagnostic methodologies being pursued by several companies. In what ways is your tool different/better compared to the available methods?

2) Your work is based on small cohorts of patients. HCMV can be responsible for a disproportionate number of TCR clonotypes and can have a stable profile over many years, even by the standards of other chronic infections. So, we wonder whether an attractive feature of this tool, which is its ability to glean information from data on small cohorts will carry over beyond the HCMV setting. Would one need a larger cohort for most disease applications? Is there a theoretical way to estimate what part of the repertoire (in terms of unique clones and their frequencies) has to be devoted to an infectious disease in order for it to be resolved in a given cohort size (maybe a scaling rule)?

3) A general concern is whether this tool will have broad impact. You used TCR sequence data sets of HLA-matched patients, which had already been shown to have T cell expansion to a defined peptide/HLA combination. The approach works if one has HLA-matched subjects, because one would assume in general that different HLA alleles will present a different spectrum of peptides, and that TCR Vgene – HLA interactions impact ligand specificity. Further, a significant element to the approach only works retrospectively; without prior knowledge of the TCR specificity one does not know whether the identified outlier populations are disease-relevant. For example, in Figure 2, there are many outlier sequences (red), which may be noise, may be disease relevant or may be consequences of HLA-biased selection that was not accounted for by the authors computationally-defined normal TCRb distribution. Can the method be useful for prospective studies, or learn new biology in retrospective studies? Perhaps one way to address this question is the following:

Is there any other viral infection for which comparable data is available? For instance, it would be nice if data on viruses like HSV or EBV were available along with data related to acute infection or vaccination (where a method such as this, if applicable, would be incredibly valuable).

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Method for identification of condition-associated public antigen receptor sequences" for further consideration at eLife. Your revised article has been favorably evaluated by Arup Chakraborty (Senior editor), a Reviewing editor, and 1 reviewer.

You present an important statistical tool for analyzing TCR sequence diversity. We will soon get to a point where entire human TCR repertoires will be sequenced and then the question is how do we understand this data to assess immune health – your method is an advance that allows us to take steps toward this goal, and is therefore a valuable resource/tool. However, you often confuse TCR with clonotype in the writing of your paper. A TCR β chain sequence is not a TCR clonotype. Throughout the manuscript, you describe what you are analyzing as TCR clonotypes when you are actually referring to independent TCR β sequences. An independent TCR α chain or a TCR β chain sequence is NOT a receptor. Once this point is consistently clarified in the manuscript, your paper will be accepted.

Reviewer #1:

Re-review:

The major concern I have with the method, writing and the inability of this platform to provide a framework for prospective studies is the author's miss-statements (I assume they know better) that a TCR β chain sequence is not a TCR clonotype. Throughout the manuscript, they authors describe what they are analyzing as TCR clonotypes when they are actually referring to independent TCR β sequences. An independent TCRa chain or a TCRb chain sequence is NOT a receptor. The CDR3beta sequence does not equal or even mark a particular T cell clone.

Problem statement that this creates:

Example in the Introduction: the authors argue that there are "Several mechanisms leading to the repertoire overlap. The first mechanism is convergent recombination[…]" This first statement or option is a product of sequencing just the TCRb chain, i.e., there are common rearrangements that occur do to V-J or V-D-J rearrangements that have no or few N-region additions. The second and third "options" are indeed TCRa + TCRb clonotypes, comparing the first possibility to 2 and 3 is comparing apples to oranges.

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 Britanova et al., (2014). 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 Quigley et al., (2010), 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 Pogorelyy et al., (2017). The third and most interesting mechanism for sharing receptor sequences is convergent selection, in response to a common antigen.

- The terminology makes these types of statements nonsensical. Identical TCR clonotypes will absolutely have the identical antigen-specificity.

Issues that miss-representing TCRb sequences for TCR clonotypes in the general approach being promoted:

The incorrect use of TCR clonotype leads to some particularly difficult to imagine scenarios. For examples, in the reviewer response and within the manuscript, "our method can be used to narrow down the potential candidates for further experimental validation of responsive receptors." (i.e., identify the antigen being recognized). Their idea regarding this approach is to use: "Functional tests (like cultivation with peptides, or cytokine secretion assays) are the ultimate way to confirm specificity of these predicted clonotypes."

Within the cover letter: "We also are careful to note that this method should be used to identify antigen-specific candidates that need to be further verified by other methods. Nevertheless, we believe that identifying candidates for these experiments is extremely useful."

- However, the sequencing and analyses method of the manuscript does not identify the V α chain, with only half of the receptor there is nothing to study further. Thus, there is no method provided that one could identify interesting receptors for prospective studies.

In summary, this is a method to identify outlier populations, for which retrospective data can be superimposed to create a snapshot of an immune response that has already been described using for example pMHC tetramers or other methods.

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

Author response

Your manuscript describes a computational pipeline that identifies particular TCR sequences that are over represented in groups of individuals. In effect, it identifies the outlier TCR sequences populations from a computationally determined normal distribution. You go on to show that, based on TCR deep sequencing datasets, this method can retrospectively identify TCRβ chains that are over represented following CMV infections, and TCRβ sequences that have reactivity to defined autoantigens in type 1 diabetes. The manuscript is theoretically interesting, useful, timely, and clearly written. However, we have some significant concerns about the general applicability of this tool and its value for prospective, rather than retrospective, studies. These concerns are listed below and need to be addressed in a revision.

1) Methods to identify outlier populations are understood and statistical approaches to achieve this goal are available. Indeed, a related approach to this question is the basis for lymphoid tumor diagnostic methodologies being pursued by several companies. In what ways is your tool different/better compared to the available methods?

Current approaches for lymphoid tumor and minimal residual disease (MRD) diagnostics rely on longitudinal tracking of certain malignant rearrangements before and after treatment. These technologies use the abundance of the malignant clones in the organism before and after the treatment, which is measured directly with NGS, qPCR or flow cytometry. The sequences of the malignant clones are not linked to a known function. In our method we use a sequence-based theoretical estimate of a TCR clone’s abundance in the population to identify candidates for functional receptors that are convergently selected in different patients in response to a specific disease, which is a different task than the direct longitudinal tracking of clonotypes. However, there is a possible lesson to be learned from the overlap between the two methods: tracking malignant clones that have a high recombination probability may not be reliable because one could find an independently recombined healthy clone with the same sequence and interpret it as malignant. Conversely, the TCR recombination probability may be used to identify malignant TCR clones not suitable for longitudinal tracking due to their high probability to be recombined multiple times (Nazarov et al., 2016). If such clonotypes are found after chemotherapy in lymphoma patients, they could be surviving malignant clonotypes as usually assumed, but alternatively they could be different healthy clones with exactly the same independently recombined TCR sequence. In this case it is better to use different MRD markers and the antigen receptor generation probability could help to identify this situation.

We added the following text to clarify the features of our method (Discussion section):

“Antigen receptor sequencing currently has little clinical applications. One of the most important ones is diagnostics and tracking of malignant T-cell and B-cell clones in lymphomas, where it allows for directly measuring the abundances of certain clones at different timepoints. Our method allows for a sequence-based theoretical prediction of T-cell abundances at the population level, and for the identification of T-cell clones associated with infectious and autoimmune conditions. Extensive databases of condition-associated clones can provide a means of disease diagnostics and extend the clinical utility of antigen receptor repertoire sequencing technologies.”

2) Your work is based on small cohorts of patients. HCMV can be responsible for a disproportionate number of TCR clonotypes and can have a stable profile over many years, even by the standards of other chronic infections. So, we wonder whether an attractive feature of this tool, which is its ability to glean information from data on small cohorts will carry over beyond the HCMV setting. Would one need a larger cohort for most disease applications?

Data from published studies suggests that CMV is not a unique condition both in terms of the response magnitude and the number of donors sharing disease specific clones. Faham et al., have shown that 30–40 percent of donors with ankylosing spondylitis share a certain aminoacid TCR sequence, which is absent in healthy controls. This sharing number is larger than the fraction of donors sharing a CMV-specific clone in our current study. We would have liked to analyze different datasets, but we could not access other datasets with suitable cohort sizes and sequencing depth (see below for a discussion of this point).

One can also extend this method to acute infections. It is known that CMV-related clonotypes can occupy a large fraction of the repertoire for a long time [Khan et al., 2004]. However, in an acute viral infection model of yellow fever immunization it was shown that YF-reactive clones that are undetected in the repertoire before immunization can occupy up to 10 percent of the repertoire at the peak of the response (Miller et al., 2008, DeWitt et al., 2015). The response to acute infections is much less stable and strongly decreases on short timescales, but one may overcome this limitation by collecting samples at the peak of the immune response to an infection.

We added the following sentences into the text to clarify this point (Discussion section):

“The cohort size necessary for the identification of antigen-specific clonotypes with our method varies (see “Designing the experiment” subsection in Methods). It depends on the strength and diversity of the response to the given antigen. […] T-cell response to herpesviruses is also not unique in terms of public clonotype involvement – in ankylosing spondylitis (Faham et al., 2017), 30–40 percent of patients share a certain TCR aminoacid sequence, which is more than the fraction of patients sharing CMV-specific clonotypes that we analysed in this study.”

Is there a theoretical way to estimate what part of the repertoire (in terms of unique clones and their frequencies) has to be devoted to an infectious disease in order for it to be resolved in a given cohort size (maybe a scaling rule)?

This is a very interesting question, which we now address with a new section and discussion. In short, the answer depends on both the cohort size and the sequencing depth. To simplify the experiment design, we added an additional pipeline into the software, which is able to tell if a clone of given abundance in the population of patients would be resolved in a given cohort size.

We added an additional subsection in the Materials and methods section (“Designing the experiment”) and a new figure (see Figure 4), demonstrating examples of such analysis.

3) A general concern is whether this tool will have broad impact. You used TCR sequence data sets of HLA-matched patients, which had already been shown to have T cell expansion to a defined peptide/HLA combination. The approach works if one has HLA-matched subjects, because one would assume in general that different HLA alleles will present a different spectrum of peptides, and that TCR Vgene – HLA interactions impact ligand specificity. Further, a significant element to the approach only works retrospectively; without prior knowledge of the TCR specificity one does not know whether the identified outlier populations are disease-relevant. For example, in Figure 2, there are many outlier sequences (red), which may be noise, may be disease relevant or may be consequences of HLA-biased selection that was not accounted for by the authors computationally-defined normal TCRb distribution. Can the method be useful for prospective studies, or learn new biology in retrospective studies?

Our study may seem retrospective, because we chose examples where the right answer (public clonotypes with proved specificity) is known. This was necessary to validate the output of the algorithm. In 3 experiments out of 4 the correct (meaning validated in the source study) clonotype was the most significant signal in our analysis. We also are careful to note that this method should be used to identify antigen-specific candidates that need to be further verified by other methods. Nevertheless, we believe that identifying candidates for these experiments is extremely useful.

Functional tests (like cultivation with peptides, or cytokine secretion assays) are the ultimate way to confirm specificity of these predicted clonotypes. Clonal expansion of identified clonotypes at the inflammation site, presence in the repertoire of activated/memory T-cells, or absence in the control cohort repertoire, may provide additional evidence for condition-association of a given candidate clone.

We added the following sentence into the main text (Results section):

“Additional information, like the expansion of the identified clone in the inflammation site, the presence of the same clonotype in the repertoire of activated or memory T-cells, or absence in a cohort of healthy controls, could provide additional evidence for functional association of identified candidates with a given condition.”

Additionally, there is an interest in purely retrospective studies. In retrospective studies (where the specificity of TCR clones is already defined using, for example, MHC-multimers) it may be useful to predict the abundance of the known antigen-specific clones in the general population. For instance, if a researcher isolated and sequenced antigen-specific T-cells from a patient, some of these T-cell clones could be private and may not be found in other people, while some of them may be public, and easily found in many more yet unstudied patients. Such public clones are much more valuable for diagnostics, and they can be identified using only their TCR aminoacid sequence in such retrospective studies.

We added the following sentences to the main text to clarify this point (Discussion section):

“This method may also be useful in the analysis of known antigen-specific TCR clonotypes. The typical source of such TCR sequences are MHC-multimer positive cells isolated from one or a few donors (Shugay et al., 2017., Tickotsky et al., 2017). Some of these antigen-specific clonotypes are private, and are hard to find in other patients, providing limited diagnostic value. Our method is able to distinguish these clones from publicly responding clonotypes that are likely to be shared by many patients using only their CDR3 amino acid sequence.”

Perhaps one way to address this question is the following:

Is there any other viral infection for which comparable data is available? For instance, it would be nice if data on viruses like HSV or EBV were available along with data related to acute infection or vaccination (where a method such as this, if applicable, would be incredibly valuable).

To our knowledge, only one study on response to acute infection is available: DeWitt et al., 2015. We have tried to analyze this dataset but unfortunately, only 9 donors are available in this study, which is not enough to get significant results with our method. The authors also do not report sequences of significantly expanded clonotypes, so it is not possible to validate the output of the algorithm.

As a rule of thumb, assuming currently used typical sequencing depth, our method is applicable for identification of shared TCR sequences in vaccination and acute infection data, but cohorts of at least 20–30 donors are required for such analysis.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

You present an important statistical tool for analyzing TCR sequence diversity. We will soon get to a point where entire human TCR repertoires will be sequenced and then the question is how do we understand this data to assess immune health – your method is an advance that allows us to take steps toward this goal, and is therefore a valuable resource/tool. However, you often confuse TCR with clonotype in the writing of your paper. A TCR β chain sequence is not a TCR clonotype. Throughout the manuscript, you describe what you are analyzing as TCR clonotypes when you are actually referring to independent TCR β sequences. An independent TCR α chain or a TCR β chain sequence is NOT a receptor. Once this point is consistently clarified in the manuscript, your paper will be accepted.

We agree that our terminology was not specific enough and may confuse the readers. We have updated the manuscript to be more specific about single α or β chains and complete α/β TCR clonotypes. See our point-by-point response and implemented change below

Reviewer #1:

Re-review:

The major concern I have with the method, writing and the inability of this platform to provide a framework for prospective studies is the author's miss-statements (I assume they know better) that a TCR β chain sequence is not a TCR clonotype. Throughout the manuscript, they authors describe what they are analyzing as TCR clonotypes when they are actually referring to independent TCR β sequences. An independent TCRa chain or a TCRb chain sequence is NOT a receptor. The CDR3beta sequence does not equal or even mark a particular T cell clone.

We completely agree with the reviewer and have updated the manuscript to clarify this point.

We have edited the first paragraph of the Introduction:

“These receptors consist of two hypervariable chains, the α and β chains in the case of TCR. Progress in high throughput sequencing now allows for deep profiling of TCRalpha and TCRbeta chain repertoires, by establishing a near-complete list of unique receptor chain sequences, or ‘’clonotypes’’, present in a sample. Most sequencing data available corresponds to TCRbeta, but the same principles discussed below apply to TCRalpha repertoires or to paired α β repertoires.”

We also specified when we talk about single T-cell receptor chains specifically throughout the text.

Problem statement that this creates:

Example in the Introduction: the authors argue that there are "Several mechanisms leading to the repertoire overlap. The first mechanism is convergent recombination…" This first statement or option is a product of sequencing just the TCRb chain, i.e., there are common rearrangements that occur do to V-J or V-D-J rearrangements that have no or few N-region additions. The second and third "options" are indeed TCRa + TCRb clonotypes, comparing the first possibility to 2 and 3 is comparing apples to oranges.

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 Britanova et al., (2014). 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 Quigley et al., (2010), 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 Pogorelyy et al., (2017). The third and most interesting mechanism for sharing receptor sequences is convergent selection, in response to a common antigen.

- The terminology makes these types of statements nonsensical. Identical TCR clonotypes will absolutely have the identical antigen-specificity.

We agree with the reviewer that case 2 (physical sharing of cells) implies the sharing of complete α/β TCR sequence. In case 3 (selection by affinity to the same antigen), the situation varies for different epitopes and in some cases tetramer-sorted T-cells of the same specificity have been reported to have the same TCRbeta but very different TCRalpha, and vice versa (see examples of paired sequences from (Dash et al., 2017), deposited to vdjdb.cdr3.net)). Therefore, the sharing of single receptor chain sequences due to selection by affinity to the same antigen is generally more likely than sharing of complete α/β TCR.

We have edited the second paragraph of the introduction to clarify these points in the following way:

“Comparison of sequenced repertoires has revealed that in any pair of individuals, large numbers of TCRbeta sequences have the same amino acid sequence Venturi et al., (2011). […] 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 β or α chain sequences can be found in different individuals Miles et al., (2011); Dash et al., (2017); Glanville et al., (2017).”

Issues that miss-representing TCRb sequences for TCR clonotypes in the general approach being promoted:

The incorrect use of TCR clonotype leads to some particularly difficult to imagine scenarios. For examples, in the reviewer response and within the manuscript, "our method can be used to narrow down the potential candidates for further experimental validation of responsive receptors." (i.e., identify the antigen being recognized). Their idea regarding this approach is to use: "Functional tests (like cultivation with peptides, or cytokine secretion assays) are the ultimate way to confirm specificity of these predicted clonotypes."

Within the cover letter: "We also are careful to note that this method should be used to identify antigen-specific candidates that need to be further verified by other methods. Nevertheless, we believe that identifying candidates for these experiments is extremely useful."

We agree with reviewer that the mere TCR β chain or the mere TCR α chain does not determine complete receptor specificity. The TCRbeta sequence is definitely not enough to determine receptor specificity in vitro. However, if a β-chain or an α-chain shared in patient cohort is identified as significant by our approach it is likely to be a part of T-cell receptor with the same specificity in each patient. T-cells from patients with the largest expansion of identified TCRbeta or TCRalpha clonotypes could then be used for single cell RNA sequencing to identify partner receptor chains or for various functional assays with candidate antigens.

We added sentences to the Discussion section on the contribution of both receptor chains to TCR specificity and possible assays to identify second receptor chain and possible antigen, see below.

- However, the sequencing and analyses method of the manuscript does not identify the V α chain, with only half of the receptor there is nothing to study further. Thus, there is no method provided that one could identify interesting receptors for prospective studies.

We agree that complete antigen-specific α/β T-cell receptor sequences would be much more informative and practically useful for prospective studies. Unfortunately, technologies for paired sequencing of α and β TCR chains (i.e. single-cell RNAseq) are still relatively low-throughput and expensive, comparing to the conventional RepSeq study targeting only single receptor chains. TCRalpha and TCRbeta pairing information is also impossible to retrieve from the bulk RNA-seq data. There are no published large datasets of complete paired α/β TCR sequences for cohorts of patients with the same condition. However, if such datasets were to become available in the future, our approach could be used directly to identify candidate complete receptor sequences there. The idea would be the same: the presence of shared TCR α chains, or β chains, or paired α β clonotypes with low recombination probability in many patients is the consequence of clonal expansion after recognition of the same antigenic peptide, and paired TCRbeta-TCRalpha datasets would allow us to retrieve the partner α chain (for TCRbeta condition-associated clonotypes, and vice versa for TCRalpha analysis) from each donor, yielding several complete α/β candidate receptors.

We also believe that even single condition-associated TCRbeta or TCRalpha chain sequences could be useful for diagnostic purposes, as was demonstrated by (Emerson et al., 2017) for CMV-associated TCRbeta clonotypes.

We have added the following text to the Discussion section to emphasize that single receptor chain does not fully determine TCR specificity. We also reference recently emerged methods to identify the second receptor chain, and to determine specificity of the complete α/β TCR, which might be useful in combination with our approach:

“Our approach can be used on other hypervariable receptor chains (TCR α, BCR heavy and light chains), as well as other species (mice, fish, etc.). […] As paired sequencing becomes more widespread, our method can be extended to the analysis of full paired TCR by applying the exact same analysis using the joint recombination probability of α β clonotypes.”

In summary, this is a method to identify outlier populations, for which retrospective data can be superimposed to create a snapshot of an immune response that has already been described using for example pMHC tetramers or other methods.

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

Article and author information

Author details

Mikhail V Pogorelyy

Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Moscow, Russia

Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0773-1204
Anastasia A Minervina

Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Moscow, Russia

Contribution
Data curation, Software, Investigation

Competing interests
No competing interests declared
Dmitriy M Chudakov
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Moscow, Russia
2. Center for Data-Intensive Biomedicine and Biotechnology, Skolkovo Institute of Science and Technology, Moscow, Russia
3. Central European Institute of Technology, Brno, Czech republic
Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration, Writing—review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0430-790X
Ilgar Z Mamedov

Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Moscow, Russia

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration, Writing—review and editing

Competing interests
No competing interests declared
Yuri B Lebedev
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Moscow, Russia
2. Biological Faculty, Moscow State University, Moscow, Russia
Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration, Writing—review and editing

Contributed equally with
Thierry Mora and Aleksandra M Walczak

For correspondence
lebedev_yb@mx.ibch.ru

Competing interests
No competing interests declared
Thierry Mora

Laboratoire de Physique Statistique, CNRS, Sorbonne University, Paris-Diderot University, École Normale Supérieure, Paris, France

Contribution
Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing

Contributed equally with
Yuri B Lebedev and Aleksandra M Walczak

For correspondence
tmora@lps.ens.fr

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5456-9361
Aleksandra M Walczak

Laboratoire de Physique Theorique, CNRS, Sorbonne University, École Normale Supérieure, Paris, France

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing

Contributed equally with
Yuri B Lebedev and Thierry Mora

For correspondence
awalczak@lpt.ens.fr

Competing interests
Reviewing editor, eLife

"This ORCID iD identifies the author of this article:" 0000-0002-2686-5702

Funding

Russian Science Foundation (15-15-00178)

Dmitriy M Chudakov
Ilgar Z Mamedov
Yuri B Lebedev

European Research Council (724208)

Aleksandra M Walczak

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

Acknowledgements

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

Reviewing Editor

Arup K Chakraborty, Massachusetts Institute of Technology, United States

Publication history

Received: October 24, 2017
Accepted: March 12, 2018
Accepted Manuscript published: March 13, 2018 (version 1)
Version of Record published: March 28, 2018 (version 2)

Copyright

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.