clonotypeR: Identify and analyse B and T cell receptors at a high throughput

Author: Charles Plessy <plessy@riken.jp>
Date:   1 May 2013

clonotypeR is a R package and accompanying scripts to identify and analyse clonotypes from high-throughput T cell receptors sequence libraries. clonotypeR is suited to process and organise very large number of clonotypes, in the order of millions, typically produced by Roche 454 instruments, and to prepare these sequences for differential expression analysis with the typical transcriptomics tools as well as for statistical analysis using existing R packages.

The home page of clonotypeR is http://clonotyper.branchable.com/.

clonotypeR's workflow

Typically, the user receives the output of a next-generation sequencer and runs some shell commands that are not part of the clonotypeR R package, but that are distributed with it on http://clonotyper.branchable.com/.

Currently, clonotypeR provides only support for mouse sequences, by providing pre-formatted sequences from V and J segments. Support of human sequences is in preparation. See http://clonotyper.branchable.com/references/README for more information.

The workflow presented here summarises the different commands to run. Other examples are available on line at http://clonotyper.branchable.com/doc/workflow/.

This example analysis assumes a unix system (Linux, Mac OS, …)

Test data

ClonotypeR come with test data in the R package. While it is not deep enough to discuss about biology, you can use it to familiarise yourself with the commands or test them.

The R package is loaded as usual.

library(clonotypeR)

The data is a table of 120 clonotypes in the extdata folder of the package. The command read_clonotypes will parse it in a data frame. The clonotypes are arbitrarily assigned to three libraries called A, B, and C. The read_clonotypes comments determines at load time if the peptidic sequence has a stop codon or is frame-shifted, and records the information in the unproductive column.

clonotypes <- read_clonotypes(system.file('extdata', 'clonotypes.txt.gz', package = "clonotypeR"))
summary(clonotypes)
##  lib                V            J          score            mapq       
##  A:40   TRAV14-2     : 1   TRAJ31 :35   Min.   : 45.0   Min.   :  0.00  
##  B:40   TRAV14-3     :23   TRAJ49 :10   1st Qu.:215.0   1st Qu.: 31.00  
##  C:40   TRAV14D-3/DV8: 1   TRAJ57 :10   Median :226.0   Median : 37.50  
##         TRAV14N-1    :83   TRAJ34 : 9   Mean   :221.6   Mean   : 36.62  
##         TRAV14N-2    :11   TRAJ27 : 8   3rd Qu.:233.0   3rd Qu.: 40.00  
##         TRAV14N-3    : 1   TRAJ23 : 7   Max.   :247.0   Max.   :218.00  
##                            (Other):41                                   
##      read               dna                qual               pep           
##  Length:120         Length:120         Length:120         Length:120        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  unproductive    ambiguous      
##  Mode :logical   Mode :logical  
##  FALSE:102       FALSE:120      
##  TRUE :18                       
##                                 
##                                 
##                                 
## 

The clonotype_table command counts how many times a given clonotype is found in each library. It can also count simpler features, in particular V and J segments, or any combination of them.

head(clonotype_table(levels(clonotypes$lib), data=clonotypes))
##                                A B C
## TRAV14-2 AASAHYGSSGNKLI TRAJ32 1 0 0
## TRAV14-3 AAADYSNNRLT TRAJ7     1 0 0
## TRAV14-3 AAISNYNQGKLI TRAJ23   1 0 0
## TRAV14-3 AAMNQGGSAKLI TRAJ57   1 1 0
## TRAV14-3 AASAVDYANKMI TRAJ47   1 0 0
## TRAV14-3 AATVDSNYQLI TRAJ33    1 0 0
head(clonotype_table(levels(clonotypes$lib), "V", data=clonotypes))
##                A  B  C
## TRAV14-2       1  0  0
## TRAV14-3       6  7  5
## TRAV14N-1     23 24 24
## TRAV14N-2      2  3  4
## TRAV14N-3      1  0  0
## TRAV14D-3/DV8  0  0  1
head(clonotype_table(levels(clonotypes$lib), "J", data=clonotypes))
##        A  B  C
## TRAJ13 1  1  0
## TRAJ15 1  0  1
## TRAJ23 5  0  1
## TRAJ26 2  2  0
## TRAJ27 2  2  3
## TRAJ31 9 12 12

