Estimating T-cell repertoire diversity: limitations of classical estimators and a new approach

1The human T-cell receptor (TCR) repertoire—the range of different TCRs expressed—plays a vital role in host defence. By recombination, random insertion, deletion and substitution, the small set of genes that encode the T-cell receptor has the potential to create between 1015 and 1020 TCR clonotypes (a clonotype is a population of T cells that carry an identical TCR). However, the actual diversity of a person's TCR repertoire cannot possibly lie in this range. There are only an estimated 1013 cells in the human body, and many clonotypes are of high abundance due to strong selection forces (for example, thymic education or antigen specificity). The actual, or realized, diversity of the human TCR repertoire remains unknown. The term ‘diversity’ is commonly used to mean either the number of classes (also known as ‘species richness'), or the degree of dispersion among those classes. In this study, we use the term ‘species' to refer to a single TCR clonotype, and ‘diversity’ to refer to the number of TCR clonotypes.

2TCRs are heterodimers and fall into two classes: TCR-αβ and TCR-γδ; γδ T cells constitute 1–10% of the T-cell repertoire. A variable (V), joining (J) and constant region (C) constitute the TCR α- and γ-chains. The TCR β- and δ-chains are also made up of a V, J and C region, with an additional diversity (D) region. One segment from each region is recombined, with additional nucleotide additions and/or deletions, to generate each rearranged TCR (figure 1). This recombination generates high T-cell diversity and enables the recognition of millions of antigens. While V(D)J gene rearrangement is believed to be random, some clonotypes are produced more commonly than others, leading to unequal frequencies of naive T-cell clonotypes and to ‘public’ clonotypes, i.e. clonotypes shared between people. This unequal frequency distribution is believed to be due to a process known as convergent recombination, whereby certain nucleotide sequences can be produced using a greater variety of recombination events; certain amino acid sequences can be made by a greater number of nucleotide triplets; and certain TCRs require fewer nucleotide insertions, deletions or substitutions.

3The third complementarity-determining region 3 (CDR3) of both the TCR α- and TCR β-chains straddles the V(D)J junction (figure 1b), the primary site of antigen contact. The CDR3 is the region most affected by recombination, and the CDR3 region of the β-chain accounts for most of the variation within a person's T-cell repertoire. Antigenic cross-reactivity of T cells results in a discrepancy between structural diversity (the number of different nucleotide or amino acid TCR combinations in the host) and functional diversity (the number of different antigens recognized by the T-cell repertoire).

4TCR diversity is associated with the effective control of viral infections and other pathogens. The number of clonotypes observed in the blood in one person has been reported to decrease with age, viral challenge, immunization and as a result of immune suppression, for example after haematopoietic stem cell transplantation (HSCT). TCR diversity has also been positively associated with autoimmunity in both mice and humans. Accurate quantification of diversity is important to assess the extent of immune convergence (sharing of clonotypes between people).

5Species diversity is also important in many systems outside T-cell immunology, for example, in estimating the repertoire of antibody classes, assessing the size of the metagenome in microbial communities and measuring the rate of evolution of quasi-species of a pathogenic virus. The original motivation for estimating diversity comes from population ecology, where the question of how many species there are in a given population gives rise to the ‘unseen species problem’: how many species are present, but unobserved, in the population of interest? Typically, there is a nonlinear relationship between the number of individuals (e.g. a T cell, a microbe) and the number of ‘species' (e.g. a clonotype or viral variant), and so diversity cannot usually be estimated through linear scaling.

6Estimating the diversity of the T-cell repertoire is difficult for many reasons. First, the repertoire is highly diverse. Given the number of T cells, (assumed to be of the order of 1012), a diversity of (say) 107 clonotypes is unlikely to be directly observed owing to the limited volume of blood that can be taken from a person at any one time, and to the heavy-tailed frequency distributions with highly non-uniform clonotype abundances.

7Second, the precise relationship between the diversity of different TCR-α and TCR-β sequences and the actual TCR functional diversity is unclear. Most recent studies focus on the CDR3 region, because it is the most variable region and because it is short enough to be captured in a single sequence read. However, a T-cell receptor consists of pairings between either α and β chains or γ and δ chains; this adds a further level of diversity that is not routinely captured by many sequencing approaches. Furthermore, the relationship between TCR sequence and three-dimensional structural diversity and functional diversity are not fully understood.

8Third, laboratory techniques that give absolute and unbiased estimates of clonotype frequency are technically challenging. Early studies measured TCR diversity qualitatively, where different clonotypes were identified visually as discrete bands on genomic southern blots. Other approaches used flow cytometry to measure the average observed frequency of each clonotype, reasoning that if this frequency was low then the population was more diverse.

