TCR Analysis

Introduction

scRepertoire is designed to take filter contig outputs from the 10x Genomics Cell Ranger pipeline, process that data to assign clonotype based on two TCR or Ig chains and analyze the clonotype dynamics. The latter can be separated into 1) clonotype-only analysis functions, such as unique clonotypes or clonal space quantification, and 2) interaction with mRNA expression data using Seurat or SingleCellExperiment packages.

Load libraries

Load Seurat Object

1. Load contigs

scRepertoire functions using the filtered_contig_annotations.csv output from the 10x Genomics Cell Ranger. This file is located in the ./outs/ directory of the VDJ alignment folder. To generate a list of contigs to use for scRepertoire:

• load the filtered_contig_annotations.csv for each of the samples.

• make a list in the R environment.

2. Combining Contigs into Clones

contig_list <- list(L1, L2, L3_B, L4_B,L5,L6, L7)
contig_list_with_PBMC <- list(L1, L2, L3_B, L4_B,L5,L6, L7, PBMC)

combined.TCR <- combineTCR(contig_list,
                           samples = c("L1", "L2", "L3_B", "L4_B", "L5", "L6", "L7"),
                           removeNA = FALSE,
                           removeMulti = FALSE,
                           filterMulti = FALSE)

combined.TCR_with_PBMC <- combineTCR(contig_list_with_PBMC, 
                           samples = c("L1", "L2", "L3_B", "L4_B", "L5", "L6", "L7", 
                                       "PBMC"),
                           removeNA = FALSE, 
                           removeMulti = FALSE, 
                           filterMulti = FALSE)

#  exportClones(combined.TCR, 
#               write.file = TRUE,
#               dir = "../TCR_Analysis/",
#              file.name = "clones.csv")
# 
# exportClones(combined.TCR_with_PBMC, 
#              write.file = TRUE,
#              dir = "../TCR_Analysis/",
#              file.name = "clones_with_PBMC.csv")

head(combined.TCR_with_PBMC[[1]])

NA

3. Basic Clonal Visualizations

3.1. clonalQuant

#clonalQuant
 clonalQuant(combined.TCR, 
            cloneCall="strict", 
            chain = "both", 
            scale = TRUE)

clonalQuant(combined.TCR_with_PBMC, 
            cloneCall="strict", 
            chain = "both", 
            scale = TRUE)

clonalQuant(combined.TCR, cloneCall = "gene", group.by = "sample", scale = TRUE)

clonalQuant(combined.TCR_with_PBMC, cloneCall = "gene", group.by = "sample", scale = TRUE)

3.2 clonalAbundance

#clonalAbundance

clonalAbundance(combined.TCR_with_PBMC, 
                cloneCall = "gene",  palette = "Zissou 1",
                scale = FALSE)+ 
  scale_fill_manual(values = hcl.colors(8,"geyser"))

clonalAbundance(combined.TCR_with_PBMC, 
                cloneCall = "gene", palette = "Zissou 1",
                scale = TRUE)

3.3 clonalLength

#clonalLength
clonalLength(combined.TCR_with_PBMC, 
             cloneCall="nt", 
             chain = "both") 

clonalLength(combined.TCR_with_PBMC, 
             cloneCall="aa", 
             chain = "both") 

#TRA
clonalLength(combined.TCR_with_PBMC, 
             cloneCall="nt", 
             chain = "TRA", 
             scale = TRUE) 

clonalLength(combined.TCR_with_PBMC, 
             cloneCall="aa", 
             chain = "TRA", 
             scale = TRUE) 

#TRB
clonalLength(combined.TCR_with_PBMC, 
             cloneCall="nt", 
             chain = "TRB", 
             scale = TRUE) 

clonalLength(combined.TCR_with_PBMC, 
             cloneCall="aa", 
             chain = "TRB", 
             scale = TRUE) 

3.4 clonalCompare

# clonalCompare
clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L1", "L2"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L3_B", "L4_B"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L5", "L6","L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L5","L6"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L5","L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L6","L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L1", "L2", "L3_B", "L4_B", "L5", "L6", "L7", "PBMC"), 
                  cloneCall="aa", 
                  graph = "alluvial")

