Starting work with scRepertoire.
April 8, 2021
Package
scRepertoire 1.4.0
1 Introduction
scRepertoire is designed to take filter contig outputs from the 10x Genomics Cell Ranger pipeline, processes 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, SingleCellExperiment or Monocle 3 packages.
1.1 Loading Libraries
suppressMessages(library(scRepertoire))1.2 Loading and Processing Contig Data
scRepertoire comes with a data set derived from T cells derived from three patients with renal clear cell carcinoma in order to demonstrate the functionality of the R package. More information on the data set can be found at preprint 1 and preprint 2. The samples consist of paired peripheral-blood and tumor-infiltrating runs, effectively creating 6 distinct runs for T cell receptor (TCR) enrichment. We can preview the elements in the list by using the head function and looking at the first contig annotation. Here notice the barcode is labeled as PX_P_############# - this refers to Patient X (PX) and Peripheral Blood (P).
The data, contig_list, is 6 filtered_contig.csv file outputs from Cell Ranger that made into a list.
data("contig_list") #the data built into scRepertoire
head(contig_list[[1]])## barcode is_cell contig_id high_confidence length
## 1 AAACCTGAGAGCTGGT TRUE AAACCTGAGAGCTGGT-1_contig_1 TRUE 705
## 2 AAACCTGAGAGCTGGT TRUE AAACCTGAGAGCTGGT-1_contig_2 TRUE 502
## 3 AAACCTGAGCATCATC TRUE AAACCTGAGCATCATC-1_contig_1 TRUE 693
## 4 AAACCTGAGCATCATC TRUE AAACCTGAGCATCATC-1_contig_2 TRUE 567
## 5 AAACCTGAGCATCATC TRUE AAACCTGAGCATCATC-1_contig_5 TRUE 361
## 6 AAACCTGAGTGGTCCC TRUE AAACCTGAGTGGTCCC-1_contig_1 TRUE 593
## chain v_gene d_gene j_gene c_gene full_length productive cdr3
## 1 TRB TRBV20-1 TRBD1 TRBJ1-5 TRBC1 TRUE TRUE CSASMGPVVSNQPQHF
## 2 TRB None None TRBJ1-5 TRBC1 FALSE None None
## 3 TRB TRBV5-1 TRBD2 TRBJ2-2 TRBC2 TRUE TRUE CASSWSGAGDGELFF
## 4 TRA TRAV12-1 None TRAJ37 TRAC TRUE TRUE CVVNDEGSSNTGKLIF
## 5 TRB None None TRBJ1-5 TRBC1 FALSE None None
## 6 TRB TRBV7-9 TRBD1 TRBJ2-5 TRBC2 TRUE TRUE CASSPSEGGRQETQYF
## cdr3_nt reads umis raw_clonotype_id
## 1 TGCAGTGCTAGCATGGGACCGGTAGTGAGCAATCAGCCCCAGCATTTT 16718 6 clonotype96
## 2 None 6706 3 clonotype96
## 3 TGCGCCAGCAGCTGGTCAGGAGCGGGAGACGGGGAGCTGTTTTTT 26719 11 clonotype97
## 4 TGTGTGGTGAACGATGAAGGCTCTAGCAACACAGGCAAACTAATCTTT 18297 6 clonotype97
## 5 None 882 4 clonotype97
## 6 TGTGCCAGCAGCCCCTCCGAAGGGGGGAGACAAGAGACCCAGTACTTC 11218 6 clonotype98
## raw_consensus_id
## 1 clonotype96_consensus_1
## 2 None
## 3 clonotype97_consensus_2
## 4 clonotype97_consensus_1
## 5 None
## 6 clonotype98_consensus_1Some workflows will have the additional labeling of the standard barcode. Before we proceed, we will use the function stripBarcode() in order to avoid any labeling issues down the line. Importantly, stripBarcode() is for removing prefixes on barcodes that have resulted from other pipelines.
No need for stripBarcode function, if the barcodes look like: + AAACGGGAGATGGCGT-1 + AAACGGGAGATGGCGT
In terms of using stripBarcode(), please think about the following parameters.
column + The column in which the barcodes are present
connector + The character that is connecting the barcode with the prefix
num_connects + the levels of barcode prefix, where X_X_AAACGGGAGATGGCGT-1 == 3, X_AAACGGGAGATGGCGT-1 = 2.
for (i in seq_along(contig_list)) {
contig_list[[i]] <- stripBarcode(contig_list[[i]], column = 1, connector = "_", num_connects = 3)
}You can see now the barcode in column 1, we have removed the P## prefixes.
2 Combining the Contigs
As the output of CellRanger are quantifications of both the TCRA and TCRB chains, the next step is to create a single list object with the TCR gene and CDR3 sequences by cell barcode. This is performed using the combineTCR(), where the input is the stripped contig_list. There is also the relabeling of the barcodes by sample and ID information to prevent duplicates.
cells + T-AB - T cells, alpha-beta TCR + T-GD - T cells, gamma-delta TCR
removeNA + TRUE - this is a stringent filter to remove any cell barcode with an NA value in at least one of the chains + FALSE - the default setting to include and incorporate cells with 1 NA value
removeMulti + TRUE - this is a stringent filter to remove any cell barcode with more than 2 immune receptor chains + FALSE - the default setting to include and incorporate cells with > 2 chains
filterMulti + TRUE - Isolated the top 2 expressed chains in cell barcodes with multiple chains + FALSE - the default setting to include and incorporate cells with > 2 chains
combined <- combineTCR(contig_list,
samples = c("PY", "PY", "PX", "PX", "PZ","PZ"),
ID = c("P", "T", "P", "T", "P", "T"), cells ="T-AB")The output of combineTCR() will be a list of contig data frames that will be reduced to the reads associated with a single cell barcode. It will also combine the multiple reads into clonotype calls by either the nucleotide sequence (CTnt), amino acid sequence (CTaa), the gene sequence (CTgene) or the combination of the nucleotide and gene sequence (CTstrict). The analogous function for B cells, combineBCR() functions similarly with 2 major caveats: 1) Each barcode can only have a maximum of 2 sequences, if greater exists, the 2 with the highest reads are selected. 2) The strict definition of clonotype (CTstrict) is based on the v gene and >85% normalized hamming distance of the nucleotide sequence. The hamming distance is calculated across all BCR sequences recovered, regardless of the run.
3 Other Processing Functions
3.1 Adding Additional Variables
What if there are more variables to add than just sample and ID? We can add them by using the addVariable() function. All we need is the name of the variable you’d like to add and the specific character or numeric values (variables). As an example, here we add the batches in which the samples were processed and sequenced.
example <- addVariable(combined, name = "batch",
variables = c("b1", "b1", "b2", "b2", "b2", "b2"))
example[[1]][1:5,ncol(example[[1]])] # This is showing the first 5 values of the new column added## [1] "b1" "b1" "b1" "b1" "b1"3.2 Subsetting Contigs
Likewise we can remove specific list elements after combineTCR() using the subsetContig() function. In order to subset, we need to identify the vector we would like to use for subsetting (name) and also the variable values to subset (variables). Below you can see us isolate just the 4 sequencing results from PX and PY.
subset <- subsetContig(combined, name = "sample",
variables = c("PX", "PY"))4 Visualizing Contigs
cloneCall + “gene” - use the genes comprising the TCR/Ig + “nt” - use the nucleotide sequence of the CDR3 region + “aa” - use the amino acid sequence of the CDR3 region + “gene+nt” - use the genes comprising the TCR/Ig + the nucleotide sequence of the CDR3 region. This is the proper definition of clonotype.
Important to note, that the clonotype is called using essentially the combination of genes or nt/aa CDR3 sequences for both loci. As of this implementation of scRepertoire, clonotype calling is not incorporating small variations within the CDR3 sequences. As such the gene approach will be the most sensitive, while the use of nt or aa moderately so, and the most specific for clonotypes being gene+nt. Additionally, the clonotype call is trying to incorporate both loci, i.e, both TCRA and TCRB chains and if a single cell barcode has multiple sequences identified (i.e., 2 TCRA chains expressed in one cell). Using the 10x approach, there is a subset of barcodes that only return one of the immune receptor chains, the unreturned chain is assigned an NA value.
The first function to explore the clonotypes is quantContig() to return the total or relative numbers of unique clonotypes. scale + TRUE - relative percent of unique clonotypes scaled by total size of the size of the clonotype repertoire + FALSE - Report the total number of unique clonotypes
quantContig(combined, cloneCall="gene+nt", scale = TRUE)
Within each of the general analysis functions, there is the ability to export the data frame used to create the visualization. To get the exported values, use exportTable = TRUE. It will return the data frame used to make the graph, instead of the visual output.
quantContig_output <- quantContig(combined, cloneCall="gene+nt",
scale = TRUE, exportTable = TRUE)
quantContig_output## contigs values total scaled
## 1 2690 PY_P 3208 83.85287
## 2 1495 PY_T 3119 47.93203
## 3 823 PX_P 1068 77.05993
## 4 918 PX_T 1678 54.70799
## 5 1143 PZ_P 1434 79.70711
## 6 749 PZ_T 2768 27.05925The other option here is to be able to define the visualization by data classes. Here we used the combineTCR() to define the ID variable as part of the naming structure. We can the group to specifically use a column in the data set to organize the visualization.
quantContig(combined, cloneCall="gene", group = "ID", scale = TRUE)
We can also examine the relative distribution of clonotypes by abundance. Here abundanceContig() will produce a line graph with a total number of clonotypes by the number of instances within the sample or run. Like above, we can also group this by vectors within the contig object using the group variable in the function
abundanceContig(combined, cloneCall = "gene", scale = FALSE)
abundanceContig(combined, cloneCall = "gene", group = "ID", scale = FALSE)
As you can see the peripheral blood sample derived from patient 1 is a relative extreme outlier. Another method to examine the relative abundance is to look at the density by using the scale call in the function.
abundanceContig(combined, group = "ID", scale = TRUE)
Lastly on the basic visualization side, we can look at the length distribution of the CDR3 sequences by calling the lengtheContig() function. Importantly, unlike the other basic visualizations, the cloneCall can only be “nt” or “aa”. Due to the method of calling clonotypes as outlined above, the length should reveal a multimodal curve, this is a product of using the NA for the unreturned chain sequence and multiple chains within a single barcode.
lengthContig(combined, cloneCall="aa", chains = "combined") 
Or we can visualize the individual chains of the immune receptors by selecting chains = “single”. Notably this will remove the NA component of combined clonotypes, so visualize is only the sequences recovered in the filtered contig annotation file from Cell Ranger.
lengthContig(combined, cloneCall="nt", chains = "single") 
We can also look at clonotypes between samples and changes in dynamics by using the compareClonotypes() function.
samples + Can be used to isolate specific samples based on the name of the list element
graph + “alluvial” - graph imaged below + “area” - graph by area of the respective clonotype
number + The top number of clonotypes to graph, this will be calculated based on the frequency of the individual sample. This can also be left blank.
clonotypes + Can be used to isolate specific clonotype sequences, ensure the call matches the sequences you would like to visualize.
compareClonotypes(combined, numbers = 10, samples = c("PX_P", "PX_T"),
cloneCall="aa", graph = "alluvial")
5 More Advanced Clonal Analysis
After we have completed the basic processing and summary functions in scRepertoire, we can begin to explore the clonotypes of the single-cell data in more detail.
5.1 Clonal Space Homeostasis
By examining the clonal space, we are effectively looking at the relative space occupied by clones at specific proportions. Another way to think about this would be thinking of the total immune receptor sequencing run as a measuring cup. In this cup, we will fill liquids of different viscosity - or different number of clonal proportions. Clonal space homeostasis is asking what percentage of the cup is filled by clones in distinct proportions (or liquids of different viscosity, to extend the analogy). The proportional cut points are set under the cloneType variable in the function and can be adjusted, at baseline the bins are as follows:
cloneTypes + Rare = .0001 + Small = .001 + Medium = .01 + Large = .1 + Hyperexpanded = 1
clonalHomeostasis(combined, cloneCall = "gene")
clonalHomeostasis(combined, cloneCall = "aa")
5.2 Clonal Proportion
Like clonal space homeostasis above, clonal proportion acts to place clones into separate bins. The key difference is instead of looking at the relative proportion of the clone to the total, the clonalProportion() function will rank the clones by total number and place them into bins.
The split represents ranking of clonotypes by copy or frequency of occurrence, meaning 1:10 are the top 10 clonotypes in each sample. The default bins are under the split variable in the function and can be adjusted, but at baseline they are as follows.
split + 10 + 100 + 1000 + 10000 + 30000 + 100000
clonalProportion(combined, cloneCall = "gene") 
clonalProportion(combined, cloneCall = "nt") 
5.3 Overlap Analysis
If you are interested in measures of similarity between the samples loaded into scRepertoire, using clonalOverlap() can assist in the visualization. Two methods currently can be performed in clonalOverlap() 1) overlap coefficient and 2) Morisita index. The former is looking at the overlap of clonotypes scaled to the length of unique clonotypes in the smaller sample. The Morisita index is more complex, it is an ecological measure of the dispersion of individuals within a population, incorporating the size of the population.
clonalOverlap(combined, cloneCall = "gene+nt", method = "morisita")
Another recent addition to scRepertoire is the ability to cluster the samples by the clone size distribution using clonesizeDistribution() adapted from the powerTCR R package. Please read and cite the respective citation if using this function to analyze the similarities of sample clone size distributions. In this function, method refers to the method for which the hierarchical clustering will be based on.
clonesizeDistribution(combined, cloneCall = "gene+nt",
method="ward.D2")
## NULL5.4 Diversity Analysis
Diversity can also be measured for samples or by other variables. Diversity is calculated using four metrics: 1) Shannon, 2) inverse Simpson, 3) Chao1, and 4) Abundance-based Coverage Estimator (ACE). With the former two generally used to estimate baseline diversity and Chao/ACE indices used to estimate the richness of the samples. New implementation of this function include downsampling with 100 boot straps (n.boots) using the minimum number of unique clonotypes, as a more robust diversity estimate.
clonalDiversity(combined, cloneCall = "gene", group = "samples", n.boots = 100)
clonalDiversity(combined, cloneCall = "gene", group = "ID")
5.5 Clustering Clonotypes
The nucleotide or amino acid sequences of the chains can be used to cluster clonotypes, by examining the edit distance of the sequences. This approach is underlying the combineBCR() function, but now can be applied to the T cell receptors. It will add to the end of each list element a cluster by generating a network connected by the similarity in sequence. This network is directed by the threshold variable, where 0.85 is the normalized mean edit distance.
Clonotype Cluster by distance:
- TCRA:LD.1.TRAV29DV5
- [chain:] + :LD + [number] + [vchain]
LD denotes if the cluster was called using the normalized Levenshtein distance
sub_combined <- clusterTCR(combined[[1]], chain = "TCRA",
sequence = "aa", threshold = 0.85)
sub_combined <- as.data.frame(sub_combined)
counts_TCRA <- table(sub_combined$TCR1)
counts_TCRAcluster <- table(sub_combined$TCRA_cluster)
plot(counts_TCRA, axes = FALSE)
plot(counts_TCRAcluster, axes=FALSE)
- TCRA:LD.1.TRAV29DV5
- [chain:] + :LD + [number] + [vchain]
LD denotes if the cluster was called using Levenshtein distance
sub_combined <- clusterTCR(combined[[1]], chain = "TCRA", sequence = "aa", threshold = 0.85)
sub_combined[[1]]$TCRA_cluster[1:20]## [1] NA "TCRA.1" NA "TCRA:LD.1" "TCRA:LD.2" "TCRA:LD.3"
## [7] "TCRA.2" "TCRA:LD.4" "TCRA.3" NA NA NA
## [13] NA "TCRA.4" NA "TCRA.5" "TCRA:LD.5" "TCRA:LD.6"
## [19] "TCRA.4" "TCRA:LD.7"6 Interacting with mRNA Expression
As mentioned previously, this data set is derived from work performed in the laboratory of Weizhou Zhang. For the purposes of the vignette, we have randomly sampled 100 cells from the the fully integrated Seurat object to minimize the size of the package. A full version of the Seurat object is available via GitHub under the Getting Data header. We will use both Seurat and SingleCellExperiment (SCE) with scater to perform the further visualizations in tandem.
library(Seurat)
library(scater)
screp_example <- get(data("screp_example"))
sce <- suppressMessages(UpdateSeuratObject(screp_example))
sce <- as.SingleCellExperiment(screp_example)
#Seurat Format
DimPlot(screp_example)
##Single Cell Experiment Format
plotUMAP(sce, colour_by = "seurat_clusters")
Here you can see we have 12 total clusters (C1-12), which we have labeled as such for simplicity. We can also get a little more granular information on the number of cells by using the table() function.
table(screp_example$seurat_clusters)##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 15 20 8 11 11 12 11 4 4 2 2 0Next we can take the clonotypic information and attach it to our Seurat object using the combineExpression() function. Importantly, the major requirement for the attachment is matching contig cell barcodes and barcodes in the row names of the meta data of the Seurat or SCE object. If these do not match, the attachment will fail. Based on ease, we suggest you make the changes to the Seurat object row names.
We can call (cloneCall) the 4 variations of clonotypes: 1) genes, 2) CDR3 amino acid sequence, 3) CDR3 nucleotide sequence, or 4) genes and CDR3 nucleotide sequence. The attaching function will also calculate the frequency of the clonotype based on the groupBy variable. If blank, groupBy will calculate frequencies of clonotypes by individual run, but because we have 6 samples of paired peripheral and tumor T cells, we are actually going to use the groupBy variable to call “sample” in order to calculate frequencies across both the peripheral blood and tumor T cells of the same patient.
In order to categorize the frequency, we have the variable proportion which if TRUE allows for the relative proportion or when FALSE will use absolute frequency to define clonotype groups cloneTypes acts as a bin to place labels. As a default, cloneTypes is set to equal cloneTypes=c(Rare = 1e-4, Small = 0.001, Medium = 0.01, Large = 0.1, Hyperexpanded = 1). However, below you can see an example of using total frequency as expansion assignements.
screp_example <- combineExpression(combined, screp_example, cloneCall="gene", groupBy = "sample", proportion = FALSE, cloneTypes=c(Single=1, Small=5, Medium=20, Large=100, Hyperexpanded=500))
sce <- combineExpression(combined, sce, cloneCall = "gene", groupBy = "sample")We first want to look at the distribution of peripheral versus tumor T cells. We can use the same color scheme as the rest of the scRepertoire package by calling the object colorblind_vector using the following hex codes.
colorblind_vector <- colorRampPalette(c("#FF4B20", "#FFB433", "#C6FDEC", "#7AC5FF", "#0348A6"))
DimPlot(screp_example, group.by = "Type") + NoLegend() +
scale_color_manual(values=colorblind_vector(2))
We can also look at the composition of each cluster by comparing the proportion of the cluster comprising peripheral blood versus tumor T cells. We can do this by first forming a table of the cluster and type of cells, then scaling the rows of the table by the total number of cells sequenced.
table <- table(screp_example$Type, Idents(screp_example))
table[1,] <- table[1,]/sum(table[1,]) #Scaling by the total number of peripheral T cells
table[2,] <- table[2,]/sum(table[2,]) #Scaling by the total number of tumor T cells
table <- as.data.frame(table)
table$Var2 <- factor(table$Var2,
levels = c("C1", "C2", "C3", "C4", "C5", "C6",
"C7", "C8", "C9", "C10", "C11", "C12"))
ggplot(table, aes(x=Var2, y=Freq, fill=Var1)) +
geom_bar(stat="identity", position="fill", color="black", lwd=0.25) +
theme(axis.title.x = element_blank()) +
scale_fill_manual(values = c("#FF4B20","#0348A6")) +
theme_classic() +
theme(axis.title = element_blank()) +
guides(fill=FALSE)
Now we can look at the distribution of the clonotype bins by first ordering the clonoType as a factor, this prevents the coloring from being in alphabetical order. Next we use the DimPlot() function call in Seurat with our scale_color_manual additional layer.
slot(screp_example, "meta.data")$cloneType <- factor(slot(screp_example, "meta.data")$cloneType,
levels = c("Hyperexpanded (100 < X <= 500)", "Large (20 < X <= 100)",
"Medium (5 < X <= 20)", "Small (1 < X <= 5)",
"Single (0 < X <= 1)", NA))
DimPlot(screp_example, group.by = "cloneType") 
scale_color_manual(values = c(rev(colorblind_vector(5))), na.value="grey")## <ggproto object: Class ScaleDiscrete, Scale, gg>
## aesthetics: colour
## axis_order: function
## break_info: function
## break_positions: function
## breaks: waiver
## call: call
## clone: function
## dimension: function
## drop: TRUE
## expand: waiver
## get_breaks: function
## get_breaks_minor: function
## get_labels: function
## get_limits: function
## guide: legend
## is_discrete: function
## is_empty: function
## labels: waiver
## limits: NULL
## make_sec_title: function
## make_title: function
## map: function
## map_df: function
## n.breaks.cache: NULL
## na.translate: TRUE
## na.value: grey
## name: waiver
## palette: function
## palette.cache: NULL
## position: left
## range: <ggproto object: Class RangeDiscrete, Range, gg>
## range: NULL
## reset: function
## train: function
## super: <ggproto object: Class RangeDiscrete, Range, gg>
## rescale: function
## reset: function
## scale_name: manual
## train: function
## train_df: function
## transform: function
## transform_df: function
## super: <ggproto object: Class ScaleDiscrete, Scale, gg>plotUMAP(sce, colour_by = "cloneType")
6.1 clonalOverlay
Using the dimensional reduction graphs as a reference, we can also generate an overlay of the position of clonally expanded cells using clonalOverlay(). Select the reduction for the visualization, default is “PCA” and the freq.cutpoint or lowest clonal frequency or proportion to generate the contour plot. We can modify the contours by selecting the number of bins or number of contours drawn. clonalOverlay() can be used to look across all cells or faceted by a meta data variable using facet. As we facet, the overall dimensional reduction will be maintained, while the contour plots will adjust based on the facet variable. Coloring of the dot plot is based on the active identity of the single-cell object. This visualization was authored by Dr. Francesco Mazziotta from Johns Hopkins and inspired by Drs. Carmona and Andreatta and their work with ProjectTIL.
clonalOverlay(screp_example, reduction = "umap",
freq.cutpoint = 30, bins = 10, facet = "Patient") +
guides(color = FALSE)
We can also look at the clonotypes by calling specific sequences in the highlightClonotypes() below. In order to highlight the clonotypes, we first need to use the cloneCall the type of sequence we will be using and then the specific sequences themselves using sequence. Below you can see the steps to highlight the two most prominent sequences “CAVNGGSQGNLIF_CSAEREDTDTQYF” with a frequency = 482 (clonotype 1) and “NA_CATSATLRVVAEKLFF” with a frequency = 287 (Clonotype2).
screp_example <- highlightClonotypes(screp_example, cloneCall= "aa",
sequence = c("CAVNGGSQGNLIF_CSAEREDTDTQYF", "NA_CATSATLRVVAEKLFF"))
Seurat::DimPlot(screp_example, group.by = "highlight")
6.2 occupiedscRepertoire
We can also look at the count of cells by cluster assigned into specific frequency ranges, by using the occupiedscRepertoire() function and selecting the x.axis to display cluster or other variables in the meta data of the single cell object.
occupiedscRepertoire(screp_example, x.axis = "cluster")
6.3 alluvialClonotypes
After the metadata has been modified, we can look at clonotypes across multiple categories using the alluvialClonotypes() function. To understand the basic concepts of this graphing method, I’d highly recommend reading this post, essentially we are able to use the plots to examine the interchange of categorical variables. Because this function will produce a graph with each clonotype arranged by called stratifications, this will take some time depending on the size of your total cells. To expedite, we will actually subset the seurat object before using alluvialClonotypes().
alluvialClonotypes(screp_example, cloneCall = "gene",
y.axes = c("Patient", "cluster", "Type"),
color = "TRAV12-2.TRAJ42.TRAC_TRBV20-1.TRBJ2-3.TRBD2.TRBC2") +
scale_fill_manual(values = c("grey", colorblind_vector(1)))
alluvialClonotypes(sce, cloneCall = "gene",
y.axes = c("Patient", "seurat_clusters", "Type"),
color = "seurat_clusters") 
6.4 getCirclize
Like alluvial graphs, we can also visualize the interconnection of clusters using the chord diagrams from the circlize R package. The first step is getting the data frame output to feed into the chordDiagram() function in circlize, which can be done using getCirclize(). This will calculate the relative number of clonotypes shared based on the groupBy variable using the product of combineExpression().
library(circlize)
library(scales)
circles <- getCirclize(screp_example, groupBy = "cluster")
#Just assigning the normal colors to each cluster
grid.cols <- hue_pal()(length(unique(Idents(screp_example))))
names(grid.cols) <- levels(slot(screp_example, "active.ident"))
#Graphing the chord diagram
chordDiagram(circles, self.link = 1, grid.col = grid.cols)
6.5 Diversity of single-cells using Startrac
From the excellent work by Lei Zhang, et al in Lineage tracking reveals dynamic relationships of T cells in colorectal cancer, the authors introduce new methods for looking at clonotypes by cellular origins and cluster identification. Their startrac software has been incorporated into scRepertoire. If you are using the output of this specific function, please cite their excellent work.
In order to use the StartracDiversity() function, you will need to include the product of the combinedSeurat() function. The second requirement is the a column header in the meta data of the Seurat object that has tissue of origin. In the example data,type corresponds to the column “Type”, which includes the “P” and “T” classifier. The indices can be subseted for a specific patient or examined overall using the by variable. Importantly, the function uses only the strict definition of clonotype of the genes and the CDR3 nucleotide sequence.
The indices output includes:
+ expa - Clonal Expansion
+ migr - Cross-tissue Migration
+ tran - State Transition
StartracDiversity(screp_example, type = "Type", sample = "Patient", by = "overall")## [2021-10-26 18:46:29] initialize Startrac ...
## [2021-10-26 18:46:29] calculate startrac index ...
## [2021-10-26 18:46:29] calculate pairwise index ...
## [2021-10-26 18:46:29] calculate indices of each patient ...
## [2021-10-26 18:46:29] collect result
## [2021-10-26 18:46:29] return
For users that would like greater ability to use the meta data in the Seurat objects to perform the analysis that scRepertoire provides, there is also the option of using the expression2List() function that will take the meta data and output the data as a list by cluster.
combined2 <- expression2List(screp_example, group = "cluster")
combined3 <- expression2List(sce, group = "cluster")6.6 1. Clonal Diversity
clonalDiversity(combined2, cloneCall = "nt")
clonalDiversity(combined3, cloneCall = "nt")
6.7 2. Clonal Homeostasis
clonalHomeostasis(combined2, cloneCall = "nt")
clonalHomeostasis(combined3, cloneCall = "nt")
6.8 3. Clonal Proportion
clonalProportion(combined2, cloneCall = "nt")
clonalProportion(combined3, cloneCall = "nt")
6.9 4. Clonal Overlap
clonalOverlap(combined2, cloneCall="aa", method="overlap")
clonalOverlap(combined3, cloneCall="aa", method="overlap")
7 Conclusion
If you have any questions, comments or suggestions, feel free to visit the github repository or email me.
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scales_1.1.1 circlize_0.4.13
## [3] scater_1.22.0 scuttle_1.4.0
## [5] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
## [7] Biobase_2.54.0 GenomicRanges_1.46.0
## [9] GenomeInfoDb_1.30.0 IRanges_2.28.0
## [11] S4Vectors_0.32.0 BiocGenerics_0.40.0
## [13] MatrixGenerics_1.6.0 matrixStats_0.61.0
## [15] SeuratObject_4.0.2 Seurat_4.0.5
## [17] scRepertoire_1.4.0 ggplot2_3.3.5
## [19] BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.22
## [3] tidyselect_1.1.1 htmlwidgets_1.5.4
## [5] grid_4.1.1 BiocParallel_1.28.0
## [7] Rtsne_0.15 munsell_0.5.0
## [9] ScaledMatrix_1.2.0 codetools_0.2-18
## [11] ica_1.0-2 future_1.22.1
## [13] miniUI_0.1.1.1 withr_2.4.2
## [15] colorspace_2.0-2 highr_0.9
## [17] knitr_1.36 ggalluvial_0.12.3
## [19] ROCR_1.0-11 tensor_1.5
## [21] listenv_0.8.0 labeling_0.4.2
## [23] GenomeInfoDbData_1.2.7 polyclip_1.10-0
## [25] farver_2.1.0 parallelly_1.28.1
## [27] vctrs_0.3.8 generics_0.1.1
## [29] xfun_0.27 R6_2.5.1
## [31] doParallel_1.0.16 ggbeeswarm_0.6.0
## [33] rsvd_1.0.5 VGAM_1.1-5
## [35] bitops_1.0-7 spatstat.utils_2.2-0
## [37] DelayedArray_0.20.0 assertthat_0.2.1
## [39] promises_1.2.0.1 beeswarm_0.4.0
## [41] gtable_0.3.0 beachmat_2.10.0
## [43] globals_0.14.0 goftest_1.2-3
## [45] rlang_0.4.12 GlobalOptions_0.1.2
## [47] splines_4.1.1 lazyeval_0.2.2
## [49] spatstat.geom_2.3-0 BiocManager_1.30.16
## [51] yaml_2.2.1 reshape2_1.4.4
## [53] abind_1.4-5 httpuv_1.6.3
## [55] tools_4.1.1 bookdown_0.24
## [57] cubature_2.0.4.2 ellipsis_0.3.2
## [59] spatstat.core_2.3-0 jquerylib_0.1.4
## [61] RColorBrewer_1.1-2 ggridges_0.5.3
## [63] Rcpp_1.0.7 plyr_1.8.6
## [65] sparseMatrixStats_1.6.0 zlibbioc_1.40.0
## [67] purrr_0.3.4 RCurl_1.98-1.5
## [69] rpart_4.1-15 deldir_1.0-6
## [71] pbapply_1.5-0 viridis_0.6.2
## [73] cowplot_1.1.1 zoo_1.8-9
## [75] ggrepel_0.9.1 cluster_2.1.2
## [77] magrittr_2.0.1 data.table_1.14.2
## [79] magick_2.7.3 scattermore_0.7
## [81] SparseM_1.81 lmtest_0.9-38
## [83] RANN_2.6.1 truncdist_1.0-2
## [85] fitdistrplus_1.1-6 gsl_2.1-7
## [87] patchwork_1.1.1 mime_0.12
## [89] evaluate_0.14 xtable_1.8-4
## [91] shape_1.4.6 gridExtra_2.3
## [93] compiler_4.1.1 tibble_3.1.5
## [95] KernSmooth_2.23-20 crayon_1.4.1
## [97] htmltools_0.5.2 mgcv_1.8-38
## [99] later_1.3.0 tidyr_1.1.4
## [101] powerTCR_1.14.0 DBI_1.1.1
## [103] MASS_7.3-54 Matrix_1.3-4
## [105] permute_0.9-5 parallel_4.1.1
## [107] evd_2.3-3 igraph_1.2.7
## [109] pkgconfig_2.0.3 plotly_4.10.0
## [111] spatstat.sparse_2.0-0 foreach_1.5.1
## [113] vipor_0.4.5 bslib_0.3.1
## [115] stringdist_0.9.8 XVector_0.34.0
## [117] stringr_1.4.0 digest_0.6.28
## [119] sctransform_0.3.2 RcppAnnoy_0.0.19
## [121] vegan_2.5-7 spatstat.data_2.1-0
## [123] rmarkdown_2.11 leiden_0.3.9
## [125] uwot_0.1.10 DelayedMatrixStats_1.16.0
## [127] evmix_2.12 shiny_1.7.1
## [129] lifecycle_1.0.1 nlme_3.1-153
## [131] jsonlite_1.7.2 BiocNeighbors_1.12.0
## [133] viridisLite_0.4.0 fansi_0.5.0
## [135] pillar_1.6.4 lattice_0.20-45
## [137] fastmap_1.1.0 httr_1.4.2
## [139] survival_3.2-13 glue_1.4.2
## [141] png_0.1-7 iterators_1.0.13
## [143] stringi_1.7.5 sass_0.4.0
## [145] BiocSingular_1.10.0 dplyr_1.0.7
## [147] irlba_2.3.3 future.apply_1.8.1