9Greater precision was achieved with spectratyping, where the number of different CDR3 lengths is used as a proxy for the number of clonotypes. The degree to which the frequency distribution of CDR3 lengths deviates from normality is used as a metric of clonal expansion and thus of reduced diversity (because of limited lymphocyte capacity). Although this inference seems reasonable, expansion of some clonotypes does not imply the extinction of other clonotypes, merely their reduced relative frequency. Spectratyping produces incomplete sequence information without further subcloning of the PCR product which is low-throughput and labour-intensive.

10High-throughput sequencing (HTS) allows greater sequencing depth and significantly more accurate quantification of TCR clonotype abundance, albeit at a greater expense than spectratyping. However, HTS is still subject to PCR bias and sequencing error, with the consequences that clonotype abundances can be drastically distorted and that non-existent clonotypes can be recorded, thus falsely increasing the observed diversity.

115′ rapid amplification of cDNA ends (RACE) is reported to suffer from markedly less bias than other HTS approaches. Nevertheless, 5′ RACE (and unbiased sequencing more generally) is unlikely to be sufficient for diversity estimation. Diversity estimation usually makes use of two quantities: the relative abundances of observed species, and the extent to which each species is repeatedly observed in the sample. If PCR amplification is unbiased, then relative abundances will be preserved but the degree of repetition in the sample will not.

12Because all T cells within a sample of blood will not usually be detected in a single sequencing experiment, many researchers have used ‘exhaustive sequencing’, i.e. the library is sequenced with the greatest possible depth, to maximize the number of reads per clonotype. It can then be justifiably concluded that further sequencing of the same library would not yield greater observed diversity. It is therefore tempting to conclude that the sample of blood contains a complete census of clonotypes in the periphery. However, such a conclusion would be false.

13The principle that exhaustive sequencing does not capture full repertoire diversity was demonstrated by Warren et al.. The authors exhaustively sequenced a library derived from a peripheral blood sample. However, upon sequencing a second library derived from the same blood sample, they found that 75% of the sequences returned were new, i.e. not contained in the first library. Furthermore, sequencing data were obtained from a second independent blood sample, and only 13% of the clonotypes observed in the second sample were observed in the first. This indicates that exhaustive sequencing of a single sample is incapable of capturing diversity, regardless of the apparent degree of repetition of species provided. That is, a saturating relationship between the number of reads and the number of clonotypes does not imply that there is a saturating relationship between the number of T cells and the number of clonotypes. The limiting factor is the number of TCRs present in the sample, not the extent of amplification or depth of sequencing. We have observed similar effects of ‘false repetition’ and ‘false saturation’ in our work: figure 2 shows apparent saturation of the number of new clonotypes observed as the number of sequence reads increases. However, the number of clonotypes in the full data (2.5 × 104) is a drastic underestimate of TCR diversity, where between 105 and 106 distinct CDR3 sequences have been directly observed. Finally, it has been noted that exhaustive sequencing of either or both of the TCR α and β chains is insufficient to capture the full repertoire of a person.

14Recent approaches have used DNA barcoding or amplicon length discrimination to resolve the problems of PCR amplification. Under DNA barcoding, a clonotype is identified by its nucleotide or amino acid sequence, but a second identifier is assigned to each individual short DNA sequence through the addition of a random DNA sequence label. Thus, the combination of a given clonotype nucleotide sequence and a given random label is unique. This allows identical T cells to be distinguished from identical sequence reads, and so preferential amplification is irrelevant. For example, if there are two amplicons that have identical CDR3 sequences and identical labels, it can be concluded that both amplicons have been derived from a single DNA sequence. The resulting data therefore consist of absolute—not relative—clonotype abundances, which are required for any abundance-based estimator. Furthermore, DNA barcoding can be extended to correct for sequencing error.

15Another factor that prevents absolute quantification of TCR abundance is the sequencing of cDNA rather than genomic DNA, since a single T cell may express multiple mRNA copies. Therefore, cDNA is not suitable for diversity estimation.

16Even where data collection involves considerable sequencing depth, and where unbiased data have been obtained, estimators of the number of unseen clonotypes will need to be employed because of limits on blood volume that can be taken from donors. Several estimators of species richness (i.e. the number of species) developed in ecology have been applied to estimate TCR diversity, treating each clonotype as a ‘species'. Such estimators fall into two broad categories: parametric estimators, where the shape of the species frequency distribution is assumed to follow some analytical form, and non-parametric estimators that make no such assumptions, and thus population frequencies cannot be inferred. Since the true numbers of species or clonotypes are unknown, it is difficult to validate estimators of diversity, and so in common with ecological populations, it is often unclear which estimator should be used. We compare diversity estimators below and in table 1.