# clonalCompare(combined.TCR_with_PBMC, 
#               top.clones = 10,
#               highlight.clones = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF"),
#               relabel.clones = TRUE,
#               samples = c("L1", "L2"), 
#               cloneCall="aa", 
#               graph = "alluvial")
# 
# 
# clonalCompare(combined.TCR_with_PBMC, clones = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF"),
#     relabel.clones = TRUE, samples = c("L1", "L2"), cloneCall = "aa", graph = "alluvial")

3.5 clonalScatter

#clonalScatter

clonalScatter(combined.TCR_with_PBMC, 
              cloneCall ="gene", 
              x.axis = "L1", 
              y.axis = "L2",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR_with_PBMC, 
              cloneCall ="gene", 
              x.axis = "L3_B", 
              y.axis = "L4_B",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR_with_PBMC, 
              cloneCall ="gene", 
              x.axis = "L5", 
              y.axis = "L6",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR_with_PBMC, 
              cloneCall ="gene", 
              x.axis = "L5", 
              y.axis = "L7",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR_with_PBMC, 
              cloneCall ="gene", 
              x.axis = "L6", 
              y.axis = "L7",
              dot.size = "total",
              graph = "proportion")

4. Visualizing Clonal Dynamics

4.1 clonalHomeostasis

# Visualizing Clonal Dynamics

# clonalHomeostasis

clonalHomeostasis(combined.TCR_with_PBMC, 
                  cloneCall = "gene")

clonalHomeostasis(combined.TCR_with_PBMC, 
                  cloneCall = "gene",
                  cloneSize = c(Rare = 0.001, Small = 0.01, Medium = 0.1, Large = 0.3, Hyperexpanded =1))

4.2 clonalProportion

# clonalProportion

clonalProportion(combined.TCR_with_PBMC, 
                 cloneCall = "gene") 

clonalProportion(combined.TCR_with_PBMC, 
                 cloneCall = "nt",
                 clonalSplit = c(1, 5, 10, 100, 1000, 10000)) 

5. Summarizing Repertoires

# Summarizing Repertoires

#percentAA

percentAA(combined.TCR, 
          chain = "TRB", 
          aa.length = 20)

Warning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generated

# positionalEntropy
# We can also quantify the level of entropy/diversity across amino acid residues along the cdr3 sequence. positionalEntropy() combines the quantification by residue of percentAA() with the diversity calls in clonalDiversity().
# 
# method
# 
# “shannon” - Shannon Diversity
# “inv.simpson” - Inverse Simpson Diversity
# “norm.entropy” - Normalized Entropy

positionalEntropy(combined.TCR_with_PBMC, 
                  chain = "TRB", 
                  aa.length = 20)

Warning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generatedWarning: invalid factor level, NA generated

# positionalProperty
# Like positionalEntropy(), we can also examine a series of amino acid properties along the cdr3 sequences using positionalProperty(). Important differences from the above function for positionalProperty() is dropping NA values as they would void the mean calculation. positionalProperty() also display a ribbon with the 95% confidence interval surrounding the mean value for the selected properties.
# 
# method
# 
# “Atchley” - Atchley Factors
# “Kidera” - Kidera Factors
# “stScales” - stScales Vectors
# “tScales” - tScales Vectors
# “VHSE” - Vectors of Hydrophobic, Steric, and Electronic properties

# positionalProperty(combined.TCR_with_PBMC[c(1,2)], 
#                   chain = "TRB", 
#                   aa.length = 20, 
#                   method = "Atchley") + 
#   scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])
# 
# positionalProperty(combined.TCR_with_PBMC[c(3,4)], 
#                   chain = "TRB", 
#                   aa.length = 20, 
#                   method = "Atchley") + 
#   scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])
# 
# positionalProperty(combined.TCR_with_PBMC[c(5,6)], 
#                   chain = "TRB", 
#                   aa.length = 20, 
#                   method = "Atchley") + 
#   scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])
# 
# positionalProperty(combined.TCR_with_PBMC[c(5,7)], 
#                   chain = "TRB", 
#                   aa.length = 20, 
#                   method = "Atchley") + 
#   scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])
# 
# 
# positionalProperty(combined.TCR_with_PBMC[c(6,7)], 
#                   chain = "TRB", 
#                   aa.length = 20, 
#                   method = "Atchley") + 
#   scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])
# 
# positionalProperty(combined.TCR_with_PBMC[c(5,6,7)], 
#                   chain = "TRB", 
#                   aa.length = 20, 
#                   method = "Atchley") + 
#   scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])


