# Cdr3 Spectratyping Analysis Essay

By In 1

## Abstract

Motivation: The effectiveness of vertebrate adaptive immunity depends crucially on the establishment and maintenance of extreme diversity in the antigen receptor repertoire. Spectratype analysis is a method used in clinical and basic immunological settings in which antigen receptor length diversity is assessed as a surrogate for functional diversity. The purpose of this paper is to describe the systematic derivation and application of statistical methods for the analysis of spectratype data.

Results: The basic probability model used for spectratype analysis is the multinomial model with n, the total number of counts, indeterminate. We derive the appropriate statistics and statistical procedures for testing hypotheses regarding differences in antigen receptor distributions and variable repertoire diversity in different treatment groups.

We then apply these methods to spectratype data obtained from several healthy donors to examine the differences between normal CD4+ and CD8+ T cell repertoires, and to data from a thymus transplant patient to examine the development of repertoire diversity following the transplant.

Availability:http://www.duke.edu/~kepler/spa.html

Contact:kepler@duke.edu

## 1 INTRODUCTION

### 1.1 Biology of antigen receptors

The immune system of all jawed vertebrates has two major divisions, denoted as the innate and adaptive systems, which are distinguished from each other by the absence or presence, respectively, of randomly rearranged antigen receptors. Pathogenic microorganisms display enormous molecular variability compared with their vertebrate hosts, owing to their short generation times and large populations. The strategy of the adaptive immune system, therefore, is to counter this extraordinary variability by using somatic diversification, which occurs at time scales typical of individual somatic cell turnover, rather than the much longer time scales of host germline turnover, or generations. These antigen receptors are the T cell receptors (TCR) and immunoglobulins (Ig), borne by T and B cells, respectively. For an excellent introduction to the biology of the immune system, see Janeway et al. (2005).

Both T cell and B cell receptors are encoded by gene segments that must be rearranged on the chromosome to produce a complete productive gene. This process involves the stochastic selection of gene segments from each of two or three libraries (depending on the type of chain in question) and the further stochastic selection of specific recombination points given these segments, as well as diversity derived from additional non-templated nucleotides. The region containing these highly variable junctions is the third of three complementarity-determining regions (CDRs) that are seen crystallographically to contact antigen. The sources of TCR diversity are naturally broken down hierarchically into gene segment family (library), segment within family, CDR3 length and CDR3 nucleotide diversity.

The immune response to infection involves the manifold expansion of a small number of T cell clones, and consequently, a decrease in the TCR repertoire diversity. Severe immune disruptions, either caused by genetic lesions disrupting normal developmental processes or acquired during life as a result of lymphoma or leukemia, result in loss of antigen receptor diversity or of one or more antigen receptors altogether.

Transplantation of bone marrow or of primary lymphoid organs is performed to establish or re-establish the missing populations of T or B cells. These interventions are often followed up by spectratype analyses to monitor the progress of the diversification of the affected cellular populations. Spectratype analysis (Cochet et al., 1992; Pannetier et al., 1993; Pannetier et al., 1997) provides information on antigen receptor diversity at the level of CDR3 length. The point is not that length diversity itself is of particular relevance (though it might be), but that length heterogeneity is very likely to be representative of overall sequence heterogeneity.

Spectratyping has proven to be a valuable method for the monitoring of antigen receptor repertoire diversity subsequent to lymphoid transplantation and to infection (Cochet et al., 1992; Pannetier et al., 1993; Pannetier et al., 1997; Bousso et al., 2000; Sarzotti et al., 2003; Markert et al., 2003). The purpose of this paper is to present a set of statistical methods that bring further power and flexibility to the use of this assay. Previous work on bringing statistical methods to spectratype analysis include Collette and Six (2002) and Collette et al. (2003) who have developed statistical methods and implemented them in an Excel package, and Gorochov et al. (1998). These methods, while clearly valuable, are based on ad hoc measures of spectratype differences. Our aim in the present paper, in contrast, is to start from first principles and derive the statistical measures and tests most natural for the assay itself.

The methods developed in this paper can be applied to both TCR and Ig; for the sake of clarity and simplicity, we will specifically describe the procedures and background in terms of TCR.

The paper is structured as follows. First, we describe the molecular biological assay itself and present examples of the data generated. Next, we derive the probability density functions (pdfs) and compute their relevant characteristics. We next develop two fundamental statistical procedures and illustrate their use with data from healthy human subjects and from human thymus transplant patients. Finally, we discuss the mathematical relationship between the TCR repertoire diversity and the primary statistic, d, derived in our analysis.

### 1.2 Spectratype analysis

The spectratype assay begins with the collection of a peripheral blood sample from the subject, and the isolation of CD4 or CD8 T cells from it, although total CD3 or PBMC can be used. PCR is used to specifically replicate the variable-length region (CDR3) of the rearranged TCR variable region beta chain (TCRBV). Primers specific to individual TCRBV families (or family subsets) are used to provide independent spectratypes for each. The resulting mixture of CDR3 replicons is size-separated by electrophoresis, and quantified by densitometry (Fig. 1).

Spectratypes are presented as histograms of the number of T cells bearing receptors versus receptor length for each of the 24 TCR V beta families tested. While the utility of these histograms for the monitoring of TCR diversity is widely accepted, the theoretical relationship between these histograms and the underlying receptor diversity has not been described. Spectratype histograms have largely been analyzed by subjective classification. We have undertaken this study to allow the quantitative and objective analysis of spectratype data, and to render the relationship between spectratype and receptor diversity unambiguous.

## 2 METHODS

### 2.1 Experimental methods

All of the data discussed in this paper were collected from human subjects as follows.

The immunoscope analysis of PCR-amplified products is performed as described by Pannetier et al. (1997). Briefly, RNA is extracted from PBMC samples (210 × 106 cells/sample) using Triazol (Life Technology, Gaithersburg, MD) and reverse transcribed to single-stranded cDNA with AMV reverse transcriptase using an oligo (dT) primer according to the manufacturer's protocol (Promega, Madison, WI). The newly synthesized cDNA is then used as a template for 24 PCRs. The PCRs are carried out in 20 μl volumes by standard procedures (Pannetier et al., 1997), using 24 Vbeta primers and one Cbeta 1 primer. The samples are subjected to 40 cycles of denaturation (25 s at 94°C), annealing (45 s at 60°C) and elongation (45 s at 72°C). After the last cycle, a final elongation step (5 min at 72°C) is performed. The PCR products are visualized on a 1.5% agarose gel by ethidium bromide staining before using 2 μl of the amplified products for a run-off elongation reaction with a fluorescent (G-FAM) C beta 2 primer, as described previously (Pannetier et al., 1997; Bousso et al., 2000; Sarzotti et al., 2003; Markert et al., 2003). The elongation products are then run on an ABI 3100 Genetic Analyzer. Fluorescence-labeled size markers (Applied Biosystems, CA), are loaded with the run-off products. After analysis on an automated sequencer (Applied Biosystems) size determination of the run-off products and the analysis of the CDR3 region products are performed using the GeneScan software.

### 2.2 Computational methods

Software for the numerical computations was written in Fortran90. Web-based implementation is described by M. He, J.K. Tomfohr, B.H. Devlin, M. Sarzotti, M.L. Markert and T.B. Kepler (submitted for publication)Additional routine statistical analyses were performed using Splus 6.1 (Insightful, Inc.)

## 3 APPROACH

In this section, we derive the pdfs required to devise the statistical methods for data analysis, identify the statistics that will most directly address the issues of interest in these analyses, and compute the expected values and sampling distributions of these statistics.

### 3.1 Derivation of probability distributions

In order to analyze spectratype data appropriately, we need to determine the relevant probability distributions. The densitometric intensity of the peak at any given electrophoretic displacement is, ideally, proportional to the number of TCR with CDR3 of the corresponding length. If we could simply count these TCR, a multinomial model would be appropriate. In our case, however, we do not measure the absolute number of T cells in the sample. In this subsection, we derive a family of density functions obtained from the multinomial, but in which n, corresponding to the total number of T cells sampled, is uncertain.

#### 3.1.1 Multinomial distribution with n uncertain

To establish notational conventions, suppose m is a L-dimensional multinomial random variable with population parameter vector q. Then, the probability mass function (pmf) for m is given bywhere Cn(m) is the multinomial coefficientfor

$${\sum }_{i=1}^{L}{m}_{i}=n$$

and zero otherwise.

${f}_{M}\left(\hbox{ m }\right|n,\hbox{ q })={C}_{n}(\hbox{ m }){\displaystyle \prod _{i=1}^{L}}{q}_{i}^{{m}_{i}},$

${C}_{n}\left(\hbox{ m }\right)\equiv \frac{\Gamma (n+1)}{{\prod }_{i=1}^{L}\Gamma ({m}_{i}+1)}$

#### 3.1.2 Density functions for relative frequencies

The absolute number mi of T cells with CDR3 length i is not measured; only the relative frequencies rm/n are. For n sufficiently large, the relative frequencies may be treated as continuous random variables with pdf determined by a transformation of variables (Gardiner, 2004)

${f}_{r}\left(\hbox{ r }\right|n,\hbox{ q })={n}^{L-1}{C}_{n}(n\hbox{ r }){\displaystyle \prod _{i=1}^{L}}{q}_{i}^{{r}_{i}n}$

(1)

