2.1 Immprint scores

To discriminate between the autologous and heterologous scenarios, one can simply count the number of unique nucleotide receptor sequences shared between the two samples, which we call $$ \mathcal{S}$$. Samples coming from the same individual should have more receptors in common because T-cells are organized in clones of cells carrying the same TCR. By contrast, $$ \mathcal{S}$$ should be low in pairs of samples from different individuals, in which sharing is due to rare convergent recombinations. Appropriately setting a threshold to jointly minimize the rates of false positives and false negatives (Fig 1B), we can use $$ \mathcal{S}$$ as a classifier to distinguish autologous from heterologous samples. $$ \mathcal{S}$$ is not normalized for sequencing depth and values of $$ \mathcal{S}$$ should not be compared between samples of different size.

The $$ \mathcal{S}$$ score can be improved upon by exploiting the fact that some receptors are much more likely than others to be generated during V(D)J-recombination, with variations in generation probability (P_gen, [13–15]) spanning 15 orders of magnitude. Public sequences (with high P_gen) are likely to be found in multiple individuals [16], while rare sequences (low P_gen) are unlikely to be shared by different individuals, and thus provide strong evidence for the autologous scenario when found in both samples. To account for this information, we define the score:

2.2 Measuring error rates

We tested the classifiers based on the $$ \mathcal{S}$$ and $$ \mathcal{I}$$ scores on TCRβ RepSeq datasets from 656 individuals [17]—labeled according to their cytomegalovirus (CMV) serological status. Sequences were downsampled to mimic experiments where M₁ = M₂ = M cells were sampled and their receptors sequenced. The frequency of a particular clonotype in the sample (the proportion of cells expressing a particular beta chain) was estimated using the read counts of unique TCRβ sequences, and the mean values of $$ \mathcal{S}$$ and $$ \mathcal{I}$$ computed with a procedure designed to correct for the limited diversity of the sampled repertoire relative to the full repertoire, see Methods 4.3. Similar results may be obtained when M₁ and M₂ are different (see Methods). The clones most often shared between two autologous samples are also the most clonally expanded—and hence are probably antigen experienced. We verified that the sequence statistics of those expanded clonotypes did not differ from generic ones (S1 Fig).

In Fig 1C, we plot the mean value of $$ \mathcal{S}$$ (over many draws of M = 5000 cells) for each individual (autologous scenario, in blue) and between pairs of different individuals (heterologous scenario, in yellow). The two scenarios are clearly discernable under both scores. This result holds for pairs of monozygotic twins obtained from a distinct dataset [18] (green dots), consistent with previous reports that twins differ almost as much in their repertoires as unrelated individuals [18–20]. Heterologous scores (yellow dots) vary little, and may be bounded from above by a theoretical prediction (dashed line) based on a model of recombination [21] (see Methods). On the other hand, autologous scores (blue dots) show several orders of magnitude of variability across individuals. These variations stem from the clonal structure of the repertoire, and correlate with measures of diversity (S2 Fig), which is known to vary a lot between individuals and correlates with age [22], serological status, and infectious disease history [23, 24]. To explore the worst case scenario of discriminability, hereafter we will focus on the individual with the lowest autologous $$ \mathcal{S}$$ found in the dataset.

The sampling process introduces an additional source of variability within each individual. Two samples of blood from the same individual do not contain the exact same receptors, and the values of $$ \mathcal{S}$$ and $$ \mathcal{I}$$ is expected to vary between replicates. Examples of these variations are shown in Fig 1D. The blue (respectively yellow) curves correspond to the sample distributions in the autologous (heterologous) scenario for different individuals (pairs of individuals). The distribution of $$ \mathcal{S}$$ is well-approximated by a Poisson distribution, while $$ \mathcal{I}$$ follows approximately a compound distribution of a normal and Poisson distributions (see Methods for details). Armed with these statistical models of variations, we can predict upper bounds for the false negative and false positive rates. As seen from the detection error trade-off (DET) graph Fig 1E, the Immprint classifier performs very well for a few thousand receptors with an advantage for $$ \mathcal{I}$$.

With 10, 000 cells, corresponding to ∼10 μL of blood, Immprint may simultaneously achieve a false positive rate of < 10⁻¹⁶ and false negative rate of < 10⁻⁶, allowing for the near-certain identification of an individual based on the $$ \mathcal{I}$$ score in pairwise comparisons against the world population ∼10¹⁰. When a large reference repertoire has been collected (M₁ = 1, 000, 000, corresponding to ∼1mL of blood), an individual can be identified with just 100 cells (S3 Fig).

