What is ECLIPSE? Why should I use it to analyze my scTCR-seq data?

Traditional analysis of single-cell TCR sequencing (scTCR-seq) data defines T cell clones based on the paired CDR3 sequences of the TCRα and TCRβ chains. However, as >50% of cells in some samples possess only 1 or >2 TCR chains, this approach is ill-suited to interpret the clonotypes of a substantial proportion of cells.

ECLIPSE builds on VDJdive to address this limitation of standard scTCR-seq analysis. By analyzing chain pairings across a sample, it is able to resolve clonal ambiguity, allowing for the prediction of chains lost to dropout and the deciphering of true biological expression of extra TCR chains from technical artifacts. With less missing data and larger clone sizes, users are then able to make best use of their scTCR-seq data. The ECLIPSE workflow is also quite streamlined, typically utilizing just 2 core functions (ECLIPSE –> tcrArtifactDetect). After completing these two functions, users are ready for convenient downstream analysis with the scRepertoire package.

What are the inputs to ECLIPSE? Is there any pre-processing required?

ECLIPSE has two main inputs. The first is all or filtered contig files that can be generated by running cellranger vdj or cellranger multi on scTCR-seq FASTQ files. The second is a Seurat object containing the cells in the contig files. Alternatively, users may use SingleCellExperiment objects (see the supplementary information section for details).

How is ECLIPSE different from other Bioconductor and non-Bioconductor scTCR-seq analysis packages?

ECLIPSE has similarities with other scTCR-seq analysis packages, but it also has many differences. Like many other packages, it defines individual TCR chains by the CDR3 amino acid sequence and uses contig files from cellranger for its input. However, in contrast to other packages, for cells that don’t have a clear, unambiguous clonotype (i.e. 1 TCRα and 1 TCRβ), ECLIPSE examines chain pairings from other cells in the same sample to predict the most likely clonotype. Conventional packages (such as scRepertoire) would instead filter the cell out, leave the cell with missing chains, or only select the highest expressed chain of each type.

Another core distinction of ECLIPSE from other packages is it is specifically designed to annotate clones with 3 TCR chains and discern true biological expression of 3 TCR chains from sequencing artifacts such as doublets or TCR contamination. Standard packages either do not allow clones to have 3 TCR chains or do not distinguish true biological expression from sequencing artifacts.

ECLIPSE shares some features with another Bioconductor scTCR-seq analysis package, VDJdive. VDJdive implemented the EM algorithm to predict the clonotypes of cells with ambiguous clonotypes. ECLIPSE uses this implementation; however, there are some key differences between the packages. First, ECLIPSE utilizes custom contig filtering in order to preserve as many TCR chains as possible, while VDJdive uses the filtered contigs from cellranger. ECLIPSE has hard assignment of cells to clonotypes (i.e. each cell is assigned a single clonotype), whereas VDJdive has soft assignment (i.e. 1 cell can be split into multiple different clonotypes). Lastly, ECLIPSE allows clones to possess 2 or 3 TCR chains and integrates with scRepertoire for downstream analysis and visualization, while VDJdive permits clones to have only 2 TCR chains and has its own visualization/analysis functionalities. Other differences exist, but these are the most important ones to understand.

Contig files

The first input is contig files from cellranger multi or cellranger vdj, specifically all_contig_annotations.csv or filtered_contig_annotations.csv files. We strongly recommend using all_contig_annotations files as the filtered contig files remove some biologically true contigs (ex. TCRα chains that use the TRDV1 gene segment, see here) as well as non-productive contigs. While these non-productive chains may not influence T cell antigen specificity, their presence in the contig files allows for enhanced fidelity of clonal tracking by ECLIPSE.

To feed the contig files into ECLIPSE, there are two options:

# Finding where the vignette data for ECLIPSE is stored
vignette_dir <- system.file("extdata/vignette", package = "ECLIPSE")

# Generates the vector of folder paths and excludes the current directory
tcr_folder_paths <- list.dirs(path = vignette_dir, full.names = TRUE)[-1]
print(tcr_folder_paths)
#> [1] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D3B"
#> [2] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D3L"
#> [3] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D5B"
#> [4] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D5L"