17There are important qualitative differences between T-cell repertoires and ecological populations in the uniformity of sampling. In ecological populations, data collection is frequently not random. For example, while placement of quadrats may be random, all of the individuals present in that quadrat are counted, leading to clustering of data. Also, the probability of detection varies between species as it is influenced by colour, physical size, noise emission, geographical distribution, movement, variety of habitats and relationship to other species. By contrast, in samples of T cells derived from blood, it is reasonable to assume that individual T cells have the same probability of detection; this assumption is less justifiable in solid tissue, as for example, lesions are non-randomly sampled.

18In many ecological populations (e.g. plants, arthropods), the actual counting of individuals present in the sample is more straightforward than for populations of T cells, where sequencing introduces biases and where it is difficult to distinguish sequencing errors from rare species. The frequent implicit assumption that sequencing data comprised individuals that are equally detectable is often inappropriate. The probability that a given sequence read is recorded is conditional on two events: first, the probability that the T cell is sampled from blood, which is equal among T cells; and second, the probability that an amplicon from a T cell is amplified, which is not equal across all CDR3 sequences. This problem does not arise in ecological datasets.

19The use of diversity indices developed in ecology that are used in T-cell repertoires is not restricted to species richness estimators. Similarity indices such as the Jaccard, Morisita-Horn, analysis of similarity (ANOSIM) and dispersion metrics such as Simpson's diversity index, the Shannon entropy and Renyi entropy have been used to compare the TCR diversities between different people or between different T-cell phenotypes. Many of the difficulties that arise in applying ecological species richness estimators to T-cell repertoires also confound the measurement of the extent of dispersion or similarity between repertoires, and ecological indices should be used with caution when analysing TCR repertoires.

20One of the most commonly used estimators is Chao1 or its bias-corrected form (Chao1-bc). These estimators have been used to estimate TCR diversity in mice, and humans, making use of an amendment to the estimator that takes account of the maximum upper bound of diversity.

21The abundance-based coverage estimator (ACE), which has been suggested as best practice and is commonly used in ecology, has been used to estimate repertoire diversity in transgenic mice in the contexts of T-cell differentiation, and TCR specificity and self-recognition. However, Hsieh et al. note that ACE is based on the probability that uniform sampling would produce the observed frequency distribution.

22Weinstein et al. used a capture–recapture approach to estimate the size of the antibody repertoire in zebrafish, and this approach was extended in Glanville et al. to estimate antibody diversity in humans. The latter study also used a technique that allows sequencing of reads long enough to span all three CDR regions, which would allow more direct data on T-cell repertoires to be collected. No validation of the capture–recapture method was performed in either study.

23Non-parametric abundance-based species richness estimators have been validated using ecological populations that have been extensively sampled and where approximate species richness is assumed to be known. However, the accuracy of these ecological estimators has been questioned in immunological populations. We recently compared the performance of widely used non-parametric species richness estimators (the Chao1bc, ACE, Bootstrap and Good-Turing estimators) from population ecology when applied to immunological and microbiological systems. We considered three distinct sets of data: the clonal distribution of cells naturally infected with human T-lymphotropic virus type-1 (HTLV-1), operational taxonomic units (OTUs) of Bifidobacteria in the gastrointestinal tract of infants, and T-cell receptor repertoires. In the case of HTLV-1, a ‘species' is a clone, defined as a population of infected cells that share a genomic site of proviral integration.

24For each set of data, we found that all estimators were biased by sample size (figures 3 and 4). This is problematic as estimates of species richness would increase if, for example, greater blood volumes were drawn or technique sensitivity was improved. Furthermore, there was strong evidence that the estimators underestimated diversity. Firstly, the estimators frequently produced estimates from subsamples that were lower than the diversity of the full observed sample. Secondly, in almost all cases, only a small number of unseen ‘species' was predicted in addition to those observed. Such estimates are implausible in the HTLV-1 and T-cell repertoire datasets where there is such a vast potential diversity.

25The incidence-based coverage estimator was used to estimate the diversity of regulatory T cells in transgenic mice, although no validation of this estimator was performed. Qi et al. used multiple replicate libraries to compute the Chao2 estimator, resulting in an estimate of 108 clonotypes. The Chao2 estimator makes use of incidence (i.e. presence or absence) data, as opposed to abundance data, across different replicates. The estimator therefore provides a method of avoiding the distorted abundances due to PCR amplification mentioned above.