ClonotypeR provides other functions for further analysis. yassai_identifier calculates a unique identifier using the V, J, peptidic and nucleotidic information, following the work of Yassai et al.

# Unique identifier
head(yassai_identifier(clonotypes))
## [1] "aAYt.22A14N2A49L10"     "aYt.2A14N2A49L10"       "aGRt.34A14N1A40L11"    
## [4] "aADTIt.1221A14N1A27L11" "AEMn.411A14N1A23L11"    "sAHy.11A14-2A32L14"

unique_clonotypes and common_clonotypes are typically used when comparing libraries.

clonotypes <- clonotype_table(levels(clonotypes$lib), data=clonotypes)

# First six clonotypes of library C
head(unique_clonotypes("C", data=clonotypes))
## [1] "TRAV14N-1 AASDTNTGKLT TRAJ27" "TRAV14N-1 AASNSNNRIF TRAJ31" 
## [3] "TRAV14N-2 AAYTGYQNFY TRAJ49"  "TRAV14-3 AASYGSSGNKLI TRAJ32"
## [5] "TRAV14N-1 AALNSNNRIF TRAJ31"  "TRAV14N-1 AAPSNTNKVV TRAJ34"
# Count clonotypes found in library A, and B or C.
length(common_clonotypes(group1="A", group2=c("B","C"), data=clonotypes))
## [1] 5
# Matrix of numbers of common clonotypes
common_clonotypes(data=clonotypes)
##    A  B  C
## A 26  2  3
## B  2 27  4
## C  3  4 27

With deeper data, a typical follow-up would be to identify differentially represented clonotypes between libraries, for instance with the edgeR package, or to calculate distance between libraries, for instance with the vegan package.

Example data (3 x 2,000 reads)

The data provided on-line at http://clonotyper.branchable.com/example_data/ is a sub-sample of three sequence libraries of mouse T cell receptors α (2,000 reads each) made on the 454 Titanium or the 454 junior platforms. The original libraries will be deposited in public databanks after publication in a peer-reviewed journal.

These example libraries are called A, B and C, and are in FASTQ format, with entries like the following (the sequence was truncated for the convenience of the display).

@HKTLYLP01B0MTM
gactGTCCATCTTCCTTTTATCGGACACTGAAGTATGGATATCAGAAGTGCAgggccttcccacgggaacg
+
IIIIIIIIIIIHHFF::::G&gt;IIIGGGIIIIIIIIIGGIIIIIIFEBDCDC&lt;//-5522------

Detection of V segments

Run the command clonotypeR detect A.fastq in the same directory as a copy of the file A.fastq.

The result is stored in a temporary directory called extraction_files, that will be created if it does not already exist.

clonotypeR detect compares the sequences to the reference V segments using BWA, and produces output like the following.

[bsw2_aln] read 2000 sequences/pairs (843395 bp)...
[samopen] SAM header is present: 167 sequences.
[main] Version: 0.6.2-r126
[main] CMD: bwa bwasw -t8 /usr/share/clonotypeR/references/V-C/index A.fastq
[main] Real time: 1.099 sec; CPU: 8.225 sec

This indicates that 2,000 reads have been processed, representing 843,395 base pairs in total. There were 167 reference V segments, and the version number of BWA was 0.6.2-r126. The whole process took less than 10 seconds.

Process the example libraries B and C similarly with the commands clonotypeR detect B.fastq and clonotypeR detect C.fastq.

Extraction of CDR3 regions

Run the command clonotypeR extract A in the same directory as where you ran clonotypeR detect A.fastq. The result is a table stored in a directory called clonotypes, that will be created if it does not already exist.

The output is quite voluminous, and indicates which V / J combinations are being found, like on the following.

TRAV14-3    233
    TRAJ61  0
    TRAJ60  0
    TRAJ59  0
    TRAJ58  1
    TRAJ57  39
    TRAJ56  2
    TRAJ55  0

The format of the table is explained in the manual page of the function read_clonotypes() of the R package.

For each library (A, B and C), one file is available in the clonotypes directory. With BWA 0.6.2-r126, the following numbers of clonotypes are found.

  1072 clonotypes/A.tsv
   924 clonotypes/B.tsv
   689 clonotypes/C.tsv