# Generates the vector of file paths
contig_paths <- list.files(path = vignette_dir,
                           pattern = "all_contig", 
                           recursive = TRUE, 
                           full.names = TRUE)
contig_paths
#> [1] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D3B/all_contig_annotations.csv.gz"
#> [2] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D3L/all_contig_annotations.csv.gz"
#> [3] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D5B/all_contig_annotations.csv.gz"
#> [4] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D5L/all_contig_annotations.csv.gz"

Seurat object meta.data

4 arguments of ECLIPSE must be provided column names from the Seurat meta.data: original_barcode, batch, donor, and group. We will create these columns now.

# Existing cell barcodes. They aren't the original!
head(rownames(vignette_object[[]]))
#> [1] "Sar3B_ACCAGTAGTAATCACC-1" "Sar3B_GAGGTGAGTGACCAAG-1"
#> [3] "Sar3B_CAGTAACCACCTCGGA-1" "Sar3B_CCTTCGACAGAAGCAC-1"
#> [5] "Sar3B_ACTTGTTCAAGGTTCT-1" "Sar3B_AATCCAGTCGAGAACG-1"

# Making a column with these barcodes where there is no prefix
vignette_object$cell_barcode <- str_replace(rownames(vignette_object[[]]),
                                            pattern = "^.*_", 
                                            replacement = "")
head(vignette_object$cell_barcode) %>% as.vector()
#> [1] "ACCAGTAGTAATCACC-1" "GAGGTGAGTGACCAAG-1" "CAGTAACCACCTCGGA-1"
#> [4] "CCTTCGACAGAAGCAC-1" "ACTTGTTCAAGGTTCT-1" "AATCCAGTCGAGAACG-1"

# Looking at the vector of file_paths
contig_paths
#> [1] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D3B/all_contig_annotations.csv.gz"
#> [2] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D3L/all_contig_annotations.csv.gz"
#> [3] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D5B/all_contig_annotations.csv.gz"
#> [4] "/tmp/Rtmp0ERlbs/Rinst1ea847c3d1f8/ECLIPSE/extdata/vignette/D5L/all_contig_annotations.csv.gz"

# Since normal tissue from sample 15 is first (S15N), it is 1.
# Tumor tissue from 15 is next (S15T), so it is 2, etc...
vignette_object[[]] <- vignette_object[[]] %>% 
                       mutate(index = case_when(
                                                Sample == "D3B" ~ 1,
                                                Sample == "D3L" ~ 2,
                                                Sample == "D5B" ~ 3,
                                                Sample == "D5L" ~ 4))
    
# Checking the column
table(vignette_object$index)
#> 
#>   1   2   3   4 
#> 250 250 250 250

# Show donor column
table(vignette_object$donor)
#> 
#>  D3  D5 
#> 500 500

# Show Sample column
table(vignette_object$Sample)
#> 
#> D3B D3L D5B D5L 
#> 250 250 250 250

Description of function

Now we are ready to run the ECLIPSE function to analyze our scTCR-seq data. Briefly, this function does the following:

1. Reading in contig data
2. Filtering contigs (e.g. must be in Seurat, high_confidence, not duplicated, TCRα or TCRβ, have a valid CDR3, etc…)
3. Calling findThreeChainClones to find clones with biological expression of 3 total TCR chains and modifying the contigs to address this
4. Utilizing the implementation of the EM algorithm in VDJdive’s clonoStats() function. This function analyzes chain pairings across a sample to determine the probability of each potential chain pairing for cells possessing only 1 or >2 total TCR chains
5. Determining the final clonotype (stored in a column called CTaa). If it is highly likely the cell is in 1 clone, then that clone is used. If no clonotype is predicted with high confidence, then clonal prediction is not performed, and the final clonotype is instead determined by scRepertoire::combineTCR() using parameters passed to ECLIPSE from the combineTcrArgs argument (see supplemental information section)
6. Formatting the object for downstream analysis with scRepertoire

For a full list of columns added by ECLIPSE and their description, please refer to the supplemental information section.

Running ECLIPSE

To run ECLIPSE, there are two options, depending on whether we generated a vector of folders containing contigs or contig file paths. If we want to use a vector of folder paths, we can use the following code:

# Option 1: Running with paths of folders containing contigs
object_after_ECLIPSE <- ECLIPSE(contig_type = "all",
                                contig_folder_paths = tcr_folder_paths,
                                seurat_object = vignette_object,
                                batch = "index",
                                donor = "donor",
                                group = "Sample",
                                original_barcode = "cell_barcode",
                                verbose = FALSE,
                                format = "None")

Note that we use the contig_folder_paths argument, omit the contig_file_paths argument, and contig_type = "all". Together, this tells ECLIPSE to search for all_contig_annotations.csv files in the folders that we provide. Additionally, this dataset stores missing data as "None", so we use the format = "None" argument. Newer datasets typically store missing data as "", so this argument should usually be omitted.

If we want to provide TCR contig file paths, we can use the contig_file_paths argument, omit the contig_folder_paths argument, and contig_type = "manual":

# Option 2: Running with contig paths
object_after_ECLIPSE <- ECLIPSE(contig_type = "manual",
                                contig_file_paths = contig_paths,
                                seurat_object = vignette_object,
                                batch = "index",
                                donor = "donor",
                                group = "Sample",
                                original_barcode = "cell_barcode",
                                format = "None")
#> 
#> Reading TCR contig data in...
#> QC Filtering:
#>  - 2883 initial contigs found before filtering
#>  - 2883 contigs were from true cells (coming from 1000 cells)
#>  - 24 contigs were removed for not being high confidence as defined by Cell Ranger
#>  - 98 orphan contigs were found with no observed pairing. Cells containing these contigs are unlikely to have a clonotype predicted
#>  - 554 contigs were removed for other reasons (ex. irregular CDR3, TCRγ/δ contig, duplicated, etc...)
#>  - After all filtering, 2305 contigs from 1000 true cells remained
#> 
#> Finding clones with extra chains...
#> Running donor 1/2 (D3):
#> 2 clones with an extra TCRα chain were detected, comprising 2% of cells (10 total cells)
#> 1 clones with an extra TCRβ chain were detected, comprising 8.4% of cells (42 total cells)
#> 
#> Running donor 2/2 (D5):
#> 2 clones with an extra TCRα chain were detected, comprising 4.4% of cells (22 total cells)
#> 1 clones with an extra TCRβ chain were detected, comprising 1% of cells (5 total cells)
#> 
#> Running VDJdive::clonoStats() for clone EM assignment...
#> Extracting results from VDJdive...
#> Finalizing clonal assignments...

#> 21 cells were detected in clones that appear to be MAIT or iNKT cells
#> 67.6% of cells in the Seurat object were assigned a clonotype using high confidence prediction from VDJdive::clonoStats()
#> 100% of cells in the Seurat object were assigned a clonotype in the end (either from prediction or traditional analysis)
#> These final clonal annotations are stored in CTaa
#> Summary of the top 5 largest clones:
#> # A tibble: 5 × 4
#> # Groups:   CTaa [5]
#>   CTaa                                   Sample Clone_Size Percent_of_Repertoire
#>   <chr>                                  <chr>       <int>                 <dbl>
#> 1 CAMKTSYDKVIF_CASTPGRGWRKLFF:::CASSEGD… D3B            42                  16.8
#> 2 CAMSANSGGYQKVTF:::CAPEAGSGNTPLVF_CASS… D5L            19                   7.6
#> 3 CALTTDSWGKLQF_CSVEWDRGRSETQYF          D3B            10                   4  
#> 4 CAVKRGGTSYGKLTF_CASSGLAGVYNEQFF        D3L            10                   4  
#> 5 CAVGVTGTASKLTF_CASSPTTGMAHEKNIQYF      D5L             8                   3.2

As seen above, by default (verbose = TRUE), ECLIPSE will print some key stats during filtering as well as a QC plot that shows how well clonal prediction went. For this experiment, most cells (shown in orange) received their final clonotype by prediction. Many of these cells also were fully assigned to one clone by the EM algorithm (i.e. y = 1), so overall the analysis went quite well.

Let’s take a look at the final clonotypes that ECLIPSE decided on, now stored in the CTaa column of the Seurat meta.data:

# Final ECLIPSE clonotype
head(object_after_ECLIPSE$CTaa) %>% as.vector()
#> [1] "CAVSLLTGGGNKLTF_CASRGQDHYGYTF"   "CAAHTGTASKLTF_CASSLTIFTQYF"     
#> [3] "CALKGNYKLTF_CASSQPGGGGAYEQYF"    "CAVIDSNYQLIW_CASSYSLAGADTQYF"   
#> [5] "CALPSSGGGGADGLTF_CASSPTGPYQPQHF" "CAAWDSVNDYKLSF_CASSYQGSKETQYF"

We can also compare these final clonotypes to what traditional analysis would have determined, now stored in the scR_CTaa column:

# Clonotype from traditional analysis
head(object_after_ECLIPSE$scR_CTaa) %>% as.vector()
#> [1] "CAVSLLTGGGNKLTF_CASRGQDHYGYTF"   "CAAHTGTASKLTF_CASSLTIFTQYF"     
#> [3] "CALKGNYKLTF_CASSQPGGGGAYEQYF"    "CAVIDSNYQLIW_CASSYSLAGADTQYF"   
#> [5] "CALPSSGGGGADGLTF_CASSPTGPYQPQHF" "CAAWDSVNDYKLSF_CASSYQGSKETQYF"

ECLIPSE also provides clone sizes and clonal proportions in the clonalFrequency and clonalProportion columns, respectively, with matched columns having a scR_prefix to show what would happen with traditional analysis. For a full list of columns added by ECLIPSE, please refer to the supplemental information section.

Using ECLIPSE without a Seurat object

library(SingleCellExperiment)
library(Seurat)

sce_after_ECLIPSE <- sceECLIPSE(sce_object = sce_obj,
                                contig_type = "all",
                                contig_folder_paths = tcr_folder_paths,
                                batch = "index",
                                donor = "donor",
                                group = "Sample",
                                original_barcode = "cell_barcode",
                                verbose = TRUE,
                                format = "None")
#> 
#> Reading TCR contig data in...
#> QC Filtering:
#>  - 2883 initial contigs found before filtering
#>  - 2883 contigs were from true cells (coming from 1000 cells)
#>  - 24 contigs were removed for not being high confidence as defined by Cell Ranger
#>  - 98 orphan contigs were found with no observed pairing. Cells containing these contigs are unlikely to have a clonotype predicted
#>  - 554 contigs were removed for other reasons (ex. irregular CDR3, TCRγ/δ contig, duplicated, etc...)
#>  - After all filtering, 2305 contigs from 1000 true cells remained
#> 
#> Finding clones with extra chains...
#> Running donor 1/2 (D3):
#> 2 clones with an extra TCRα chain were detected, comprising 2% of cells (10 total cells)
#> 1 clones with an extra TCRβ chain were detected, comprising 8.4% of cells (42 total cells)
#> 
#> Running donor 2/2 (D5):
#> 2 clones with an extra TCRα chain were detected, comprising 4.4% of cells (22 total cells)
#> 1 clones with an extra TCRβ chain were detected, comprising 1% of cells (5 total cells)
#> 
#> Running VDJdive::clonoStats() for clone EM assignment...
#> Extracting results from VDJdive...
#> Finalizing clonal assignments...

#> Meta data contains an 'ident' column and will likely result
#>                   in errors downstream.
#> Meta data contains an 'ident' column and will likely result
#>                   in errors downstream.
#> Meta data contains an 'ident' column and will likely result
#>                   in errors downstream.
#> 21 cells were detected in clones that appear to be MAIT or iNKT cells
#> 67.6% of cells in the Seurat object were assigned a clonotype using high confidence prediction from VDJdive::clonoStats()
#> 100% of cells in the Seurat object were assigned a clonotype in the end (either from prediction or traditional analysis)
#> These final clonal annotations are stored in CTaa
#> Summary of the top 5 largest clones:
#> # A tibble: 5 × 4
#> # Groups:   CTaa [5]
#>   CTaa                                   Sample Clone_Size Percent_of_Repertoire
#>   <chr>                                  <chr>       <int>                 <dbl>
#> 1 CAMKTSYDKVIF_CASTPGRGWRKLFF:::CASSEGD… D3B            42                  16.8
#> 2 CAMSANSGGYQKVTF:::CAPEAGSGNTPLVF_CASS… D5L            19                   7.6
#> 3 CALTTDSWGKLQF_CSVEWDRGRSETQYF          D3B            10                   4  
#> 4 CAVKRGGTSYGKLTF_CASSGLAGVYNEQFF        D3L            10                   4  
#> 5 CAVGVTGTASKLTF_CASSPTTGMAHEKNIQYF      D5L             8                   3.2