# vizGenes
# A visualization of the relative usage of genes of the TCR or BCR, using vizGenes(). There is some functional crossover between vizGenes() and two functions below called percentGenes() and percentVJ(). But vizGenes() is more adaptable to allow for comparisons across chains, scaling, etc.

vizGenes(combined.TCR_with_PBMC, 
         x.axis = "TRBV",
         y.axis = NULL,
         plot = "barplot",  
         scale = TRUE)

vizGenes(combined.TCR_with_PBMC, 
        x.axis = "TRAV",
        y.axis = NULL,
        plot = "barplot",  
        scale = TRUE)

vizGenes(combined.TCR[c(1,2,3,4,5,6,7)], 
         x.axis = "TRBV",
         y.axis = "TRBJ",
         plot = "heatmap",  
         scale = TRUE)

vizGenes(combined.TCR[c(1,2,3,4,5,6,7)], 
         x.axis = "TRAV",
         y.axis = "TRAJ",
         plot = "heatmap",  
         scale = TRUE)

# For the P17 patient samples, what if we are interested in chain pairings, we can look at TRBV and TRAV at the same time using them as inputs to x.axis and y.axis.
vizGenes(combined.TCR[c(1,2)], 
         x.axis = "TRBV",
         y.axis = "TRAV",
         plot = "heatmap", 
         scale = FALSE)

# percentGenes
# Quantify the proportion of V or J gene usage with percentGenes(). Like percentAA(), we select the chain of interest and then indicate the gene of interest with the gene parameter. Two major limitations of percentGenes() are, 1) the function quantifies only V or J genes, and 2) the quantification of the genes are limited to all the V or J genes seen across the samples, not all possible V or J genes.

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRB", 
             gene = "Vgene")

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRA", 
             gene = "Vgene")

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRB", 
             gene = "Jgene")

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRA", 
             gene = "Jgene")

# percentVJ
# Quantify the proportion of V and J gene usage with percentVJ(). Like percentGenes(), this function will quantify the percentage of V and J paired together across individual repertoires. The output can be visualized using a heatmap or as input for further dimensional reduction.

percentVJ(combined.TCR_with_PBMC, 
          chain = "TRB")

percentVJ(combined.TCR_with_PBMC, 
          chain = "TRA")

# percentKmer
# Another quantification of the composition of the CDR3 sequence is to define motifs by sliding across the amino acid or nucleotide sequences at set intervals resulting in substrings or kmers.
# 
# motif.length
# 
# Numerical value for the length of the kmer.
# top.motifs
# 
# Display the most variable genes determined by mean absolute deviation (MAD).
percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "aa",
            chain = "TRB", 
            motif.length = 3, 
            top.motifs = 25)

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "aa",
            chain = "TRA", 
            motif.length = 3, 
            top.motifs = 25)

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "nt",
            chain = "TRB", 
            motif.length = 3, 
            top.motifs = 25)

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "nt",
            chain = "TRA", 
            motif.length = 3, 
            top.motifs = 25)

6. Comparing Clonal Diversity and Overlap

# clonalDiversity


clonalDiversity(combined.TCR_with_PBMC, 
                cloneCall = "gene")

#clonalSizeDistribution


clonalSizeDistribution(combined.TCR_with_PBMC, 
                       cloneCall = "aa", 
                       method= "ward.D2")

#clonalOverlap

clonalOverlap(combined.TCR_with_PBMC, 
              cloneCall = "strict", 
              method = "morisita")

clonalOverlap(combined.TCR_with_PBMC, 
              cloneCall = "strict", 
              method = "raw")

7. Combining Clones and Single-Cell Objects

# Combining Clones and Single-Cell Objects

#Getting a sample of a Seurat object
scRep_example <- get(data("All_samples_Merged"))

Warning: data set 'All_samples_Merged' not found

#Define color palette 
colorblind_vector <- hcl.colors(n=7, palette = "inferno", fixup = TRUE)


scRep_example <- combineExpression(combined.TCR, 
                                   scRep_example, 
                                   cloneCall="gene", 
                                   group.by = "sample", 
                                   proportion = FALSE, 
                                   cloneSize=c(Single=1, Small=5, Medium=20, Large=100, Hyperexpanded=500))

DimPlot(scRep_example, group.by = "cloneSize")