26To validate their approach, the authors created an in silico Zipf distribution of clonotype frequencies from which to sample. They took random samples of varying sizes and found that the estimated diversity accurately estimated the number of clonotypes in their in silico distribution. Although indirect and using only one in silico distribution, this validation suggests that their method holds promise. However, we have applied Chao2 to HTLV-1 and microbial OTU data, and we again observed a bias with sample size, as seen with the other non-parametric estimators we tested (figures 3 and 4).

27Robins et al. frequently observed as many as approximately 105 TCR clonotypes using single-molecule DNA sequencing, and employed a method originally devised by Efron & Thisted and amended in Ionita-Laza et al. to estimate a peripheral blood diversity of 3 to 4 × 106 clonotypes, including 1 × 106 antigen-experienced T-cell clonotypes, where the latter is approximately one order of magnitude higher than estimated previously. Their method assumes that individual T-cell clonotypes enter the sample according to a Poisson process with clonotype-specific rates, which are inferred from the observed clonotype abundances. The method predicts the number of new sequences that would be observed in a subsequent sample. Hence their method does not merely provide an estimate of TCR diversity, but also the relationship between sample size and diversity. Therefore, the authors were able to validate their method. While this validation was direct, in that observation was compared with the predicted number of additional clonotypes, it was limited to only a single additional sample.

28Several recent studies have made use of the class of Poisson abundance models (PAMs). Sepúlveda et al. noted that species frequency data come from a multivariate hypergeometric distribution (i.e. a multinomial distribution where samples are taken without replacement). Because the size of a sample is dwarfed by the size of the total population (and therefore sampling does not drastically alter clonotype relative abundances), these authors approximated the multivariate hypergeometric distribution using a Poisson distribution. Incorporation of a varying sampling rate for clonotypes of varying frequencies leads to the class of PAMs. They applied their method to previously published data on mice with different phenotypes, and evaluated the consistency of their method by excluding clonotypes above successively higher cut-off frequencies. Worryingly, there was wide variation in diversity estimates across all phenotypes depending on the specific PAM used. Rempala et al. focused on one such model, the bivariate Poisson-lognormal distribution, and concluded that under-sampling in their repertoire datasets is more severe (and thus the population is more diverse) than would be estimated using the Good-Turing estimator. Other extensions of the class of PAMs have been developed that estimate the similarity between populations in the presence of unseen clonotypes.

29Using a compound Poisson process model used originally to estimate gene capture diversity, Wang et al. estimated TCR diversity in the context of T-cell fate and differentiation. This method can also model the relationship between the number of clonotypes and the number of T cells. They estimated approximately 106 unique TCRα and TCRβ CDR3 nucleotide sequences, approximately one third of that predicted by Arstila et al.. Their method was validated previously in the context of gene capture diversity using in silico distributions, choosing lognormal, exponential and gamma distributions of varying diversities. It is unclear how this validation translates to T-cell immunology.

30In addition to the capture–recapture approach used in Weinstein et al., Klarenbeek et al. fitted multiple Poisson mixture models to HTS data to estimate β-chain diversity in the CD4+ and CD8+ T-cell compartments. Extending the distribution fitted to the observed data to model the number of unseen clonotypes, the authors estimated that the memory compartment consists mainly of unexpanded clones and is far more diverse than thought previously (only 2 and 3–10 times less diverse than the naive repertoire in CD4+ and CD8+ T-cell compartments, respectively). Their estimates are also remarkable in that they predict that more than 90% of memory clonotypes are relatively small.

31Power laws have been used to model the form of the T-cell repertoire. An advantage of this method is that the fitted parameters are relatively easy to interpret. It can be shown that one parameter quantifies the proportion of the repertoire occupied by clonotypes of a single T cell, and the other provides a measure of dispersion. Power law characterizations of the T-cell repertoire could be extended to estimate the number of unseen clonotypes in a similar manner to Klarenbeek et al. by extending the modelled distribution.

32Parametric approaches are often evaluated using goodness of fit to the observed data, for example using χ2-tests or Akaike's information criterion (AICc). While these methods are useful for comparative purposes, they do not validate the resulting model's accuracy. A major limitation of all parametric approaches is that the estimated diversity is dependent on the assumed form of the clonotype distribution.

33We developed an estimator named DivE which uses rarefaction curves (figure 5). Similar to a species accumulation curve, an individual-based rarefaction curve is created by cumulating the number of species as the number of observed individuals (e.g. a T cell) increases, in a single resample. Species counts are averaged over multiple resamples of the data to obtain the expected number of species as a function of the number of individuals. Sample-based rarefaction curves plot the expected number of species against the number of samples.