where nL−1 is the value of the determinant of the Jacobian transformation matrix.
We assume that for each i, nri is large enough for Stirling's approximation (Weisstein, 2004ahttp://mathworld.wolfram.com/StirlingsApproximation.html) to provide an adequate representation of the gamma function, and use it to simplify the multinomial coefficientUsing this expression in Equation (1) we obtainwhere δ is the Dirac delta and we use the dotted index convention: r· ≡ ∑iri, and D(r;q) is the Kullback–Leibler divergence or relative entropy (Kullback and Leibler, 1951),r and q both obey the probability constraint r· = q· = 1.

$\begin{array}{c}log{C}_{n}\left(n\hbox{ r }\right)\equiv log\Gamma (n+1)-{\displaystyle \sum _{i=1}^{L}}log\Gamma ({r}_{i}n+1)\\ ={\displaystyle \sum _{i=1}^{L}}\left\{-n{r}_{i}log{r}_{i}-\frac{1}{2}log{r}_{i}\right\}\\ -\frac{L-1}{2}log2\pi n+O\left(\frac{1}{n}\right).\end{array}$

${f}_{r}\left(\hbox{ r }\right|n,\hbox{ q })=\frac{{n}^{(L-1)/2}{\hbox{ e }}^{-nD(\hbox{ r };\hbox{ q })}}{\sqrt{{\left(2\pi \right)}^{L-1}{\prod }_{i=1}^{L}{r}_{i}}}\delta (r-1),$

$D(\hbox{ r }:\hbox{ q })={\displaystyle \sum _{i=1}^{L}}{r}_{i}log\frac{{r}_{i}}{{q}_{i}}.$