### Check the original columns in colData
colnames(colData(sce_obj))
#> [1] "orig.ident"   "Cell_barcode" "nCount_RNA"   "donor"        "Sample"      
#> [6] "cell_barcode" "index"        "ident"

### Check that correct columns are added to colData
colnames(colData(sce_after_ECLIPSE))
#>  [1] "orig.ident"                  "Cell_barcode"               
#>  [3] "nCount_RNA"                  "donor"                      
#>  [5] "Sample"                      "cell_barcode"               
#>  [7] "index"                       "ident"                      
#>  [9] "highConfChains"              "scR_CTgene"                 
#> [11] "scR_CTnt"                    "scR_CTaa"                   
#> [13] "scR_CTstrict"                "scR_clonalProportion"       
#> [15] "scR_clonalFrequency"         "scR_cloneSize"              
#> [17] "cloneCallSource"             "CTaa"                       
#> [19] "clonalProportion"            "clonalFrequency"            
#> [21] "cloneSize"                   "cell_unconventional_subset" 
#> [23] "clone_unconventional_subset"

Helper function: findThreeChainClones

While not discussed above, ECLIPSE contains a forth function that may may be of use. This function was designed to help the ECLIPSE function find clones with biological expression of 3 TCR chains, but can be also used to see which clones were called with 3 TCR chains and how that occurred. To see these clones we can use the code below. First, we generate QC filtered contigs using the returnFilteredContigs argument of ECLIPSE:

postqc_contigs <- ECLIPSE(contig_type = "all",
                          contig_folder_paths = tcr_folder_paths,
                          seurat_object = vignette_object,
                          batch = "index",
                          donor = "donor",
                          group = "Sample",
                          original_barcode = "cell_barcode",
                          verbose = FALSE,
                          returnFilteredContigs = TRUE,
                          format = "None")

Next, we filter on a sample we are interested in, and pass the contigs onto findThreeChainClones using mode = "show" and chain = "a" to show clones with 2 TCRα chains:

d3_contigs <- postqc_contigs %>% filter(donor == "D3")

three_chain_clones <- findThreeChainClones(d3_contigs,
                                           chain = "a",
                                           mode = "show",
                                           printStats = TRUE)
#> Returning a list. The first table shows all beta chains predicted to have a double alpha pairing. The second table shows how/why many of these clones were called as such. Please note, the second table is pre-filtering, so not all clones may been in the first table since some may have been filtered out. In this context, the beta chain is the monotypic chain, and the alpha is the bitypic chain.

This provides a list containing:

1. A table of all clones with 2 TCRα chains after filtering

three_chain_clones[[1]]
#> # A tibble: 2 × 2
#>   monotypic_chain predicted_bitypic_pairing      
#>   <chr>           <chr>                          
#> 1 CSVFESDTEAFF    CAESYLNSGNTPLVF;CVVNWEIGFGNVLHC
#> 2 CSVLGQGNEQYF    CAIWGADKLIF;CVVNPVTDSWGKLQF

2. Pre-filtering table that shows how some of these clones were called

three_chain_clones[[2]]
#>   total_cells cell_counts_per_bitypic_pairing monotypic_chain
#> 1           6                           3 2 1    CSVFESDTEAFF
#> 2           4                           2 1 1    CSVLGQGNEQYF
#>         predicted_bitypic_pairing                              B1
#> 1 CAESYLNSGNTPLVF;CVVNWEIGFGNVLHC CAESYLNSGNTPLVF;CVVNWEIGFGNVLHC
#> 2     CAIWGADKLIF;CVVNPVTDSWGKLQF     CAIWGADKLIF;CVVNPVTDSWGKLQF
#>                B2              B3
#> 1 CVVNWEIGFGNVLHC CAESYLNSGNTPLVF
#> 2     CAIWGADKLIF CVVNPVTDSWGKLQF

ECLIPSE additional arguments

Seurat object columns description