The AUROC estimator (Area Under the Curve of the Receiver Operating Characteristic), a typical measure of a binary classifier performance, can be used to score the quality of the classifier with a number between 0.5 (chance) and 1 (perfect classification). The $$ \mathcal{I}$$ score outperforms the $$ \mathcal{S}$$ score (Fig 1F), particularly above moderate sample sizes (M ≈ 5000). Both scores can be readily generalized to the case of paired receptors, if the pairing of the two chains is available through single-cell sequencing [25–27] or computational pairing [28], we can consider the sequence of the full receptor TCRαβ, instead of just TCRβ. The generation probability of the full TCR is then given by P_gen(α, β) = P_gen(α) × P_gen(β)[29]. Because coincidental sharing of both chains is substantially rarer than with the β chain alone, using the paired chain information greatly improves the classifier.

While this paper focuses on T-cells and TCR sequences, the structure of the B-cell receptors (BCR) repertoire is very similar to the TCR repertoire and we expect to find qualitatively similar results. As an example we use the dataset obtained in Ref. [30] to measure $$ \mathcal{S}$$ and $$ \mathcal{I}$$ for IGH chains (forming half of the BCR receptor), S4 Fig. We see that 5000 IgG+ memory B-cells are enough to reliably identify the individuals in the study. However, B-cells are 5 times less common in peripheral blood than T-cells, and somatic hypermutation tends to distort the statistics of receptors, reducing the reliability of our classifier. Hence, for practical applications, T-cells are better means of identification.

2.3 Evolution with time

The previous results used samples obtained at the same time. However, immune repertoires are not static: interaction with pathogens and natural aging modify their composition. The evolution of clonal frequencies will decrease Immprint’s reliability with time, especially if the individual has experienced immune challenges in the meantime.

To study the effect of short-term infections, we analyzed an experiment where 6 individuals were vaccinated with the yellow fever vaccine, which is regarded as a good model of acute infection, and their immune system was monitored regularly through blood draws [18]. We observe an only moderate drop in $$ \mathcal{I}$$ caused by vaccination (Fig 2A).

Fig 2

A) Evolution of $$ \mathcal{I}$$ (M = 5000) during vaccination, between a sample taken at day 0 (vaccination date) and at a later timepoint. Each color represents a different individual. Each pair timepoint/individual has two biological replicates. The dashed line represents the threshold value. B) Evolution of $$ \mathcal{I}$$ between a sample taken at year 0 and a later timepoint. The red histogram corresponds to one of the individuals sampled in [18] and the blue curves show theoretical estimates, fitted to match (τ = 0.66). C) Evolution of the (normalized) mean of $$ \mathcal{I}$$ (M = 5000) as a function of time for different values of the turnover rate τ. The dashed line represents the threshold value divided by the smallest value of $$ \mathcal{I}$$_t = 0 (M = 5000) in the data. The data points were obtained from the datasets [35] (yellow), [18] (green) and [22] (orange). Different markers indicate different individuals.

This is consistent with the fact that infections lead to the strong expansion of only a limited number of clones, while the rest of the immune system stays stable [31–34]. While other types of infections, auto-immune diseases, and cancers may affect Immprint in more substantial ways, our result suggests that it is relatively robust to changes in condition.

We then asked how stable Immprint is over long times. Using longitudinal datasets [18, 22, 35], we show in Fig 2 that the Immprint score $$ \mathcal{I}$$ is only slightly reduced for samples collected at intervals of up to three years. For longer timescale, given the lack of longitudinal datasets, we turn to mathematical models [12, 36–39] to describe the dynamics of the repertoire. Following the model of fluctuating growth rate described in Ref. [38], we define two typical evolutionary timescales for the immune system: τ, the typical turnover rate of T-cell clones, and θ, which represents the typical time for a clonotype to grow or shrink by a factor two as its growth rate fluctuates. The model predicts a power-law distribution for the clone-size distribution, with exponent −1 − τ/2θ. This exponent has been experimentally measured to be ≈ − 2 [10], which leaves us with a single parameter τ, and θ = τ/2. An example of simulated evolution of Immprint with time is shown in Fig 2B The highlighted histogram represents a data point at two years obtained from [18]. While a fit is possible for this specific individual (τ = 0.66 in Fig 2B), the τ parameter is not universal, and we expect it to vary between individuals, especially as a function of age. In Fig 2C we explore the effect on the stability of Immprint for a range of reasonable values for the clone turn-over rate τ, from 6 months to 10 years, encompassing both previous estimates of the parameter [38] and measured turnover rates for different types of T-cells [40]. While Fig 2 focuses on $$ \mathcal{I}$$, the behaviour is similar for $$ \mathcal{S}$$ (see S5 Fig). We observe that under this model, for most individuals and bar exceptional events, Immprint should conserve its accuracy for years or even decades.

