VDJtools: Unifying Post-analysis of T Cell Receptor Repertoires

1This is a PLOS Computational Biology Software paper

2The advent of high throughput sequencing (HTS) has opened a new venue for the studies of genomics of adaptive immunity that involve deep profiling of T-cell receptor (TCR) and B-cell receptor (BCR) gene repertoires encoding a myriad of antigen specificities.

3Huge volumes of complex data produced by the immune repertoire profiling have led to the development of a diverse set of software tools, which often complement each other. We and others have recently contributed several tools that handle large amounts of raw HTS data to process it into a human-readable list of clonotypes characterized by Variable (V), Diversity (D), Joining (J) segments and V-(D)-J junction sequences of receptor genes. While such processed data carry nearly exhaustive information on the sampled immune repertoire, this information yet needs to be convolved, scaled and compared across various samples to result in sound biological conclusions.

4Post-analysis of immune repertoire data is a challenging task owing to extreme diversity of TCR and BCR sequences. For example, in technically similar microbiome profiling by 16S rRNA sequencing one deals with thousands of operational taxonomic units that represent various species, while typical TCR repertoire samples may contain hundreds of thousands of clonotypes. Moreover, the species phylogeny and annotation is well developed in the field of microbiology, while immune repertoires remain poorly annotated. To illustrate this, a simple query with “16S rRNA” currently yields more than 8 million records in GenBank, while there are only 37 thousand records annotated as “T-cell receptor”. However, unsupervised methods of studying repertoires, for example based on sample overlap, could turn out very promising, as there exists a relatively limited diversity of overlapping clonotypes.

5In the light of recent advances in storage and processing of immunological big data, community-driven initiatives for immune repertoire data sharing and analysis are likely to emerge, for example VDJserver portal which is currently under development. There are several commonly used ways to survey immune repertoire information obtained from HTS, such as tracking individual clonotypes, comparing immune receptor segment usage and comparing repertoire diversity. Still those are overwhelmingly performed using in-house scripts or even manually. This is becoming a major obstacle, as comparison and annotation of samples based on data generated in other studies is critical for comprehensive analysis of immune repertoire sequencing data. In contrast, similar fields, such as metagenomics, have a plethora of such instruments.

6The VDJtools software package presented here aims at filling this gap by incorporating a comprehensive set of routines for analysis of TCR repertoire sequencing data (Fig 1). The variety of implemented algorithms range from basic statistics calculation and clonotype table filtering to advanced routines such as repertoire clustering and computationally intensive routines such as clonotype table joining.

7VDJtools can calculate basic immune repertoire statistics that were commonly used in pre-HTS era repertoire analysis. Those include in silico spectratype (the distribution of lengths of CDR3 nucleotide sequences) that was first introduced with corresponding molecular biology assay, and various Variable/Joining segment usage statistics that root in flow cytometry analysis of T- and B-cell populations.

8The framework provides means for analyzing the diversity of immune repertoires, such as normalized unique clonotype counts (with an option to account for convergent recombination), clonotype frequency distribution, as well as rarefaction curves and lower bound estimates of total repertoire diversity widely applied in ecology field. The concept of repertoire diversity is of great importance, as it reflects the ability of our immune system to effectively withstand a multitude of encountered pathogens. By applying computational methods one could estimate how the diversity is influenced by various processes, such as aging, vaccination, and infection. Diversity measures could also be used to compare the structure of T- and B-cell repertoires in samples derived from a variety of tissues and subjects.

9Advanced set of VDJtools methods includes cluster analysis of repertoire samples and clonotype tracking which have a wide range of applications. Machine learning methods such as hierarchical clusterization and multi-dimensional scaling can aid in learning T-cell antigen specificities and disease biomarker patterns from high-dimensional TCR data. Clonotype tracking is useful in studying immune repertoire dynamics associated with vaccination, autologous hematopoietic stem cell transplantation (HSCT), checkpoint inhibitors, etc., as well as in detection of minimal residual disease in lymphoid malignancies.