The files need to be concatenated before analysis in R, with the following command.

find clonotypes/ -name '*tsv' | xargs cat > clonotypes.tsv

A copy of the result file is provided in inst/extdata/clonotypes2.tsv.xz for convenience.

Data analysis in R

Load the clonotypeR library:

library(clonotypeR)

Load the data in a R object called clonotypes:

clonotypes <- read_clonotypes('clonotypes.tsv')

Alternatively, you can load the convenience copy from inst/extdata/clonotypes2.tsv.xz (see above).

clonotypes <- read_clonotypes(system.file("extdata", "clonotypes2.tsv.xz", package = "clonotypeR"))

The command summary(clonotypes) already provides useful information.

summary(clonotypes)
##  lib                  V              J            score            mapq       
##  A:1114   TRAV14-1     :1000   TRAJ31 : 403   Min.   : 44.0   Min.   :  0.00  
##  B: 996   TRAV14-2     : 247   TRAJ23 : 281   1st Qu.:214.0   1st Qu.: 31.00  
##  C: 716   TRAV14-3     : 265   TRAJ22 : 279   Median :225.0   Median : 38.00  
##           TRAV14D-3/DV8: 255   TRAJ37 : 160   Mean   :218.6   Mean   : 63.50  
##           TRAV14N-1    : 624   TRAJ34 : 147   3rd Qu.:232.0   3rd Qu.: 56.75  
##           TRAV14N-2    : 250   TRAJ40 : 112   Max.   :248.0   Max.   :223.00  
##           TRAV14N-3    : 185   (Other):1444                                   
##      read               dna                qual               pep           
##  Length:2826        Length:2826        Length:2826        Length:2826       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  unproductive    ambiguous      
##  Mode :logical   Mode :logical  
##  FALSE:2223      FALSE:2787     
##  TRUE :603       TRUE :39       
##                                 
##                                 
##                                 
## 

Identify unique clonotypes, count their sequences in the libraries A, B and C, and store the result as a table arbitrarily named abc.

abc <- clonotype_table(c('A','B','C'), data=clonotypes)
head(abc)
##                             A B C
## TRAV14-1 AAASSGSWQLI TRAJ22 1 0 0
## TRAV14-1 AACNNRIF TRAJ31    1 0 0
## TRAV14-1 AAGAKLT TRAJ39     3 0 0
## TRAV14-1 AAGGSWQLI TRAJ22   1 0 0
## TRAV14-1 AAGTNTGKLT TRAJ27  1 0 0
## TRAV14-1 AAHDTNAYKVI TRAJ30 1 0 0
summary(abc)
##        A                 B                 C           
##  Min.   : 0.0000   Min.   : 0.0000   Min.   :  0.0000  
##  1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.:  0.0000  
##  Median : 0.0000   Median : 0.0000   Median :  0.0000  
##  Mean   : 0.7669   Mean   : 0.6797   Mean   :  0.5107  
##  3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.:  0.0000  
##  Max.   :18.0000   Max.   :22.0000   Max.   :124.0000

The summary shows that the most frequent clonotype is in C. Using R index vectors, we can see that its CDR3 sequence is AASDSNNRIF and that it was not found in the other libraries.

abc[abc$C == 124, ]
##                             A B   C
## TRAV14N-1 AASDSNNRIF TRAJ31 0 0 124

The clonotype_table function can also produce a count table for and combination of V, CDR3 or J segments.

clonotype_table(c('A','B','C'), "V", data=clonotypes)
##                 A   B   C
## TRAV14-1      238 515   0
## TRAV14-2      134  65   0
## TRAV14-3       82  11 117
## TRAV14D-3/DV8 145  51   4
## TRAV14N-1      78  26 399
## TRAV14N-2      83  63  52
## TRAV14N-3     102  33   2
head(clonotype_table(c('A','B','C'), c("V","J"), data=clonotypes))
##                  A  B C
## TRAV14-1 TRAJ11  1  1 0
## TRAV14-1 TRAJ12  2  2 0
## TRAV14-1 TRAJ13  2  1 0
## TRAV14-1 TRAJ15 11 11 0
## TRAV14-1 TRAJ16  4  3 0
## TRAV14-1 TRAJ18  1 22 0