4.1 Datasets & pre-processing

We use five independant RepSeq datasets in this study: (i) genomic DNA from Peripheral blood mononuclear cells (PBMCs) from 656 healthy donors [17]; (ii) cDNA of PBMCs sampled from three pairs of twins, before and after a yellow-fever vaccination [18]; (iii), (iv) two longitudinal studies of healthy adults [22, 35];(v) cDNA dataset of IGH genes (B-cells) from 9 individuals, with multiple replicates [30].

CDR3 nucleotide sequences were extracted with MIGEC [45] (for the second dataset) coupled with MiXCR [46]. We also extract the frequency of reads from the three datasets. The non-productive sequences were discarded (out-of-frame, non-functional V gene, or presence of a stop codon). The generation probability (P_gen) was computed using OLGA [47], with the default TCRβ model. The frequency of each clone was estimated by summing the frequencies of all reads that shared the same nucleotide CDR3 sequence and identical V,J genes.

The preprocessing code is distributed on the Git repository associated with the paper. We also developed a command-line tool (https://github.com/statbiophys/immprint) that discriminates between sample origins, and a companion webapp (https://immprint.herokuapp.com).

4.2 Discrimination scores

To discriminate between the autologous and heterologous scenarios, we introduce a log-likelihood ratio test between the two possibilities:

Since only the term y₁(s) = y₂(s) = 1 (shared sequences) is different between the autologous and heterologous cases, we obtain:

Further assuming N_c p(s) ≪ 1, and p(s) = P_gen(s)q⁻¹ (where q accounts for selection [21] and P_gen(s) is the probability of sequence generation [14]), the score simplifies to Eq 1, with γ = −ln(qN_c〈f²〉) = ln(q⁻¹〈f〉/〈f²〉). The factor γ depends on unknown parameters of the model, but can be estimated assuming a power-law for the clone size distribution [48], ρ(f)∝ f⁻² extending from f = 10⁻¹¹ to f = 0.01, and q = 0.01 [21], yielding γ ≈ 12.24. Alternatively we optimized γ to minimize the AUROC, yielding γ ≈ 15 (S9 Fig). Since performance degrades quickly for larger values, we conservatively set γ = 12.

4.3 Estimating mean scores from RepSeq datasets

The sampling of M cells from blood is simulated using large repertoire datasets. In a bulk repertoire sequencing dataset, the absolute number of cells for each clonotype (cells with a specific receptor) is unknown, but the fraction of each clonotype can be estimated using the proportion of reads that are associated with this specific receptor. To estimate the autologous $$ \mathcal{S}$$ and $$ \mathcal{I}$$ of two samples of size M₁ and M₂ in the absence of true replicates, we computed their expected values from a single dataset containing N reads, from which two random subsamples of sizes M₁ and M₂ were taken. The mean value of $$ \mathcal{S}$$ is equal to $$ \langle \mathcal{S}\rangle ={\sum }_s(1-{(1-f\left(s\right))}^{M_1})(1-{(1-f\left(s\right))}^{M_2})$$, where f(s) is the true (and unknown) frequency of sequence s. A naive estimate of 〈mS〉 may be obtained by repeatly resampling subsets of sizes M₁ and M₂ from the observed repertoire, calculate $$ \mathcal{S}$$ for each draw, and average. One get the same result by replacing f(s) by $$ {\widehat{f}}_s=n\left(s\right)/N$$ in the previous formula, where n(s) is the number of s reads in the full dataset, and N = ∑_s n(s). However, this naive estimate leads to a systematic overestimate of the sharing (visible when compared with biological replicates, see S10 Fig), simply because this procedure overestimates the probability of resampling rare sequences, in particular singletons whose true frequency may be much lower that 1/N. A similar bias occurs when computing $$ \mathcal{I}$$. To correct for this bias, we look for a function h(n) that satisfies for all f:

Expanding the right-hand side of Eq 9 into 4 terms, we find that $$ h\left(n\right)=1-g_{M_1}\left(n\right)-g_{M_2}\left(n\right)+g_{M_1+M_2}\left(n\right)$$ satisfies Eq 9 provided that:

Under the change of variable x = f/(1−f), the expression becomes:

Identifying the polynomial coefficients in xⁿ on both sides yields:

These corrected estimates agree with the direct estimates using biological replicates (S10 Fig).

Similarly, $$ \langle \mathcal{S}\rangle $$ and $$ \langle \mathcal{I}\rangle $$ in heterologous samples can be estimated using:

4.4 Theoretical upper bound on heterologous scores

When the two samples were extracted from two different individuals (heterologous scenario), we can use the universality of the recombination process to give upper bounds on the values of $$ \mathcal{S}$$ and $$ \mathcal{I}$$. These bounds are represented by the dashed lines in Fig 1C). If two samples of respectively M₁ and M₂ cells, containing T₁ ≤ M₁ and T₂ ≤ M₂ unique sequences are extracted from two different individuals, the number of shared sequences between them is given by [21]:

p(s) is the probability of finding a sequence s in the blood. Following [21], we make the approximation p(s) = P_gen(s)q⁻¹, where the q = 0.01 factor is the probability that a generated sequence passes selection. Then 〈p(s)〉 can be estimated from the mean over generated sequences. Similarly, we can estimate $$ \mathcal{I}$$ as

4.5 Error rate estimates

To make the quantitative predictions shown in Fig 1, we need to constrain the tail behavior of the distributions of $$ \mathcal{S}$$ and $$ \mathcal{I}$$, for the two scenarios.

The $$ \mathcal{S}$$ statistic can be rewritten as a sum of Bernouilli variables over all possible sequences, each with a parameter corresponding to its probability of being present in both samples, either in the autologous or the heterologous case. Therefore $$ \mathcal{S}$$ is a Poisson binomial distribution, a sum of independent Bernouilli variables with potentially different parameters. The variance and tails of that distribution are bounded by those of the Poisson distribution with the same mean, denoted by m_a for the autologous case, and m_h for the heterologous case (S11 Fig).

Thanks to that inequality, the rates of false negatives and false positives for a given threshold r are bounded by:

Similarly, $$ \mathcal{I}$$ can be viewed as a sum of $$ \mathcal{S}$$ independent random variables, all following the distribution of ln(1/P_gen)−γ. However, this distribution differs in the two scenarios. Sequences shared between more than one donor have an higher probability of being generated, their ln(P_gen) distribution has higher mean and smaller variance (S12 Fig).

The sum is composed of a relatively large number of variables in most realistic scenarios. Hence, we rely on the central limit theorem to approximate it by a normal distribution, of mean and variance proportional to $$ \mathcal{S}$$. Explicitly:

The AUROC are computed based on these estimates, by numerically integrating the true positive rate $$ P(\mathcal{S},\mathcal{I}<r|\text{heterologous})$$ with respect to the false negative rate $$ P(\mathcal{S},\mathcal{I}<r|\text{autologous})$$ as the threshold r is varied.

4.6 Modeling the evolution of autologous scores

We use the model of Ref. [38] to describe the dynamics of individual T-cell clone frequencies f under a fluctuating growth rate reflecting the changing state of the environment and the random nature of immune stimuli:

With the change of variable x = ln(f), these dynamics simplify to a simple Brownian motion in log-frequency: ∂_t x = −τ⁻¹+ θ^−1/2 η(t). In that equation, τ appears as the decay rate of the frequency, while θ is the timescale of the noise, interpreted as the typical time it takes for the frequency to rise or fall by a logarithmic unit owing to fluctuations. Considering a large population of clones, each with their independent frequency evolving according to Eq 24, and a source term at small f corresponding to thymic exports, one can show that the steady-state probability density function of f follows a power-law [38], ρ(f) ∝ f^−α, with exponent α = 1+ 2θ/τ. α was empirically found to be ≈2 in a wide variety of immune repertoires [10, 48–50], implying 2θ ≈ τ. The turn-over time τ is unknown, and was varied from 0.5 years to 10 years in the simulations.

We simulated the evolution of human TRB repertoires by starting with the empirical values of the frequencies of each observed clones, $$ f(s,0)=\widehat{f}(s,0)=n(s,0)/N$$ from the analysed datasets. A sample of size M was randomly selected with respect to these frequencies, and the frequencies of the clones captured in that sample were then evolved with a time-step of 2 days using Euler-Maruyama’s method, which is exact in the case of Brownian motion. Clones with frequencies falling below 10⁻¹¹ (corresponding to a single cell in the organism) were removed. At each time t > 0, we measured the mean value of $$ \mathcal{S}$$ with the formula ∑_s(1 − (1 − f(s, t))^M) where the sum runs over the sequences captured in the initial sample.