DimPlot(scRep_example, group.by = "cloneSize") +
    scale_color_manual(values=rev(colorblind_vector[c(1,3,4,5,6)]))

NA
NA

8. Visualizations for Single-Cell Objects

# Visualizations for Single-Cell Objects
# clonalOverlay

#Adding patient information
scRep_example$Patient_origin <- substr(scRep_example$orig.ident, 1,3)

clonalOverlay(scRep_example, 
              reduction = "umap", 
              freq.cutpoint = 1, 
              bins = 10, 
              facet.by = "Patient_origin") + 
              guides(color = "none")

#clonalNetwork
#ggraph needs to be loaded due to issues with ggplot
library(ggraph)

Attaching package: 'ggraph'

The following object is masked from 'package:sp':

    geometry

clonalNetwork(scRep_example, 
              reduction = "umap", 
              group.by = "seurat_clusters",
              filter.clones = NULL,
              filter.identity = NULL,
              cloneCall = "aa")

#Examining Cluster 3 only
clonalNetwork(scRep_example, 
              reduction = "umap", 
              group.by = "seurat_clusters",
              filter.identity = 8,
              cloneCall = "aa")

shared.clones <- clonalNetwork(scRep_example, 
                               reduction = "umap", 
                               group.by = "seurat_clusters",
                               cloneCall = "aa", 
                               exportClones = TRUE)
head(shared.clones)

scRep_example <- highlightClones(scRep_example, 
                    cloneCall= "aa", 
                    sequence = c("CATGPNGSSNTGKLIF;CALSNNARLMF_CSATTGFYGYTF", 
                                 "CAYSESGGSNYKLTF_CSALAGGYTDTQYF",
                                 "CLVGETGRRALTF_CSARGDRGQPQHF",
                                 "CLVGETGRRALTF_CSARGDRGQPQHF;CATSDFKQVSNQPQHF",
                                 "CATGPNGSSNTGKLIF_CSATTGFYGYTF",
                                 "NA_CSARGDRGQPQHF"))

Seurat::DimPlot(scRep_example, group.by = "highlight") + 
  ggplot2::theme(plot.title = element_blank())

#clonalOccupy
clonalOccupy(scRep_example, 
              x.axis = "seurat_clusters")

clonalOccupy(scRep_example, 
                     x.axis = "ident", 
                     proportion = TRUE, 
                     label = FALSE)

# alluvialClones

# scRep_example$Patient_origin <- substr(scRep_example$orig.ident, 8,8)
# 
# alluvialClones(scRep_example, 
#                cloneCall = "aa", 
#                y.axes = c("Patient_origin", "orig.ident", "Cell_line_Immunophenotype"), 
#                color = c("CATGPNGSSNTGKLIF;CALSNNARLMF_CSATTGFYGYTF", 
#                                  "CAYSESGGSNYKLTF_CSALAGGYTDTQYF",
#                                  "CLVGETGRRALTF_CSARGDRGQPQHF",
#                                  "CLVGETGRRALTF_CSARGDRGQPQHF;CATSDFKQVSNQPQHF",
#                                  "CATGPNGSSNTGKLIF_CSATTGFYGYTF",
#                                  "NA_CSARGDRGQPQHF")) + 
#     scale_fill_manual(values = c("grey", colorblind_vector[3]))
# 
# 
# 
# alluvialClones(scRep_example, 
#                    cloneCall = "gene", 
#                    y.axes = c("Patient_origin", "orig.ident", "Cell_line_Immunophenotype"), 
#                    color = "ident") 



library(circlize)
library(scales)

circles <- getCirclize(scRep_example, 
                       group.by = "seurat_clusters")

#Just assigning the normal colors to each cluster
grid.cols <- hue_pal()(length(unique(scRep_example$seurat_clusters)))
names(grid.cols) <- unique(scRep_example$seurat_clusters)

#Graphing the chord diagram
chordDiagram(circles, self.link = 1, grid.col = grid.cols)

circles <- getCirclize(scRep_example, group.by = "cell_line", proportion = TRUE)

grid.cols <- scales::hue_pal()(length(unique(scRep_example@active.ident)))
names(grid.cols) <- levels(scRep_example@active.ident)

chordDiagram(circles, 
             self.link = 1, 
             grid.col = grid.cols, 
             directional = 1, 
             direction.type =  "arrows",
             link.arr.type = "big.arrow")