10An overview of 20 recently published immune repertoire studies (S1 Table) demonstrates that VDJtools can perform most of emerging post-analysis tasks therefore greatly facilitating the analysis process and removing the need to develop multiple custom scripts. Currently there are few software tools capable to perform post-analysis of immune repertoire data, all of which provide less functionality when compared to VDJtools (S2 Table). Moreover, in contrast to VDJtools which can handle output generated by various pre-processing software, these tools only support datasets in their internal formats.

11The study was approved by ethics committee of the Russian Children's Hospital from January 20, 2011.

12The core API of the software is implemented in Java/Groovy languages and automatically resolves all dependencies during compilation using Maven. The API includes generalized entities, such as Clonotype, Sample and SampleCollection classes, and allows storing sample metadata using MetadataTable class. The API also contains a comprehensive set of routines for computing sample-specific and cross-sample statistics, which are optimized for parallel computation. VDJtools API can be easily integrated in any software written in Java or related programming languages (e.g. Groovy, Scala and Clojure). VDJtools is an open-source software, the source code can be accessed at GitHub.

13Comprehensive software documentation is hosted at ReadTheDocs and contains basic usage guidelines (including the description of common pitfalls), a summary of implemented algorithms, as well as examples that cover some typical VDJtools usage cases. The documentation also contains step-by-step instructions for reproducing the analysis described in present paper.

14VDJtools has a command line interface that allows executing analysis routines that produce tabular and publication-ready graphical output. Tabular output can be used for post-hoc analysis in R or explored in spreadsheet software such as Excel. Plotting parameters are optimized to provide the most intuitive and comprehensive graphical representation for most usage cases while users can specify their own sample groups and factors to be visualized.

15VDJtools accepts tabular output of commonly used pre-processing software: MIGEC, MiTCR, ImmunoSEQ, IMGT/HighV-QUEST, and MiXCR. VDJtools also supports IgBlast software format. Of note, using IgBlast requires a considerable amount of parsing and post processing, as it only reports Variable segment alignment and doesn’t provide the CDR3 sequence. Moreover, vanilla IgBlast doesn’t accept FASTQ format input, does not provide clonotype assembling (grouping of sequencing reads with identical Variable segment, Joining segment and CDR3 sequence) and is not optimized for parallel computations. We have implemented all those features in our wrapper for IgBlast software, MIGMAP, that could be downloaded from. VDJtools converts all input datasets to its own internal format, which is a tab-delimited table containing abundance, CDR3 sequence, V, D and J segment names and markup of CDR3 sequence germline regions.

16An immune repertoire browser VDJviz which serves as a lightweight GUI for VDJtools was built using Play framework and VDJtools API and could be accessed at.

17Raw data for multiple sclerosis patients is deposited at SRA (PRJNA280417). Pre-processed clonotype tables can be found in a separate GitHub repository, which also contains shell scripts that can be used to reproduce the analysis.

18To demonstrate the efficiency of VDJtools, we have analyzed TCR beta repertoires for the peripheral blood samples of 13 young (6–15 years old) individuals diagnosed with multiple sclerosis (MS1-13), and 6–25 years old control group (C1-11) described in Ref.. The multiple sclerosis dataset was prepared and sequenced using the same protocol as the control one. We have also included a sample from the MS8 patient after hematopoietic stem cell transplant (MS8HSCT). The list of samples is provided in S3 Table.

19To remove quantitative biases and reduce possible impact from PCR and sequencing artifacts, we have utilized unique molecular identifiers. Analysis of raw molecular barcoded data was performed using our MIGEC software. Molecular identifier groups represented by a single read were discarded, and the remaining groups were subjected to cDNA consensus assembly and CDR3 extraction as previously described. Hereafter we will use the term T-cell receptor beta chain cDNA molecules (TRBM) for describing clonotype count units. Note that in these experiments we obtained ~0.5 mln cDNA molecules per ~1–10 mln starting T-cells, so we can assume that each TRBM roughly represents a single T cell.

20We have started our analysis by comparing the repertoire diversity of MS and control samples. To support the diversity measure choice and check for possible biases we have performed a benchmark on previously published T cell immunity aging data and additional ANOVA analysis to identify factors that bias diversity estimates (S1 Text, S1 Fig, S4 Table). We have used common diversity measures: the observed diversity (number of unique clonotypes), Chao and Efron estimates for lower bound on total species diversity, Shannon and Simpson indices, as well as extrapolated Chao estimate.