34DivE involves fitting multiple simple mathematical models, many of which are well known in ecological studies, to rarefaction curves, and to nested subsamples of these curves. Novel criteria are then used to determine the most appropriate model; as well as assessing the quality of fit to seen data these criteria also assess the quality of fit to unseen data, i.e. how well a given model can predict the full dataset from random subsets thereof. The best-performing models are then aggregated and extrapolated to a user-specified population size to produce the diversity estimate (figure 5).

35We used three methods to validate the performance of DivE. We measured the extent to which DivE could: (i) estimate the diversity of the observed dataset from subsamples; (ii) estimate from a single dataset the diversity of additional independent HTLV-1 data, obtained using separate blood samples taken in immediate succession and (iii) provide consistent estimates given samples of unequal size. In each validation, the estimator performed better than the non-parametric abundance-based estimators we tested (figure 3). We believe the principal reason that DivE performs well is that candidate models are selected on their ability to consistently predict additional rarefaction data. The additional data (i.e. the full rarefaction curve) have no influence on fitted parameter values, and so DivE not only assesses goodness of fit but also evaluates the accuracy of the model. DivE has been provided as an R package, available at http://cran.r-project.org/web/packages/DivE/index.html.

36Accurate extrapolation of rarefaction curves assumes that the sampled population is representative of the whole population to be extrapolated to. This is a reasonable assumption in the case of T-cell sampling in the blood, i.e. T cells sampled in one blood draw are likely to be representative of all T cells in the peripheral blood. However, this is a poor assumption when trying to infer the TCR diversity in the whole body as T cells sampled in the blood may not be representative of T cells in lymphoid tissue, etc. The difficulty of inferring total population diversity from estimates in the blood is not unique to DivE and will adversely affect the accuracy of all estimators.

37An alternative approach to rarefaction curve extrapolation that is based on a rigorous statistical footing has recently been developed. However, to estimate the rarefaction curve, this method requires an input of species richness (usually provided by ACE or Chao1), which is the quantity we seek to estimate. Furthermore, the authors of these papers caution that this method is not suitable for extrapolation beyond two- or threefold.

38We summarize the advantages and disadvantages of different diversity estimation methods in T-cell repertoire analysis in table 1.

39To estimate repertoire diversity it is essential to obtain unbiased data, with absolute counts of TCR clonotypes. If unbiased absolute count data are not available, neither relative abundances nor the degree of repetition of observations are credible, and so diversity estimators should not be applied. While Chao2 does not require abundance data, we have found that this estimator too is biased by sample size in immunological and microbiological data (figure 3).

40We also caution against estimating diversity using severely under-sampled data, whether due to limited sequencing depth or low blood volume. To quantify ‘under-sampling’, we previously defined a parameter based on the curvature of the observed rarefaction curve (see for further details). A linear rarefaction curve implies an implausible constant rate of species accumulation. As sampling depth increases, the rate of species accumulation should decrease as previously encountered species are repeatedly observed. Abundance-based estimators should not be applied when the rarefaction curve is close to linear.

41Recent advances in HTS combined with DNA barcoding mean that unbiased absolute count data is now increasingly available. However, because of the enormous potential diversity of the TCR repertoire and the limited amount of blood that can be drawn from a donor at any given time, there will almost certainly be unseen TCR clonotypes regardless of the precision of data collection. Therefore, estimators of diversity must be employed. Existing parametric estimators suffer from the requirement of an a priori form of the species frequency distribution. Furthermore, each non-parametric estimator we have tested, either abundance- or incidence-based, was significantly biased by sample size.

42Absolute count data allow important simplifying assumptions to be made about the relationship between the observed data and the underlying T-cell repertoire, namely that individual T cells have been sampled independently, randomly and with equal detection probabilities. These assumptions in turn allow the extrapolation of models fitted to individual-based rarefaction curves. The question of which model to fit, however, is non-trivial. DivE selects which models are most appropriate based on their ability to faithfully reproduce all observed rarefaction data from subsamples, providing a degree of robustness that we have not observed with classical non-parametric estimators. Crucially, the form of the model chosen depends on the data, and so DivE does not require a priori assumptions regarding the form of the clonotype frequency distribution, or regarding the relationship between the number of T cells and the number of TCR clonotypes.

43We have validated DivE across three independent immunological and microbiological systems. In all systems, the estimator was accurate, and considerably more so than the non-parametric estimators we examined. We believe that this estimator will be an important tool to estimate T-cell repertoire diversity.