We determine the pdf for D(r;q) by first computing its cumulant generating function (cgf) (Weisstein, 2004b, http://mathworld.wolfram.com/Cumulant-GeneratingFunction.html)

$\begin{array}{c}h\left(s\right)\equiv log\int {\hbox{ d }}^{L}r{\hbox{ e }}^{sD(\hbox{ r };\hbox{ q })}{f}_{r}\left(\hbox{ r }\right|n,\hbox{ q })\\ =\lambda (logn-log2\pi )+log\int {\hbox{ d }}^{L}r\frac{{\hbox{ e }}^{-(n-s)D(\hbox{ r };\hbox{ q })}}{\sqrt{{\prod }_{i=1}^{L}{r}_{i}}}\delta ({r}_{\cdot }-1)\\ =-\lambda log\left(1-\frac{s}{n}\right).\end{array}$

We recognize this cgf as that of a gamma random variable with shape parameter λ ≡ (L − 1)/2 and scale parameter 1/n,

${f}_{D}\left(D\right|n)=\frac{{n}^{\lambda }}{\Gamma \left(\lambda \right)}{D}^{\lambda -1}{e}^{-nD}.$

(2)

When n is uncertain but can be described as a gamma random variable with shape parameter γ and scale parameter 1/ξ, we marginalize Equation (2) over n to get

${f}_{D}\left(D\right|\gamma ,\xi )=\frac{\Gamma (\gamma +\lambda )}{\Gamma \left(\gamma \right)\Gamma \left(\lambda \right)}\frac{{\xi }^{\gamma }{D}^{\lambda -1}}{{(\xi +D)}^{\lambda +\gamma }}.$

We find that the mean and variance are

$E\left[D\right]=\frac{\lambda \xi }{\gamma -1}=\frac{1}{{n}_{0}}[1+{c}_{0}+O({c}_{0}^{2}\left)\right],$

(3)

where E is the expectation operator, andwhere n0 and c0 are the prior mean and coefficient of variation, respectively.

$\hbox{ var }\left(D\right)=\frac{\lambda {\xi }^{2}(2\gamma -3)}{{(\gamma -1)}^{2}(\gamma -2)}=\frac{1}{{n}_{0}^{2}}[2+5{c}_{0}+O({c}_{0}^{2}\left)\right],$

#### 3.1.3 Hierarchical relative multinomial distribution

The distributions we have just derived are useful for pairwise spectratype comparisons. There are situations of great interest for clinical applications, however, where the quantities of interest are the divergences of multiple observed spectratypes from some control or asymptotic spectratype. In these cases, the patient TCR repertoire can be characterized in terms of their divergence from the control repertoire, while the measured spectratype is additionally divergent from true patient spectratype due to sampling variability.

We, therefore, need to consider a hierarchical model in which there are two sampling stages. We start with a relative multinomial with population parameter vector p corresponding to the ideal, perfectly sampled healthy CDR3 length distribution. We imagine that the actual (but not yet observed) CDR3 length distributions in individual subjects are samples q from this ideal distribution. The parameter d−1 describes the completeness of this stage one sampling; the larger d−1, the more similar will be p and q. The observed spectratypes r are samples or size n (stage 2) from these individuals. The pdf for this model is

${f}_{H}\left(\hbox{ r }\right|\hbox{ p },n,{d}^{-1})=\int {\hbox{ d }}^{L}q{f}_{r}\left(\hbox{ r }\right|\hbox{ q },n){f}_{r}(\hbox{ q }|\hbox{ p },{d}^{-1}).$

(4)

We perform the integration by Laplace's method (Erdélyi, 1956). The critical point is given bywhere μ is a Lagrange multiplier enforcing the normalization of q. We solve for μ by summing on i, gettingwhere D0D(q;p). Then

$0=n({r}_{i}-{q}_{i})+{d}^{-1}\left({D}_{0}-log\frac{{q}_{i}}{{p}_{i}}\right)+\frac{1}{2}({L}_{{q}_{i}}-1).$

(5)

$0={q}_{i}\frac{\partial log{f}_{H}}{\partial {q}_{i}}=n{r}_{i}-{d}^{-1}{q}_{i}\left(1+log\frac{{q}_{i}}{{p}_{i}}\right)-\frac{1}{2}-\mu {q}_{i},$

$\mu =n-1+{D}_{o}d-\frac{L}{2}.$

The Hessian matrix, whose determinant is required for the application of Laplace's method, iswhere δij = 1 for i = j and equals zero otherwise.

${H}_{ij}\equiv \frac{{\partial }^{2}log{f}_{H}}{{\partial }_{{q}_{i}}{\partial }_{{q}_{j}}}=-{\delta }_{ij}\frac{1}{{q}_{i}}\left(n\frac{{r}_{i}}{{q}_{i}}+{d}^{-1}-\frac{1}{2{q}_{i}}\right),$

With these results and using appropriate care in applying Laplace's method to integrals with integrands with constraints, we are able to write the pdf,

$\begin{array}{c}log{f}_{H}\left(\hbox{ r }\right|\hbox{ p },n,{d}^{-1})=-n{D}_{1}-{d}^{-1}{D}_{0}-\lambda log\left(\frac{1}{n}+d\right)\\ -\lambda log2\pi -\frac{1}{2}{\displaystyle \sum _{2}^{L}}log{r}_{i},\end{array}$

(6)

where D1D(r;q).

The first hierarchical stage is the one of biomedical interest, and the parameter d, measuring the divergence of the true subject spectratype q from the control spectratype is the quantity we will be trying to estimate. The second stage arises in sampling. The inverse n−1 of the sample size, according to Equation (3), gives the mean sampling divergence. We will hereafter assume adequate sampling, nd ≫ 1. Where the sampling is inadequate, the measurements can provide no more than an approximate lower bound on d.

Two quantities that will prove to be of value in what follows areand

${D}_{*}\equiv \sum {r}_{i}log\frac{{r}_{i}}{{p}_{i}}$

${S}_{D}\equiv {\sum {r}_{i}\left(log\frac{{r}_{i}}{{p}_{i}}-{D}_{*}\right)}^{2}.$

Expanding qi in ε ≡ (dn)−1 using Equation (5) gives

${q}_{i}={r}_{i}+\epsilon {r}_{i}\left({D}_{*}-log\frac{{r}_{i}}{{p}_{i}}\right)+O\left({\epsilon }^{2}\right).$

Substituting this expression into Equation (6) gives the pdf valid under these assumptions,Define the generating functionBy definition, we haveandCombining these results, we find thatFinally, we can write down estimators for d and 1/nand

$\begin{array}{c}log{f}_{H}\left(\hbox{ r }\right|\hbox{ p },n,{d}^{-1})=-m\left({D}_{*}-\frac{{s}_{D}}{2dn}\right)-\lambda log\left(d+\frac{1}{n}\right)\\ -\lambda log2\pi -\frac{1}{2}{\displaystyle \sum _{i=1}^{L}}log{r}_{i}+O\left({\epsilon }^{2}\right).\end{array}$

$\begin{array}{c}{h}_{H}(s,t)\equiv \lambda log\left(\frac{nm}{2\pi (n+{d}^{-1})}\right)\\ +log\int {\hbox{ d }}^{L}r\frac{exp\left\{(s-{d}^{-1})\left({D}_{*}-\frac{{d}^{-1}-s}{2(n-t)}{s}_{D}\right)\right\}}{\sqrt{{\prod }_{i=1}^{L}{r}_{i}}}\\ =\lambda \left[log\left(\frac{1}{{d}^{-1}-s}+\frac{1}{n-t}\right)-log\left(d+\frac{1}{n}\right)\right].\end{array}$

$E\left[{D}_{*}-\frac{{s}_{D}}{dn}\right]=\frac{\partial {h}_{H}}{\partial s}(0,0)=\lambda \left(d-\frac{1}{n}\right)+O\left(d{\epsilon }^{2}\right)$

$E\left[{s}_{D}\right]=\frac{2}{{d}^{2}{n}^{2}}\frac{\partial {h}_{H}}{\partial t}(0,0)=2\lambda \left(d-\frac{1}{n}\right)+O\left(d{\epsilon }^{2}\right).$

$E\left[{D}_{*}\right]=\lambda \left(d+\frac{1}{n}\right).$

$\frac{1}{2\lambda }E\left[{D}_{*}+\frac{1}{2}{s}_{D}\right]=d$

$\frac{1}{2\lambda }E\left[{D}_{*}+\frac{1}{2}{s}_{D}\right]=\frac{1}{n}.$

## 4 DATA ANALYSIS

The pdfs derived above serve as the point of departure for the development of statistical methods for the analysis of spectratype data. In this section, we introduce these tests, and demonstrate their use with spectratype data collected from healthy volunteers, as well as data from thymus transplant patients.

The one-sample test is of limited utility, so in the interest of economy, we start with the tests for the comparison of parameter vectors in two or more treatment groups. Then, we describe the technique of linear modeling for the analysis of variable deviations from a given parameter vector. This latter is based on the hierarchical model.

### 4.1 Comparison of parameter vectors

The first scenario to explore is one in which there are two or more treatment groups, and the issue to be addressed is whether the distribution of CDR3 lengths is the same in all groups. The null hypothesis is that the population parameter vector q is identical in all groups. This will look very familiar, since, as we will see, the Kullback–Leibler divergence supports an additive decomposition of variability analogous to the partitioning of variance in linear models.

Denote by rijk the relative frequency of counts with CDR3 length i in the k-th member of group j, which has nj members. Define the sample mean parameter vector for group j:If the number of groups is G, the grand sample mean isDenote by rjk the vector with components rijk.

${\overline{r}}_{ij.}\equiv \frac{1}{{n}_{i}}{\displaystyle \sum _{k=1}^{{n}_{j}}}{r}_{ijk}.$

${\overline{r}}_{i\cdot \cdot }\equiv \frac{1}{n}{\displaystyle \sum _{j=1}^{G}}{n}_{i}{r}_{ij}$

Then the total divergence is partitioned as

$\begin{array}{c}{D}_{\hbox{ tot }}\equiv {\displaystyle \sum _{jk}}D({\hbox{ r }}_{jk};\hbox{ q })\\ ={\displaystyle \sum _{jk}}\left[D\right({\hbox{ r }}_{jk};{\overline{\hbox{ r }}}_{k\cdot })+D({\overline{\hbox{ r }}}_{k\cdot };\overline{\hbox{ r }}\mathrm{.})+D(\overline{\hbox{ r }}\mathrm{.};\hbox{ q })].\end{array}$

Under the null hypothesis that the groups have identical population parameter vectors, the expected values for each of these partial divergences areandBecause of Equation (2), the statistic

$f\equiv \frac{(n.-G+1){\sum }_{k}D({\overline{\hbox{ r }}}_{k\cdot };\overline{\hbox{ r }}\mathrm{.})}{(G-1){\sum }_{jk}D({\hbox{ r }}_{jk};{\overline{\hbox{ r }}}_{k\cdot })}$

(7)

is asymptotically distributed like Fisher's F with (L − 1)(G − 1) numerator and (L − 1)(n·G + 1) denominator degrees of freedom under the null hypothesis of no difference among groups.

$\begin{array}{l}E\left[D\right({\hbox{ r }}_{jk};{\overline{\hbox{ r }}}_{k\cdot }\left)\right]=\frac{L-1}{N}\left(1-\frac{G-1}{n.}\right),\\ E\left[D\right({\overline{\hbox{ r }}}_{k\cdot };\overline{\hbox{ r }}\mathrm{.}\left)\right]=\frac{L-1}{{n}_{k}N}\left(1-\frac{1}{G}\right),\end{array}$

$E\left[D\right(\overline{\hbox{ r }}\mathrm{.};\hbox{ q }\left)\right]=\frac{L-1}{n.N}.$

To emphasize the parallels between this test and the usual test in analysis of variance, note that D(k·;··) is a between-group divergence (the divergence between the k-th group mean and the grand sample mean) and D(rjk;k·) is a within-group divergence (the divergence between the j-th sample in group k, and the mean forgroup k).

The null hypothesis is rejected for values of f larger than the appropriate critical value.

Example. Comparison of CD4+ and CD8+ T cell repertoires in healthy volunteers. We collected peripheral blood samples from healthy volunteers and fractionated them into CD4+ and CD8+ T cells as described in Section 2.1. We compared spectratypes between these two subsets in each of the TCRBV families for which complete spectratype could be obtained using the f statistic computed as given in Equation (7). The null hypothesis tested is that the spectratype population parameter vectors in the two subsets are identical, and thus the divergence between two sample mean parameter vectors is attributable entirely to the same sources as the within-subset spectratype variability. Figure 2 shows the histograms corresponding to two of the TCRBV families studied, TCRBV1 and TCRBV5. TCRBV1 is judged to exhibit differences between the subsets, and TCRBV5 is not (Table 1). Nevertheless, the observed difference between the CD4 and CD8 spectratypes is consistent between the two families, with CD8 cells favoring shorter CDR3 lengths in both.

### 4.2 Variable divergence from a given population parameter vector

A second relevant scenario is one in which we have two or more spectratypes, and our interest is in describing how much they each differ from a given population parameter vector q. This acquires additional salience when q is the parameter vector corresponding to the maximally diverse TCR repertoire. In that case, as we show in the Discussion below, the divergence from q is a direct measure of the corresponding repertoire diversity.

The measure of the departure from the population parameter vector is the parameter d defined in Equation (4). The analysis following that equation shows that log d is distributed approximately normally, with the log transformation regularizing the variance. The approach we take here, then, is to estimate d as appropriate to the model under investigation and use the logs of these estimates in parameter fitting and hypothesis testing. All the machinery of normal-model statistics including analysis of variance, regression and linear modeling, more generally, are then at our disposal.

Example. Development of a diverse TCR repertoire following thymus transplantation. The immune response to infectious agents involves the selective expansion of particular T cell specificities, a perturbation that transiently reduces the diversity of the TCR repertoire. Autoimmune disease is similarly accompanied by decreased TCR diversity, as is the uncontrolled expansion of T cells that defines leukemias. In each of these cases, the return to health and to a steady-state, implies a return to the prior state of diversity. This process can be described in terms of continuous-time statistical models.

In this example, we examine the establishment of diversity over time in a Di George syndrome patient (DIG 102; Markert et al., 2004) following thymus transplantation. We use a three-parameter piecewise linear model given byIn this model, the divergence is constant until time τ, when healthy T cells passing through the transplanted thymus begin to reach effective levels. τ is thus the regression breakpoint and is to be estimated using the data. H is the Heaviside function, having value one for non-negative argument and zero otherwise. The intercept α is treated as a random effect grouped by TCRBV family—different TCRBV families will have different starting d values. The slope β quantifies the rate at which the TCR repertoire diversifies. The errors, ε are independent and normally distributed.

$log{d}_{i}=\alpha +H({t}_{i}-\tau )\beta ({t}_{i}-\tau )+{\epsilon }_{i}.$

The parameter values ± their standard errors, estimated using data from DIG102 are,

$$\widehat{\alpha }=-0.159\pm 0.150$$

,

$$\widehat{\tau }=111\pm 5.4$$

, and

$$\widehat{\beta }=-4.21\times {10}^{-2}\pm 3.86\times {10}^{-3}$$

; the regression curve and data are shown in Figure 3.

## 5 DISCUSSION

We have provided a first-principles method for the comparison of antigen receptor spectratypes under two distinct sets of circumstances. Yet spectratype data are typically regarded as providing information about the diversity of the antigen receptor repertoire. Here, we comment on the relationship between our methods and diversity per se.

### 5.1 Kullback–Leibler divergence and totalTCR diversity

TCR diversity can be decomposed hierarchically, from distribution of TCRBV family usage, to distributions of CDR3 lengths within TCRBV family, to specific DNA sequence within CDR3 length. Alternative decompositions are possible, but these are the levels that correspond to conveniently available biological assays: family identification by flow cytometry, CDR3 length by spectratype, and specific DNA sequence by nucleotide sequencing. The data obtained in each of these assays represents the relative frequency of TCR counts in a subclass conditional on the parent class. We denote these relative frequencies ri, rj|i and rk|ij, respectively.

Each functional T cell clone, defined as the set of T cells that respond to the same antigenic peptide–MHC ligands, is restricted in size by intraclonal competition for TCR-specific growth signals. Assume, as a simplification, that each functional clone has the same maximum size. Let the number of functional clones in family i with CDR3-length class j be nij.

Then the expected proportion of TCR in the i-th family is

$\hbox{ E }\left[{r}_{i}\right]\equiv {q}_{i}=\frac{{n}_{i\cdot }}{n\mathrm{.}},$

(8)

the expected proportion of TCR in CDR3-length class j conditional on being in family i is

$E\left[{r}_{j|i}\right]\equiv {q}_{j|i}=\frac{{n}_{ij}}{{n}_{i\cdot }},$

(9)

and the expected proportion of TCR with sequence k conditional on CDR3 length j and family i, is

$E\left[{r}_{k|ij}\right]\equiv {q}_{k|ij}=\frac{1}{{n}_{ij}}.$

(10)

The total TCR repertoire diversity can be quantified by the entropy,which, using Equations (8)–(10), becomes

$\begin{array}{c}S=-{\displaystyle \sum _{ijk}}{r}_{i}{r}_{j|i}{r}_{k|ij}log\frac{{r}_{i}{r}_{j|i}{r}_{k|ij}}{n\mathrm{.}{q}_{i}{q}_{j|i}{q}_{k|ij}}\\ =logn\mathrm{.}-D-{\displaystyle \sum _{i}}{r}_{i}{D}_{i}-{\displaystyle \sum _{ij}}{r}_{i}{r}_{j|i}{D}_{j|i}\end{array}$

(11)

where the D's are the Kullback–Leibler divergences at different levels:and

$S=-{\displaystyle \sum _{ijk}}{r}_{ijk}log{r}_{ijk},$

$\begin{array}{l}D\equiv {\displaystyle \sum _{i}}{r}_{i}log\frac{{r}_{i}}{{q}_{i}},\\ {D}_{i}\equiv {\displaystyle \sum _{j}}{r}_{j|i}log\frac{{r}_{j|i}}{{q}_{j|i}},\end{array}$

${D}_{j|i}\equiv {\displaystyle \sum _{k}}{r}_{k|ij}log\frac{{r}_{k|ij}}{{q}_{k|ij}}.$

Spectratype divergences correspond to Di for each TCRBV family i. The first term in Equation (11), referring to the total number of (potential) T cell specificities is invariant. The variable part of the TCR repertoire entropy is given by the negative-weighted sum of the hierarchical Kullback–Leibler divergences.

Thus, when a good approximation of the ‘true’ population parameter vector is used, the analysis of deviation from this parameter vector, as described above, yields valid estimates of one component of the TCR repertoire diversity.

## 6 CONCLUSION

T cell receptor diversity is essential to the effective functioning of the immune system; careful measurement of the TCR diversity can provide valuable information on the state of the immune system. We have started from first principles and derived statistical methods that allow spectratype data to be quantified and used for hypothesis tests and parameter estimation. The primary statistic that arises, the Kullback–Leibler divergence, is generally a measure of difference between two probability functions, and in the appropriate context is a natural measure of deviation from maximum diversity.

The techniques we have described provide the objectivity and statistical power to further open new avenues for the application of spectratype analysis in both clinical and basic research areas.

We have made these methods publicly available via the World Wide Web at cbcb.duke.edu/SpA as described in greater detail by M.He, J.K.Tomfohr, B.H.Devlin, M.Sarzotti, M.L.Markert and T.B.Kepler (submitted for publication).

Fig. 1

Schematic diagram illustrating the spectratyping assay. T cells are fractionated from a sample of peripheral blood, the mRNA is isolated and TCRBV-family PCR primers are used to amplify part of the rearranged TCR containing the complete CDR3. This heterogeneous mixture of DNA molecules is then size-separated by electrophoresis. Finally, the quantity of material deposited in each band is quantified by densitometry, and the spectratype trace is produced.

Fig. 1

Schematic diagram illustrating the spectratyping assay. T cells are fractionated from a sample of peripheral blood, the mRNA is isolated and TCRBV-family PCR primers are used to amplify part of the rearranged TCR containing the complete CDR3. This heterogeneous mixture of DNA molecules is then size-separated by electrophoresis. Finally, the quantity of material deposited in each band is quantified by densitometry, and the spectratype trace is produced.

Fig. 2

## Abstract

T and B cell repertoires are collections of lymphocytes, each characterized by its antigen-specific receptor. We review here classical technologies and analysis strategies developed to assess immunoglobulin (IG) and T cell receptor (TR) repertoire diversity, and describe recent advances in the field. First, we describe the broad range of available methodological tools developed in the past decades, each of which answering different questions and showing complementarity for progressive identification of the level of repertoire alterations: global overview of the diversity by flow cytometry, IG repertoire descriptions at the protein level for the identification of IG reactivities, IG/TR CDR3 spectratyping strategies, and related molecular quantification or dynamics of T/B cell differentiation. Additionally, we introduce the recent technological advances in molecular biology tools allowing deeper analysis of IG/TR diversity by next-generation sequencing (NGS), offering systematic and comprehensive sequencing of IG/TR transcripts in a short amount of time. NGS provides several angles of analysis such as clonotype frequency, CDR3 diversity, CDR3 sequence analysis, V allele identification with a quantitative dimension, therefore requiring high-throughput analysis tools development. In this line, we discuss the recent efforts made for nomenclature standardization and ontology development. We then present the variety of available statistical analysis and modeling approaches developed with regards to the various levels of diversity analysis, and reveal the increasing sophistication of those modeling approaches. To conclude, we provide some examples of recent mathematical modeling strategies and perspectives that illustrate the active rise of a “next-generation” of repertoire analysis.

Keywords: diversity analysis, immune receptors, next-generation sequencing, modeling, statistics, gene nomenclature, B cell repertoire, T cell repertoire

## Introduction

T and B cell repertoires are collections of lymphocytes, each characterized by its antigen-specific receptor. The resources available to generate the potential repertoires are described by the genomic T cell receptor (TR) and immunoglobulin (IG) loci. TR and IG are produced by random somatic rearrangements of V, D, and J genes during lymphocyte differentiation. The product of the V-(D)-J joining, called the complementarity determining region 3 (CDR3) and corresponding to the signature of the rearrangement, binds the antigen and is responsible for the specificity of the recognition. During their differentiation, lymphocytes are subjected to selective processes, which lead to deletion of most auto-reactive cells, selection, export, and expansion, of mature T and B cells to the periphery. Primary IG and TR repertoires are therefore shaped to generate the available peripheral or mucosal repertoires. In addition, several different functional T and B cells subsets have been identified, with differential dynamics and antigen-specific patterns. These available repertoires are dramatically modified during antigen-driven responses especially in the inflammatory context of pathogen infections, autoimmune syndromes, and cancer to shape actual repertoires. When considering the importance of efficient adaptive immune responses to get rid of infections naturally or to avoid auto-reactive damages, but also for therapeutic purposes such as vaccination or cell therapy, one realizes the relevance of understanding how lymphocyte repertoires are selected during differentiation, from ontogeny to aging, and upon antigenic challenge. However, immune repertoires of expressed antigen receptors are built by an integrated system of genomic recombination and controlled expression, and follow complex time-space developmental patterns. Thus, an efficient repertoire analysis requires both (1) methods that sample and describe the diversity of receptors at different levels for an acceptable cost and from a little amount of material and (2) analysis strategies that reconstitute the best multidimensional picture of the immune diversity from the partial information provided by the repertoire description as reviewed in Ref. (1). In the following sections, we summarize technologies developed over the past decades to describe lymphocyte repertoires and we present the growing number of analysis tools, evolving from basic to sophisticated statistics and modeling strategies with regards to the level of complexity of the data produced.

## Methods Developed to Describe the IG and TR Repertoires

B and T lymphocyte repertoires can be studied from different lymphoid tissues and at various biological levels, such as cell membrane or secreted proteins, transcripts or genes, according to the techniques used. Fluorescence microscopy or flow cytometry techniques allow to track and sort particular cell phenotypes and to quantify the expressed repertoire at the single-cell level with V subgroup-specific monoclonal antibodies. Alternatively, the IG or TR diversity may be also analyzed using proteomics methods from either the serum (for IG) or dedicated cell extracts. Finally, molecular biology techniques assess the repertoire at the genomic DNA or transcriptional levels, qualitatively and/or quantitatively.

### Analysis of IG and TR repertoires at the protein level

#### Flow cytometry single-cell repertoire analysis

The frequency of lymphocytes expressing a given IG or TR can be determined using flow cytometry when specific monoclonal antibodies are available. This technique allows for the combined analysis of the antigen receptor and of other cell surface markers. Currently, using flow cytometry, up to 13 parameters can be routinely studied at once, reaching 20 parameters with the last generation flow cytometers and 70–100 parameters with mass cytometry (2). Seminal studies in mice using specific anti-TRBV antibodies have led to the characterization of the central tolerance selection processes that occur in thymus (3–5). Later on, a comprehensive description of the human TRBV repertoire was setup (6), when monoclonal antibodies became available for most of the TRBV subgroups. Repertoire analysis with flow cytometry provides a qualitative and quantitative analyses of the variable region, often done on heterogeneous cell populations, in order to decipher, for example, selection events related to aging, perturbations, and treatments (7). However, this technology is naturally limited by the availability of specific monoclonal antibodies, and does not address more detailed issues such as junction diversity. Furthermore, polymorphism of the IG or TR genes (8, 9) may constitute a serious limitation for a systematic survey using these approaches.

#### Proteomic repertoire analysis for serum immunoglobulins

Recent developments of proteomics tools now offer sensitivity levels applicable to IG repertoire analysis. Such a description at the protein level takes into account all post transcriptional and translational modifications.

##### PANAMA-blot technology

A semi-quantitative immunoblot, called the PANAMA-blot technique (10), allows for the identification of the antibody reactivities present in collection of sera (or cell culture supernatant) against a given source of antigens (10–12). Briefly, a selected source of antigens is subjected to preparative SDS-PAGE, transferred onto nitrocellulose membranes, then incubated with the serum to be tested allowing for the revelation of the bound antibodies using an appropriate secondary antibody coupled to alkaline phosphatase. Computer-assisted analysis of the densitometric profiles allows for the rescaling and the quantitative comparison of patterns of antibody reactivity from individuals in different groups. A large amount of data is generated when testing a range of sera against various sources of antigens. Statistical analyses are included in the PANAMA-Blot approach (as described further). This global analysis helped to reveal that the IgM repertoire in mice is selected by internal ligands and independent of external antigens (13).

This method can also lead to identify IG reactivity patterns specific for a type of pathology or clinical status and has been applied to both fundamental and clinical analysis. In particular, it was used to analyze human self-reactive antibody repertoires and their potential role for down-modulating autoimmune processes (14–16).

##### Antigen micro-array chips

More recently, antigen micro-array-based technology coupled to a complex two-way clustering bioinformatics analysis was developed to evaluate the serum repertoire antibodies from diabetes-prone individuals and revealed their predictive or diagnostic value. In brief, a range of antigens (proteins, peptides, nucleotides, phospholipids…) were plated onto glass plates and incubated with sera from individuals (human diabetes patients or mice in an experimental model of diabetes). The intensity of reactivity of the serum IG for each peptide was determined and scored against the control reactivity. Clustering analysis was then implemented to determine a potential antigen signature that significantly sorts out diabetes from non-diabetes individuals. In this way, it was found that the patterns of IgG antibodies expressed early in male NOD mice can mark susceptibility or resistance to diabetes induced later and that it is different than the pattern characteristic of healthy or diabetic mice after disease induction (17). Similarly, this clustering approach was applied in humans to successfully separate human subjects that are already diabetic from healthy people (18).

### Repertoire analysis at the genomic DNA level

Other strategies that cover IG or TR repertoire analyses have been developed at the genomic DNA level. Firstly, CDR3 spectratyping studies (detailed in the following section) have been carried out at the DNA level mostly to address issues related to B or T cell development (19, 20). More recently, an original multiplex genomic PCR assay coupled to real-time PCR analysis was developed to provide a comprehensive description of the mouse T cell receptor alpha (TRA) repertoire during development (21). Although these approaches can be applied to all IG isotypes and TR, they have not been used as much as transcript CDR3 spectratyping due to sensitivity and heterozygosity issues.

Immunoglobulin or T cell receptor repertoires can also be assessed by following the diversity of rearrangement deletion circles. Since they are produced by the V-(D)-J recombination machinery when the joint signal is formed and diluted in daughter cells, they give a good representation of recently generated T or B cells. This technique has been particularly useful for describing the restoration of T cell diversity following highly active antiretroviral therapy in HIV-infected patients (22) and has been used to model thymic export (23, 24) as well as to demonstrate continued contribution of the thymus to repertoire diversity, even in older individuals (25). It also reveals that thymic output is genetically determined, and related to the extent of proliferation of T cells at DN4 stage in mice (26). However, their analysis does not provide much insight into the level of diversity since the signal joint does not vary for a given combination of genes. Therefore, the interest of such analyses is reached when combined with CDR3 spectratyping analyses to know whether a repertoire perturbation is rather attributable to newly produced T cells or peripheral T cell proliferation.

### V-(D)-J junction analysis of IG and TR transcript repertoires

Original molecular-based strategies for analyzing repertoire diversity relied on cloning and hybridization of molecular probes specific for IGHV gene subgroups first by RNA colony blot assay (27). This led to the observation that IGHV gene usage is characteristic of mouse strain and is a process of random genetic combination by equiprobable expression of IGHV genes (28). The study of selection processes revealed that the IGHV region-dependent selection determines clonal persistence of B cells (29) and that selection with age leads to biased IGHV gene expression (30).

In situ hybridization on single-cells revealed that during mouse ontogeny and early development of B cells in bone marrow, there is a non-random position-dependent IGHV gene expression, favoring D-proximal IGHV gene subgroup usage (31). Thereafter, sequencing of PCR-amplified cDNA collections were obtained from samples of interest. Although fastidious, these early studies have been useful in defining the basis of human IG and TR repertoires in terms of overall distribution, CDR3-length distribution, and V-(D)-J use (32–35), sometimes leading to the identification of new IG or TR genes. Later, more practical techniques have been developed for large-scale analysis of lymphocyte repertoires, such as quantitative PCR, micro-array, and junction length spectratyping, as described below.

#### Quantitative RT-PCR for repertoire analysis

In parallel to qualitative CDR3 spectratyping techniques (see section below), quantitative PCR strategies were developed (36). Coupling the two techniques for all V domain-C region combinations provides a complete qualitative and quantitative picture of the repertoire (37–39) described by up to 2,000 measurements per IG isotype or TR for one sample. With the development of real-time quantitative PCR, this approach opened the possibility for a more precise evaluation of repertoire diversity (39–41). Complementary tools have been also developed in order to allow normalization of spectratype analysis such as studies by Liu et al. (42) and Mugnaini et al. (43).

Matsutani et al. (44) developed another method to quantify the expression of the human TRAV and TRBV repertoires based on hybridization with gene specific primers coated plates. The cDNA from PBMC extracted RNA are ligated to a universal adaptor which allows for a global amplification of all TRAV or TRBV cDNAs. The PCR products are then transferred onto microplates coated with oligonucleotides specific for each TRAV or TRBV regions, and the amount of hybridized material is quantified. This technique was used to analyze the TR repertoire diversity of transplanted patients (45) and adapted to the study of mouse TRAV and TRBV repertoires (46). VanderBorght et al. also developed a semi-quantitative PCR-ELISA-based method for the human TRAV and TRBV repertoire analysis (38). The combined usage of digoxigenin (DIG)-coupled nucleotides and DIG-coupled reverse TRAC or TRBC primers allowed for a quantitative measurement of the amount of amplified DNA by a sandwich ELISA.

Du et al. (47) later setup a megaplex PCR strategy to characterize the antigen-specific TRBV repertoire from sorted IFNγ-producing cells after Mycobacterium infection. The clonotypic TRBV PCR products were used for Taqman probes design to quantify the expression of the corresponding clonotypes from ATLAS-amplified SMART cDNAs.

#### Direct measurement of lymphocyte diversity using micro-arrays

Another technology, similar to the one just discussed, has been developed by the group of Cascalho et al. which allows for a direct measurement of the entire population of lymphocyte-receptors. This is accomplished by hybridization of lymphocyte-receptor specific cRNA of a lymphocyte population of interest to random oligonucleotides on a gene chip; the number of sites undergoing hybridization corresponds to the level of diversity. This method was validated and calibrated using control samples of random oligonucleotides of known diversity (1, 103, 106, 109) (48, 49) and successfully demonstrated that central and peripheral diversification of T lymphocytes is dependent on the diversity of the circulating IG repertoire (49, 50). Similarly, a highly sensitive micro-array-based method has been proposed to monitor TR repertoire at the single-cell level (51).

#### CDR3 spectratyping techniques

##### Immunoscope technology

Among various techniques used to analyze the T or B cell repertoires, Immunoscope, also known as CDR3 spectratyping (52, 53) consists in the analysis of the CDR3-length usage so that antigen-specific receptor repertoires can be described by thousands of measurements. In the case of naive murine repertoires, T cell populations are polyclonal and analysis typically yields eight-peak regular bell-shaped CDR3 displays (wrongly assumed to be Gaussian), each peak corresponding to a given CDR3-length. When an immune response occurs, this regular polyclonal display can be perturbed: one can see one or several prominent peaks that correspond to the oligoclonal or clonal expansion of lymphocytes. A complete description of this technique and its applications to clinical studies has been published elsewhere (54).

In the original Immunoscope publication, Cochet et al. (55) analyzed the T cell repertoire after the immunization of mice with the pigeon cytochrome c. They provided the first description of an ex vivo follow-up of a primary T cell specific response in a mouse model. Their second paper analyzed the average CDR3-lengths as a function of TRBV-TRBJ combinations. In particular, the authors found a correlation between TRBV CDR1 and major histocompatibility (MH) haplotype (52). This group later published a large amount of original studies in various models such as lymphocyte development (40, 56–63), kinetics of antigen-specific responses (64–67), viral infection (68, 69), autoimmunity (70, 71), tumor-associated disease (72), and analysis of allogeneic T cell response and tolerance after transplantation (73). Notably, the combination of CDR3 spectratyping with flow cytometry-based IG or TR V frequency analysis provides a more comprehensive assessment, such as in Pilch et al. (74). For example, such an approach revealed the constriction of repertoire diversity through age-related clonal CD8 expansion (75). Similarly, a combination of CDR3 spectratyping, flow cytometry, and TR deletion circle analysis has allowed to define age-dependent incidence on thymic renewal in patients (76) or to evaluate the effects of caloric restriction in monkeys to preserve repertoire diversity (77). CDR3-length spectratyping was also used in other models, such as rainbow trout, to analyze TRB repertoire and its modifications induced by viral infection (78–80). While no tool such as monoclonal antibodies to T cell marker(s) was available in this model, this approach demonstrated that fish could mount specific T cell responses against virus, which could be found in all individuals (public clonotypes) or not (private clonotypes). Similar strategies, developed by other groups (81) and following the same approach in parallel, analyzed the IG repertoire in Xenopus at different stages of development, describing a more restricted IG junction diversity in the tadpole compared to the adult.

Gorski et al. (82) developed their own CDR3 spectratyping technique to analyze the complexity and stability of circulating αβT cell repertoires in patients following bone marrow transplantation as compared to normal adults. They showed that repertoire complexity of bone marrow recipients correlates with their state of immune function; in particular, individuals suffering from recurrent infections associated with T cell impairment exhibited contractions and gaps in repertoire diversity. The detailed procedure for this technique has been published in Maslanka et al. (83). A variation of this technique has been reported later by Lue et al. (84), relying on a compact glass cassette, a simpler device than the usual automated plate DNA sequencers.

##### Alternative technologies

Alternative CDR3 spectratyping techniques have been described such as single-strand conformation polymorphism (85–87) and heteroduplex analysis (88–91). These methods differ from the CDR3 spectratyping/Immunoscope technique mostly in the way PCR products are analyzed by performing non-denaturing polyacrylamide electrophoresis. The main advantage of these techniques is a more direct assessment of clonal expansion since PCR products migrate according to their conformation properties; therefore, presence of a predominant peak is strongly indicative of clonality when a smear migration pattern indicates polyclonality. However, these techniques have been less widely used probably because of the difficulty to make clear correlations between the expanded peaks across samples.

Another original alternative technique has been described by Bouffard et al. (92), analyzing products obtained after in vitro translation of PCR-amplified TR-specific products by isoelectric focusing. With this technique, clonality can also directly be assessed by looking at the obtained migration profile.

### IG/TR rearrangement sequencing: From cloning-based- to next-generation-sequencing

In order to get a better description of IG/TR diversity at the nucleotide sequence level, thus providing fine-tuned description of the actual diversity, Sanger sequencing approaches relying on bacterial cloning of rearrangements were performed in physiological conditions globally (60, 93–99) or partially to characterize particular expansions identified by other technologies such as CDR3 spectratyping (40, 59, 100–102), flow cytometry (103). They were also used in pathological/infectious conditions (104–107) sometimes leading to antigen-specific T cell TR identification and quantification through the combination of antigen-specific T cell stimulation and cytometry-based cell sorting, anchor-PCR, and bacterial cloning-based sequencing (108).

These studies pioneered the description of the repertoire and provided fruitful information regarding the extent and modification of the diversity. However, besides being time and cost-extensive, such approaches have allowed for the analysis of 102–103 sequences, far under the estimated diversity reaching 106–107 unique clonotypes in mice and humans (40, 59, 109).

In the last decade, DNA sequencing technologies have made tremendous progresses (110) with the development of so called next-generation sequencers, already reaching four generations (111). Those instruments are designed to sequence mixtures of up to millions of DNA molecules simultaneously, instead of individual clones separately. Second generation sequencers became affordable in the last 5 years and have been used for immune repertoire analysis, starting with the seminal work of Weinstein et al. (112) where the IG repertoire of Zebrafish has been described by large-scale sequencing. Consequently, exploratory works by other groups provided an overview of the complex sequence landscape of immune repertoires in humans (113–118). More recent work aimed at addressing fundamental questions such as lineage cells commitment (119–122), generation of the diversity processes (123–125), and diversity sharing between individuals (126, 127). Finally, the power of this technology has been validated in the clinic as well (128, 129).

As seen above for other technologies, combinations of approaches have been applied to NGS. Notably, deep sequencing has been used in combination with CDR3-length spectratyping by some groups to study human (130) or rainbow trout IG (131) repertoire modifications after vaccination against bacteria or viruses. In the latter, pyrosequencing performed for relevant VH/Cμ or VH/Cτ junctions identified the clonal structure of responses, and showed, for example, that public responses are made of different clones identified by (1) distinct V-(D)-J junctions encoding the same protein sequence or (2) distinct V-(D)-J sequences differing by one or two conservative amino acid changes (131) as described for public response in mammals (132, 133). These studies showed that NGS and traditional spectratyping techniques lead to remarkably similar CDR3 distributions.

With the power of such approach for genomics and transcriptomics studies in general, constant improvements are achieved to increase the sequencing deepness and read length as well as to reduce the cost, therefore offering multitude of biological explorations (139). NGS now permits a comprehensive and quantitative view of IG and TR diversity by combining and improving the sensitivity of classical approaches with accurate and large-scale sequencing. NGS has the power to identify IG or TR specific for given antigens (in combination with antigen-specific assays) and to define more complex signatures (i.e., TR sets) related to disease and/or treatment from heterogeneous T and B cell populations. Still, most of the deep sequencing efforts have been limited to only one chain of the receptor at the repertoire level (usually the β chain for TR and the heavy chain for IG). Indeed, current high-throughput approaches do not allow one to assign which combination of chains (TRA and TRB, or IGH and IGK or IGL) belong to which cell (140). A recent development by DeKosky et al. proposed a reasonably high-throughput technology to assess massively paired IG VH and VL from bulk population (141). In parallel, Turchaninova et al. (142) have proposed a similar approach for the paired analysis of the TRA and TRB chains. The parallel development of high-throughput microfluidic-based single-cell sorting will certainly push forward new developments in the field (143).

However, despite the technological advance, studies so far have mainly reported CDR3 counting and identification of major expansions. The complexity of immune repertoires is still a matter that such approach cannot completely overcome, due to the paucity of powerful analytical methods. Besides data management tools, studies are now starting to extract most of the benefit from such approach to model the immune repertoire diversity and dynamics (144), an approach that may help in understanding the interplay between cells and repertoire shaping. Accurate and powerful statistical analyses are required to manage such amount of information. Current state will be reviewed in the following sections.

## Potential and Genomic Repertoires: A Question of Ontology and Orthology

Immune repertoires sensu stricto are expressed by lymphocyte clones, each carrying a single receptor for the antigen. Such receptors comprise IG and TR in jawed vertebrates (8, 9) and VLR in Agnathans (145). The sequences of these receptors are available in databases such as GenBank or EMBL, which are difficult to use for transversal studies due to inconsistent annotation. The IMGT® information system (see below) has largely solved this problem setting standardized gene nomenclatures, ontologies and a universal numbering of the IG/TR V and C domains, thus giving a common access to standardized data from genome, proteome, genetics, two-dimensional, and three-dimensional structures (146). The accuracy and the consistency of the IMGT® data are based on IMGT-ONTOLOGY, the first, and so far, unique ontology for immunogenetics and immunoinformatics (147).

With the development of high-throughput sequencing, large numbers of new sequences of antigen receptor genes have become available, which can be classified into different categories: genomic sequences of IG or TR (in germline configuration in genome assemblies) or fragments of IG/TR transcripts, containing the CDR3 or not. Also, these datasets can be produced from species newly sequenced, as well as from new haplotypes of well-described species.

The annotation of such sequences remains an open question. Manual annotation is not applicable, and no good automated approach has been validated yet. A relevant annotation of these massive datasets will require the integration of genomic and expression data with existing standardized description charts, as offered by IMGT®. A standardized annotation is an important issue since it facilitates the re-utilization of datasets and comparison of analyses. Thus, the description of IG and TR polymorphisms, the integration of repertoire studies with structural features of antigen-specific domains, and even the usage of new genes in genetic engineering rely on a common standard for nomenclature, numbering, and annotation (147).

To take advantage of the current standards that have been established from classical sequencing data during the last 25 years, new, fast, reliable, and human-supervised annotation methods will have to be developed, integrating directly high-throughput sequence information from the increasing number of deep sequencing platforms and technologies, at different genetic levels (genome, transcriptome, clonotype repertoires). Along this line, IMGT/HighV-QUEST offers online tools to the scientific community for the analysis of long IG and TR sequences from NGS (148).

Special attention can be paid to the orthology/paralogy relationships between similar antigen receptor genes from different species. These characteristics are essential to understand the dynamics of IG and TR loci. In fact, with many important lymphocyte subsets characterized by canonical/invariant antigen receptors, such relationships are critical to transfer functional knowledge between models. Importantly, the phylogenetic analyses required to reconstitute the evolution of antigen receptor genes are based on multiple alignments, the quality of which is highly dependent on common numbering and precise annotation of sequences.

As far as immune repertoires are characterized by the diversity of receptors specifically binding antigen/pathogen motifs to initiate a defense response, they might not be limited to lymphocyte diversifying receptors, e.g., IG, TR, and VLR. The particularity of these systems is a somatic diversification combined to a clonal structure of the repertoire, each lymphocyte clone expressing the product of a recombination/hypermutation and/or conversion process. However, many other arrays of diverse receptors binding or sensing pathogens have been discovered in metazoans, in invertebrates as well as in vertebrates.

In some cases, their diversity is really “innate,” i.e., encoded in the genome as multiple genes produced by duplications. Fish NLR, finTRIMs, and NITR, primate KIR, chicken CHIR, or TLR in sea urchin, constitute good examples of such situations. While these repertoires may appear as relatively limited, polymorphism within populations, and differential expression of receptors per cell upon stimulation represent complex issues, which fall well into “traditional” repertoire approaches.

In other cases, receptors are subject to diversification processes much faster than gene duplication, which does not comply with a clonal selection pattern. The best examples are probably the DSCAM in arthropods, which hugely diversify by alternative splicing of exons encoding half-IgSF domains (149, 150), and the FREP lectins in mollusks, of which sequences are highly variable at the population level, and even between parents and offspring produced by auto-fecundation (151).

The number of such “innate” repertoires which are not expressed by clonally selected lymphocytes will likely increase with deep sequencing of new genomes/transcriptomes, as illustrated by a recent report from mussel (152). A good example of the importance of a proper structural description of key domains of receptors is provided by the extensive analysis of LRR motifs in studies on TLR evolution (153, 154). Further insights into the functions of such diverse proteins will be provided by the characterization of their expressed (available) repertoire, at different levels such as single-cells, cell populations, and animal populations.

Such analyses will require precise identification of genes and sequences as well as mutations, and a standardized approach of nomenclature and structural description will be as useful as it is for the vertebrate IG and TR sequences. Importantly, these receptors are made of a small number of structural units, such as IgSF domain or LRR domains, which suggests that standardized system(s) for sequence annotation could be developed following IMGT standards (155).

## Statistical Analysis and Modeling of Immune Repertoire Data

### Statistical repertoire analysis

The description of the repertoire modifications using flow cytometry or Immunoscope provided clear-cut and detailed insight into the clonal expansion processes during the responses against a defined antigen (64, 66). However, it is difficult to identify the relevant alterations of the repertoires in more complex situations such as pathogen infections or variable genetic backgrounds. For example, it appeared impossible to identify all significant modifications of TRB Immunoscope profiles during cerebral malaria by direct ocular comparison (107). Different methods were therefore developed to extract from IG and TR repertoire descriptions the relevant information, to encode it as numerical tables and to analyze them with statistical models.

#### CDR3 spectratype perturbation indices

Since the initial description of the CDR3 spectratyping technique, different scoring indices were developed or derived from the literature: “relative index of stimulation” (RIS) (55), “overall complexity score” (156), Reperturb (157), “complexity scoring system” (158), COPOM (159), Oligoscore (160), TcLandscape (161), “spectratype diversity scoring system” (162), Morisita-Horn index and Jaccard index (95–97), “absolute perturbance value” (163). A comparative review of such scoring strategies was published by Miqueu et al. (164).

In particular, the perturbation index Reperturb was developed by Gorochov et al. to perform TR repertoire analysis in HIV patient during progression to AIDS and under antiretroviral therapy. They could show drastic restrictions in the CD8+ T cell repertoire at all stages of natural progression that persisted during the first 6 months of treatment. In contrast, CD4+ T cell repertoire perturbations correlated with progression to AIDS with a return to a diversified repertoire in good responders to treatment (157).

Soulillou et al. refined this approach by combining the qualitative information obtained with usual CDR3 spectratyping with quantitative information of TRBV usage obtained by real-time quantitative PCR. They devised a four-dimension representation that represents TRBV subgroups, CDR3-length and percentage of TRBV use on three axis chart in addition to a color-coded representation of the CDR3 profile perturbation. Using this original approach, they were able to show that graft rejection is associated with a vigorous polyclonal accumulation of TRBV mRNA among graft-infiltrating T lymphocytes, whereas in tolerated grafts T cell repertoire is strongly altered (161, 165). Their study puts the emphasis on the importance of not only qualitative but also quantitative analysis of lymphocyte repertoires.

#### Platforms for repertoire data management and statistical analysis

Several platforms have been developed and rely mostly on CDR3 spectratyping and sequencing data, with recent developments to manage and analyze NGS data.

The ISEApeaks strategy and software were developed in order to satisfy the needs for efficient automated electrophoresis data retrieval and management (160, 166). ISEApeaks extracts peak area and length data generated by software used to determine fragment intensity and size. CDR3 spectratype raw data, consisting of peak areas and nucleotide lengths for each V-(D)-J-C combination, is extracted, smoothed, managed, and analyzed. The repertoires of different samples are gathered in a peak database and CDR3 spectratypes can be analyzed by different perturbation indices and multivariate statistical methods implemented in ISEApeaks. We have applied our ISEApeaks strategy in several studies. In an experimental model of cerebral malaria, we established a correlation between the quality of TR repertoire alterations and the clinical status of infected mice, whether they developed cerebral malaria or not (107). We contributed to the characterization of the membrane-associated Leishmania antigens (MLA) that stimulates a large fraction of naive CD4 lymphocytes. Repertoire analyses showed that MLA-induced T cell expansions used TR with various TRBV rearrangements and CDR3 lengths, a feature closer to that of polyclonal activators than of a classic antigen (167). We also revealed repertoire age-related perturbations in mice (7). ISEApeaks functions for statistical analysis was successfully applied to analyze the TR repertoire in fish as shown by our detailed analysis of the TRB repertoire of rainbow trout IELs, performed in both naive and virus-infected animals. Rainbow trout IEL TRBV transcripts were highly diverse and polyclonal in adult naive individuals, in sharp contrast with the restricted diversity of IEL oligoclonal repertoires described in birds and mammals (102). More recently, our study of the CD8+ and CD8 αβ T cell repertoire suggests different regulatory patterns of those T cell patterns in fish and in mammals (168). ISEApeaks was also used to implement a new statistically based strategy for quantification of repertoire diversity (159).

Kepler et al. described another original statistical approach for CDR3 spectratype analysis, using complex procedures for testing hypotheses regarding differences in antigen receptor distribution and variable repertoire diversity in different treatment groups. This approach is based on the derivation of probability distributions directly from spectratype data instead of using ad hoc measures of spectratype differences (169). A software (called SpA) implementing this method has been developed and made available online (170). This approach has been used in a longitudinal analysis of TRBV repertoire during acute GvHD after stem cell transplantation (171).

Another group (163) reported the development of a new software platform, REPERTOIRE, which allows handling of CDR3 spectratyping data. This software implements a perturbation index based upon an expected normal Gaussian distribution of CDR3 length profiles.

Owing to the complexity and diversity of the immune system, immunogenetics represents one of the greatest challenges for data interpretation: a large biological expertise, a considerable effort of standardization, and the elaboration of an efficient system for the management of the related knowledge were required. To answer that challenge, IMGT®, the international ImMunoGeneTics information system®(http://www.imgt.org), was created in 1989 by one of the authors (146). Overtime, it developed standards that, since 1995, have been endorsed by the World Health Organization-International Union of Immunological Societies (WHO-IUIS) Nomenclature Committee and by the WHO-International Nonproprietary Names (INN) (172–175). IMGT® comprises seven databases (sequence, gene, and structure databases), 17 online tools and more than 15,000 pages of web resources. Among the databases, IMGT/LIGM-DB, the database for nucleotide sequences (170,685 sequences from 335 species as of July 2013) and IMGT/GENE-DB, the gene database (3,081 genes and 4,687 alleles) are of great interest for repertoire analysis. Freely available since 1997, IMGT/V-QUEST is an integrated system for the standardized analysis of collections of IG and TR rearranged nucleotide sequences (176, 177). A high-throughput version, IMGT/HighV-QUEST (148), has been released in 2010 for the analysis of long IG and TR sequences from NGS using the 454 Life Sciences technology. In the same line, other analysis tools are becoming available showing the renewed interest for repertoire analyses and modeling consecutive to NGS technology developments (178–181).

Altogether, these efforts highlight the relevance of developing more efficient and powerful technologies for the evaluation of repertoire diversity. Notably, two successful French biotech companies (TcLand, Nantes; ImmunID, Grenoble) were created in the field of repertoire analysis, using different technologies. In collaboration with ImmunID, we have proposed a novel strategy for statistical modeling of T lymphocyte repertoire data obtained in humans and humanized mice. With this model, we revealed that half of the human TRB repertoire, in terms of proportion of TRBV-TRBJ combinations, is genetically determined, the other half occurring stochastically (182). In addition, the biotechnology company “Adaptive” and the “Repertoire 10K (R10K) Project” have been recently founded by researchers respectively from the Fred Hutchinson Cancer Research Center (Seattle and Washington) and the HudsonAlpha Institute (Huntsville). Both have developed platforms (immunoSEQ®, iRepertoire®) providing researchers with a global analysis of the T or B cell receptor sequence repertoires (183). However, despite the power of this technology, studies are still limited by the ability to process the complexity of the information provided. Specific software developments for the automatic treatment and annotation of IG and TR sequences and the statistical modeling of repertoire diversity can still be improved.

#### Multivariate analysis

As mentioned above, the PANAMA-Blot technique also includes statistical analysis of the data. Multi-parametric analysis was introduced to compare the global reactivity of antibodies of different individuals in different groups with a given antigenic extract. This analysis has been successfully implemented to identify reactivity patterns specific for a given pathology or clinical status (10–12, 14, 15, 184). Similarly, multi-parametric analysis was also applied to TRBV spectratype analysis in an experimental cerebral malaria model (107).

Hierarchical clustering or classification algorithms have become very popular with the growing of micro-array-based transcriptome analysis. Although still uncommon for immune repertoire analysis, such approaches have been employed to categorize large sets of repertoire data without a priori (17, 102, 107).

#### Diversity indices

The concept of immune repertoire has been devised to describe the diversity of cells involved in the immune system of an individual (1). As described above, different scoring systems were developed to assess this diversity, some are heuristics but others have been borrowed from theoretical ecology and evolution. As reviewed by Magurran (185), the Shannon entropy, introduced by Claude Shannon in 1948 for the information theory, is the most used because it not only integrates the number of different species but also the relative proportion of each of these species. In 1961, Alfred Rényi generalized this entropy to a family of functions, like Species Richness, Simpson, Quadratic, and Berger–Parker indices, for quantifying the diversity, the uncertainty or randomness of a system. Most of these indices are implemented in the free software application Estimates (http://purl.oclc.org/estimates) (186). Altogether, these diversity indices constitute a collection of tools with their own sensitivity to the variety and the relative abundances of the species that are perfectly suitable for assessing immune repertoire diversity. Indeed, the very famous index of variability proposed by Kabat and Wu (187) corresponds to the ratio of Species Richness and Berger–Parker indices. In 1990, Jores et al. showed that the resolving power of this Wu-Kabat variability coefficient can be enhanced by increasing the weight on the frequency distribution of the amino acids in the formula (188). This approach inspired Stewart et al. (189) to use the Shannon entropy to demonstrate that TR amino acid composition is significantly more diverse than that of IG. In the same way, CDR3 spectratyping data can be analyzed using the relative abundance of each peak within CDR3 length global distribution. By doing so, we adjusted the original Shannon entropy, making it reaching its maximum for a Gaussian distribution, to compare the CDR3 length diversity of splenic IgM, IgD, and IgT in infected Teleost Fish (131). Recently, the Gini index, used in ecology or economics to measure the equality of distributions, was applied to individual TR clones and compared naive and memory repertoires (190). The development of deep sequencing techniques ignited a renewed interest in IG/TR repertoire. Indeed, several studies used high-throughput analysis to describe TR repertoire of key T cell subsets in human peripheral blood (115, 126, 191). This approach assessing the repertoire diversity from the relative abundance of each species in the global distribution can be decomposed hierarchically into components attributable, respectively, to variations in TRBV-TRBJ combinations and in CDR3-length (113, 117). However, most of these studies have been limited to the counting of the observed unique clonotypes. Beside the species richness, ecology-derived indices have also been applied to assess and compare immune repertoire diversity. Föhse et al. (119) used the Morisita-Horn similarity index to compare regulatory T cell repertoires between several lymphoid organs. In addition, Simpson diversity index, associated with Shannon entropy, was used to monitor TR repertoire diversity of HIV-specific CD8 T cells during antiretroviral therapy (192) but also to quantify TR repertoire recovery in the blood after allogeneic hematopoietic stem cell transplantation (128). In the same manner, Koning et al. (193) used Shannon’s and Simpson’s indices to show the role for the peptide component of the peptide-MH1 complex on the molecular frontline of CD8+ T cell–mediated immune surveillance, by comparing the repertoire diversity of CD8+ T cell populations directed against a variety of epitopes. In parallel, using Simpson’s index as a metric allowed Johnson et al. (194) to model mathematically the naïve CD4 T cell repertoire contraction with age leading them to conclude that diversity plummet observed around the age of 70 could be correlated to cell-intrinsic mutations affecting cell division rate or death.

### Modeling strategies

Modeling approaches have a strong tradition in immunology, usually at the boundary with other disciplines such as physics (195). Before deep sequencing data was available, general design principles were proposed as desirable features of immune repertoires, with implications for the observed repertoire diversity and dynamics (196–198). Many efforts have involved the modeling of immune cell dynamics and the effects of antigens on repertoire diversity, using differential equations descriptions of the population dynamics (199–201). Recognition in the immune system is often studied both theoretically and experimentally by probing the dynamics of cells with a specific type of receptor with respect to infections (202). Alternatively one can look at the response of a small set of chosen receptors to a specific pathogenic challenge, or careful biochemical investigation of particular receptor/antigen pairs (203, 204). Much work has been devoted to systems-biology approaches to signal processing in immune cells, as reviewed in Germain et al. (205) and Emonet and Altan-Bonnet (206). Here we focus on approaches inspired by recent advances in sequencing technologies (112, 113, 115, 116, 125, 191, 207, 208) that have opened the way for data-driven modeling of the immune repertoires and interactions between receptors and antigen.

A common modeling approach for describing receptors at the amino acid level is to choose a relevant interaction parameter (e.g., chemical affinity or hydrophobicity) and assign it a simplified digit-string representation (209). These methods are extensions of the string model, which describes both receptor and epitopes as strings of length L, with values chosen from natural numbers, and quantify their interaction by the match between the two strings (197, 210, 211). Such quantitative, physically inspired descriptions of immune receptors, despite the arbitrary choice of interaction coordinates, have proven a valuable first step in statistically describing recognition in T cells (195, 212–215). Recently, lower hydrophilicity of regulatory vs. conventional T cells was suggested from CDR3 sequencing (216).

High-throughput sequencing of immune receptors raises specific challenges compared to traditional genomic sequencing. It is harder to distinguish sequencing errors from new polymorphisms, since no corresponding pre-existing sequence exists. One of the most interesting regions when studying diversity is the CDR3 with its many insertions and deletions added to the germline sequence. These regions are often hard to align to the genomic templates, or with each other (217). Therefore, extra care is needed when generating and analyzing sequence data. Not all sequencing technologies are equally good for all purposes (218): while 454 sequencing gives longer reads than Illumina it is known to have a greater probability of frameshift errors. In addition, primer-dependent PCR amplification biases require that raw sequence counts be normalized using control experiments (112) in order to accurately report clone sizes, as demonstrated by spike-in experiments (219). In TR repertoire studies, this is circumvented by using 5′RACE which provides an unbiased amplification of fully rearranged sequences, as recently demonstrated for TRB V-(D)-J transcripts (191).

Despite sequencing issues, statistical algorithms are often able to extract information from the data. Many studies of diversity focus on the V, D, and J gene usage of each rearranged sequence. Algorithms and tools have been developed to rapidly identify the V, D, and J genes for massive numbers of sequences (148, 178, 181). In many cases however, the assignment of a D gene to each sequence read is unreliable if the D region is too short owing to extensive trimming. Mora et al. (217) learned from data and analyzed statistical models of the D gene flanked by its junctions. These models are based on the principle of maximum entropy and make minimal assumptions about the mechanisms of diversity – they only rely on the observed frequencies of amino acid pairs along the sequence. These models were used to describe global features of the sequence ensemble, such as the probability distribution following Zipf’s law (220) – the observation that the probability of sequences is inversely proportional to their frequency-rank, or the observation of peaks of frequency in sequence landscape as possible signatures of past pathogenic challenges. Recently, the estimation of repertoire diversity and clonal size distribution were analyzed by Poisson abundance models (221) and simple bivariate-Poisson-lognormal (BPLN) parametric model for fitting and analyzing TR repertoire data was proposed (222). Similarly, network analysis of IG repertoire from Weinstein et al. study revealed the possibility to identify subgroups of individuals on the basis of IG network similarity (223).

The task of characterizing the CDR3 at the nucleotide level is made difficult by the fact that a deterministic assignment of the V-(D)-J recombination process is impossible, because any given sequence can be generated by many possible recombination processes. A previous study proposed a probabilistic model of nucleotide trimming of rearranged TR genes derived from a benchmark data set of TRA and TRG V-(D)-J junctions obtained by comparison to the germline genes in the IMGT® tools (224). Recently a statistical method based on the expectation-maximization algorithm was proposed to circumvent this issue and to extract the statistical properties of junctional diversity accurately from data (124). Applying it to human non-productive DNA sequences gave insight into a universal generation mechanism, reproducible from individual to individual. It was shown that each sequence could potentially be generated by the equivalent of ∼30 equally likely ways by convergent recombination. This method showed that the potential diversity of the recombination machinery was equivalent to ∼1014equally likely sequences (and a practically infinite total number of possible sequences), much more than the estimated 1012 T cells that a single human body can hold. The frequencies of the V, D, and J genes is non-uniform, even at the level of recombination, suggesting underlying physical mechanisms at work. Ndifon et al. (125) proposed a polymer model that accounts for the likelihood of connecting given genomic fragments, giving insight into the mechanistic process.

One of the ultimate goals of deep repertoire sequencing is to find signatures of the repertoire’s response to its antigenic environment. A combination of clustering methods and tree reconstruction techniques have been developed (225, 226) to identify lineages in B cells and study the response to pathogenic challenges. Statistical methods have been devised to detect and quantify the extent of antigen-driven selection acting on B cells, by analyzing the patterns of hypermutations in a Bayesian framework, with applications to deep sequencing data (227, 228).

A lot remains to be done in terms of both data-driven and small-scale models of repertoire-antigen interactions. Ultimately, a close collaboration and development of experimental techniques and models can shed light on how selection at different stages shapes the repertoire, how affinity maturation changes the diversity and the link between sequence diversity and function.

### Future prospects of biomathematical analysis of repertoire data

One of the current challenging issues in antigen-specific repertoire analysis is the development of relevant statistical analysis strategies. Biologists are usually keen on parametric tests, such as ANOVA, t-test, Fischer’s test, among others. However, such statistical methods assume that the inherent probability distribution of the observed variable follows a normal distribution. Rock et al. (229) described that the distribution of the TR diversity is far from following this distribution, thus they proposed the use of non-parametric tests. Nevertheless, different groups are dealing with this issue in order to determine the relevant way to analyze repertoire diversity data and to propose new biostatistics strategies, including principal component analysis, discriminant analysis, hierarchical clustering, specific statistics (164, 169).

In fact, the traditional use of statistics in biology aims at the falsification of a defined hypothesis, i.e., at validating significant differences between defined situations. The recent development of “systems immunology” reverses this point of view and establishes a new usage of multi-parametric statistical approaches to represent the biological data by projections and “landscapes” in the N-dimensional space of considered parameters (230). Thus, the traditional description of separate repertoires for distinct cell subsets defined from a few markers is being replaced by overlapping clouds of data, setting the limits of the different classification groups (tissue of origin, infection contexts, combination of marker expression, repertoire expression…). Moreover, repertoire diversity technologies can now be combined to complementary approaches to decipher the complexity of lymphocyte populations, such as microwell array cell culture and high-resolution imaging (231), mass cytometry (232, 233), cellular barcoding (234), intravital imaging (235, 236), single-cell gene expression (237). In addition, high-throughput repertoire descriptions will enrich mathematical and computer models of lymphocyte repertoire diversity and dynamics such as those proposed by Mehr (238), Ciupe et al. (239), or Stirk et al. (240).