21The benchmark, in which correlation with a physiological (age) and immune status (naïve T-cell count) factors was compared for various diversity estimates, has shown that best correlation can be achieved when samples are normalized to the same size (TRBM count). Correspondingly, ANOVA analysis suggests a strong sampling-related bias. Accounting for this bias is especially important in present case as the rarefaction curves are far from saturation (Fig 2A). Notably, lower bound estimates of total repertoire diversity that are especially affected by sampling bias were applied in some recent studies for the comparison of TCR repertoire diversity under uneven sample sizes.

22Using Chao1 estimate for normalized datasets that has shown the best performance together with Efron estimate (yet is far simpler to compute) in the aforementioned benchmark, one can discover that MS samples have a significantly lower diversity than the controls (Fig 2B). This suggests a substantial expansion of T-cell clones in peripheral blood of MS patients, an observation previously supported only by local measurements such as Sanger sequencing of individual T-cells and spectratyping assays. As control population is slightly older than MS group one can expect even more profound difference in case exact age matching is achieved for the control group. Still, there is no significant difference for the directly observed sample diversity (S2 Fig), which is likely due to the fact that this estimate doesn’t account for the clonotype frequency distribution in sample and thus is less sensitive.

23As there is currently no study describing an application of cluster analysis to a large set of immune repertoire datasets coming from different individuals, we have performed a benchmark of various clustering strategies using a recently published twins TCR repertoires study. We have tested the ability to distinguish TCR repertoires of identical twins from those of unrelated individuals for several commonly used similarity measures, correlation of overlapping clonotype frequencies (R), geometric mean of total frequencies of overlapping clonotypes (F), normalized number of overlapping clonotypes (D,), Jaccard and Morisita-Horn indices. Only the F similarity measure showed significant difference for both TCR alpha and beta chain datasets (S1 Text, S5 Table and S3 Fig). At the same time, it should be noted that R and D measures also proved to be useful in other experimental setups. For example, R measure accurately separated TCR alpha repertoires for the T cell subsets and tissues, as well as mutant and control mice Treg repertoires.

24We have next used cluster analysis to explore whether TCR beta repertoires of MS patients can be distinguished from healthy controls. As some samples were prepared in parallel with single-end sample barcoding, joined and then co-amplified after Illumina adapter ligation, we first checked for the possibility of cross-sample contamination (S4 and S5 Figs). It turned out that direct clustering of samples with F measure resulted in a strong co-clustering of samples prepared in the same batch. To correct for batch effect, we have selected “amino acid NOT nucleotide” clonotype intersection matching rule, i.e. matching of CDR3 amino acid, but not the nucleotide sequences.

25Hierarchical clustering with F similarity measure and “amino acid NOT nucleotide” clonotype matching rule showed some co-clustering for control but not MS datasets (Fig 3A). Further exploration with multidimensional scaling (MDS) method showed that control repertoires of healthy children are more similar to each other according to F similarity measure, while MS repertoires are all different (Fig 3B and 3C). This result is quite similar to our observations of age-related changes in TCR repertoires (our unpublished data). With aging, expansion of antigen-specific clones moves away native repertoires that are initially more close to each other due to the public clonotypes that are frequently produced in recombination. This is in line with observation of early clonal T-cell expansions in MS children (see “Estimating repertoire diversity” section above). Since those expanding T-cell clones, including potentially autoreactive ones, are predominantly private to an MS-affected person this leads to the decrease of the overlap between MS repertoires according to the clonotype size-weighted F similarity measure.

26Keeping in mind that MS was shown to have a Type I-II TCR repertoire bias, i.e. the same prominent Variable segment is used, yet only limited homology between CDR3 region is present in disease specific T-cells, we have performed hierarchical clustering of Variable segment usage profiles (Fig 3D, note that profiles are weighted by TRBM count). The resulting dendrogram distinguishes MS patients and healthy donors with 91% sensitivity and 77% specificity (P = 0.013, Fisher’s exact test for cluster—group association).

27A post-hoc testing was then performed to find out which Variable segments were more abundant in MS donors than in healthy controls (S6 Table). We have determined that 5 Variable segments had a statistically significant increase in frequency, including TRBV5-6 (1.6-fold, P = 2x10−5) and TRBV5-1 (1.5-fold, P = 5x10−4), which were previously reported to have a genetic association with MS. Of note, TRBV20-1 (1.3-fold, P = 2x10−3) which has also emerged in our results was recently shown to have no genetic association with MS in a Sicilian population carrying null allele. This suggests that the observed TRBV20-1—MS association could be either specific for Russian population or represent an indirect biomarker.

28Further we have compared TCR repertoires of blood samples taken from a single MS patient (MS8) before and after HSCT (see Fig 4). We have first tracked the clonotypes present before HSCT procedure to the post-transplantation repertoire (Fig 4A). The resulting plot clearly shows that pre-transplantation clones greatly expand (from ~25% of TRBMs to 75%) and occupy most of homeostatic space in post-HSCT repertoire. The magnitude of this effect resembles the one we previously observed in an ankylosing spondylitis patient HSCT case and in adult MS autologous HSCT study.

29Another peculiar finding is that a strong shift in Variable segment usage is observed, while no such change is present for J segment usage (Fig 4B). TRBV15 and TRBV7-8 ranking 10 and 5 replaced the top two Variable segments TRBV20-1 and TRBV29-1, while top two Joining segments TRBJ2-7 and TRBJ2-1 remained the same. This could not be attributed to CD4/CD8 balance alone, as there is strong differential Joining segment usage between those two populations. Interestingly, a significant HSCT-induced decrease was observed for TRBV5-6, TRBV5-1, TRBV5-8, TRBV7-6 and TRBV20-1 (P = 0.008, two-tailed paired T-test for log TRBM frequencies) segments that were enriched in MS patients compared to healthy controls (see previous section). The total frequency of those segments dropped from 20% of TRBMs to 14%.

30Finally, we have compared CDR3 regions of MS patients to healthy donors using a set of basic features: the length of Variable and Joining segment germline parts remaining within CDR3 region, and VJ junction (NDN) size. The length of CDR3 segment itself is a potent marker of antigen receptor reactivity. For example, longer CDR3 sequences may be more characteristic for potentially cross- and self-reactive immune receptors, while CDR3 variants with low number of randomly added “N” nucleotides are characteristic for public clonotypes, including variants specific to common pathogens such as EBV and CMV. As our analysis shows, MS patients are characterized by longer VJ junction region (Fig 5A). To check whether it is due to specific segment usage profile we have compared VJ junctions from all clonotypes of normal donors to the ones coming from clonotypes that have one of Variable segments previously shown to be over-expressed in MS patients (Fig 5B). We have found that aforementioned TRBV5-6, TRBV5-1, TRBV5-8, TRBV7-6 and TRBV20-1 are intrinsically characterized by longer VJ inserts. However, there is still a significant difference in VJ junction size between MS patients and controls for this subset of TRBV segments (Fig 5C). These results may indicate that clonal expansions in MS patients are characterized by more self-reactive T-cell clonotypes than in healthy donors. Alternatively, this could be a more general hallmark of chronic inflammation associated with MS.

31A cross-platform binary version of software in a form of executable JAR file is available from. VDJtools software is free for scientific and non-profit use. The source code is available at GitHub repositories and.

32One important aspect of VDJtools usage not mentioned in the results section is the benchmark of pre-processing software (S6 Fig) and library preparation protocols. For this purposes we plan to constantly update VDJtools so it is able to handle the output of newly developed pre-processing software.

33In future we plan extending VDJtools software to address another highly important problem in the field, the analysis of antibody repertoire. While being applicable to the analysis of BCR clonotypes, VDJtools currently doesn’t account for somatic hypermutations and therefore yet cannot offer a comprehensive analysis for the antibody repertoires. This task requires us to implement algorithms for computing statistics of hypermutation transition patterns and reconstruction of B-cell clonal lineages and visualization of hypermutation graphs. We are also looking forward for the feedback from the community to meet the demand for some exciting novel features that will surely arise in this rapidly growing field.