## ----echo=FALSE, results="hide", message=FALSE--------------------------------
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
library(BiocStyle)
library(scater)
## -----------------------------------------------------------------------------
suppressMessages(library(scRepertoire))
## ----eval=FALSE---------------------------------------------------------------
# S1 <- read.csv(".../Sample1/outs/filtered_contig_annotations.csv")
# S2 <- read.csv(".../Sample2/outs/filtered_contig_annotations.csv")
# S3 <- read.csv(".../Sample3/outs/filtered_contig_annotations.csv")
# S4 <- read.csv(".../Sample4/outs/filtered_contig_annotations.csv")
#
# contig_list <- list(S1, S2, S3, S4)
## ----eval=FALSE, tidy = FALSE-------------------------------------------------
# #Directory example
# contig.output <- c("~/Documents/MyExperiment")
# contig.list <- loadContigs(input = contig.output,
# format = "TRUST4")
#
# #List of data frames example
# S1 <- read.csv("~/Documents/MyExperiment/Sample1/outs/barcode_results.csv")
# S2 <- read.csv("~/Documents/MyExperiment/Sample2/outs/barcode_results.csv")
# S3 <- read.csv("~/Documents/MyExperiment/Sample3/outs/barcode_results.csv")
# S4 <- read.csv("~/Documents/MyExperiment/Sample4/outs/barcode_results.csv")
#
# contig_list <- list(S1, S2, S3, S4)
# contig.list <- loadContigs(input = contig.output,
# format = "WAT3R")
## ----eval = F, tidy = FALSE---------------------------------------------------
# contigs <- read.csv(".../outs/filtered_contig_annotations.csv")
#
# contig.list <- createHTOContigList(contigs,
# Seurat.Obj,
# group.by = "HTO_maxID")
## ----tidy = FALSE-------------------------------------------------------------
data("contig_list") #the data built into scRepertoire
head(contig_list[[1]])
## ----tidy = FALSE-------------------------------------------------------------
combined.TCR <- combineTCR(contig_list,
samples = c("P17B", "P17L", "P18B", "P18L",
"P19B","P19L", "P20B", "P20L"),
removeNA = FALSE,
removeMulti = FALSE,
filterMulti = FALSE)
head(combined.TCR[[1]])
## ----tidy = FALSE-------------------------------------------------------------
BCR.contigs <- read.csv("https://www.borch.dev/uploads/contigs/b_contigs.csv")
combined.BCR <- combineBCR(BCR.contigs,
samples = "P1",
threshold = 0.85)
head(combined.BCR[[1]])
## ----tidy = FALSE-------------------------------------------------------------
combined.TCR <- addVariable(combined.TCR,
variable.name = "Type",
variables = rep(c("B", "L"), 4))
head(combined.TCR[[1]])
## ----tidy = FALSE-------------------------------------------------------------
subset1 <- subsetClones(combined.TCR,
name = "sample",
variables = c("P18L", "P18B"))
head(subset1[[1]])
## -----------------------------------------------------------------------------
subset2 <- combined.TCR[c(3,4)]
head(subset2[[1]])
## ----eval = FALSE, tidy = FALSE-----------------------------------------------
# exportClones(combined,
# write.file = TRUE,
# dir = "~/Documents/MyExperiment/Sample1/"
# file.name = "clones.csv"
## ----tidy = FALSE-------------------------------------------------------------
clonalQuant(combined.TCR,
cloneCall="strict",
chain = "both",
scale = TRUE)
## -----------------------------------------------------------------------------
clonalQuant(combined.TCR,
cloneCall="gene",
group.by = "Type",
scale = TRUE)
## ----tidy = FALSE-------------------------------------------------------------
clonalAbundance(combined.TCR,
cloneCall = "gene",
scale = FALSE)
## -----------------------------------------------------------------------------
clonalAbundance(combined.TCR,
cloneCall = "gene",
scale = TRUE)
## ----tidy = FALSE-------------------------------------------------------------
clonalLength(combined.TCR,
cloneCall="aa",
chain = "both")
clonalLength(combined.TCR,
cloneCall="aa",
chain = "TRA",
scale = TRUE)
## ----tidy = FALSE-------------------------------------------------------------
clonalCompare(combined.TCR,
top.clones = 10,
samples = c("P17B", "P17L"),
cloneCall="aa",
graph = "alluvial")
## ----tidy = FALSE-------------------------------------------------------------
clonalCompare(combined.TCR,
top.clones = 10,
highlight.clones = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF"),
relabel.clones = TRUE,
samples = c("P17B", "P17L"),
cloneCall="aa",
graph = "alluvial")
## -----------------------------------------------------------------------------
clonalCompare(combined.TCR,
clones = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF"),
relabel.clones = TRUE,
samples = c("P17B", "P17L"),
cloneCall="aa",
graph = "alluvial")
## ----tidy = FALSE-------------------------------------------------------------
clonalScatter(combined.TCR,
cloneCall ="gene",
x.axis = "P18B",
y.axis = "P18L",
dot.size = "total",
graph = "proportion")
## ----tidy = FALSE-------------------------------------------------------------
clonalHomeostasis(combined.TCR,
cloneCall = "gene")
## ----tidy = FALSE-------------------------------------------------------------
clonalHomeostasis(combined.TCR,
cloneCall = "gene",
cloneSize = c(Rare = 0.001, Small = 0.01, Medium = 0.1, Large = 0.3, Hyperexpanded =
1))
## ----tidy = FALSE-------------------------------------------------------------
combined.TCR <- addVariable(combined.TCR,
variable.name = "Type",
variables = rep(c("B", "L"), 4))
clonalHomeostasis(combined.TCR,
group.by = "Type",
cloneCall = "gene")
## ----tidy = FALSE-------------------------------------------------------------
clonalProportion(combined.TCR,
cloneCall = "gene")
clonalProportion(combined.TCR,
cloneCall = "nt",
clonalSplit = c(1, 5, 10, 100, 1000, 10000))
## ----tidy = FALSE-------------------------------------------------------------
vizGenes(combined.TCR,
x.axis = "TRBV",
y.axis = NULL,
plot = "barplot",
scale = TRUE)
## ----tidy = FALSE-------------------------------------------------------------
#Peripheral Blood
vizGenes(combined.TCR[c(1,3,5,7)],
x.axis = "TRBV",
y.axis = "TRBJ",
plot = "heatmap",
scale = TRUE)
#Lung
vizGenes(combined.TCR[c(2,4,6,8)],
x.axis = "TRBV",
y.axis = "TRBJ",
plot = "heatmap",
scale = TRUE)
## ----tidy = FALSE-------------------------------------------------------------
vizGenes(combined.TCR[c(1,2)],
x.axis = "TRBV",
y.axis = "TRAV",
plot = "heatmap",
scale = FALSE)
## ----message = FALSE, tidy = FALSE--------------------------------------------
percentAA(combined.TCR,
chain = "TRB",
aa.length = 20)
## -----------------------------------------------------------------------------
positionalEntropy(combined.TCR,
chain = "TRB",
aa.length = 20)
## ----tidy = FALSE-------------------------------------------------------------
percentGenes(combined.TCR,
chain = "TRB",
gene = "Vgene")
## ----tidy = FALSE-------------------------------------------------------------
df.genes <- percentGenes(combined.TCR,
chain = "TRB",
gene = "Vgene",
exportTable = TRUE)
#Performing PCA
pc <- prcomp(df.genes)
#Getting data frame to plot from
df <- as.data.frame(cbind(pc$x[,1:2], rownames(df.genes)))
df$PC1 <- as.numeric(df$PC1)
df$PC2 <- as.numeric(df$PC2)
#Plotting
ggplot(df, aes(x = PC1, y = PC2)) +
geom_point(aes(fill =df[,3]), shape = 21, size = 5) +
guides(fill=guide_legend(title="Samples")) +
scale_fill_manual(values = hcl.colors(nrow(df), "inferno")) +
theme_classic()
## ----tidy = FALSE-------------------------------------------------------------
percentVJ(combined.TCR,
chain = "TRB")
df.genes <- percentVJ(combined.TCR,
chain = "TRB",
exportTable = TRUE)
#Performing PCA
pc <- prcomp(df.genes)
#Getting data frame to plot from
df <- as.data.frame(cbind(pc$x[,1:2], rownames(df.genes)))
df$PC1 <- as.numeric(df$PC1)
df$PC2 <- as.numeric(df$PC2)
#Plotting
ggplot(df, aes(x = PC1, y = PC2)) +
geom_point(aes(fill =df[,3]), shape = 21, size = 5) +
guides(fill=guide_legend(title="Samples")) +
scale_fill_manual(values = hcl.colors(nrow(df), "inferno")) +
theme_classic()
## ----tidy = FALSE-------------------------------------------------------------
percentKmer(combined.TCR,
cloneCall = "aa",
chain = "TRB",
motif.length = 3,
top.motifs = 25)
percentKmer(combined.TCR,
cloneCall = "nt",
chain = "TRB",
motif.length = 3,
top.motifs = 25)
## ----tidy = FALSE-------------------------------------------------------------
clonalDiversity(combined.TCR,
cloneCall = "gene",
n.boots = 20)
## -----------------------------------------------------------------------------
#Return only a subset of metrics
clonalDiversity(combined.TCR,
metrics = c("shannon", "ACE"),
cloneCall = "gene",
n.boots = 20)
## ----message=FALSE, tidy = FALSE----------------------------------------------
clonalRarefaction(combined.TCR,
plot.type = 1,
hill.numbers = 0,
n.boots = 2)
clonalRarefaction(combined.TCR,
plot.type = 2,
hill.numbers = 0,
n.boots = 2)
clonalRarefaction(combined.TCR,
plot.type = 3,
hill.numbers = 0,
n.boots = 2)
## ----tidy = FALSE-------------------------------------------------------------
clonalRarefaction(combined.TCR,
plot.type = 1,
hill.numbers = 1,
n.boots = 2)
clonalRarefaction(combined.TCR,
plot.type = 2,
hill.numbers = 1,
n.boots = 2)
clonalRarefaction(combined.TCR,
plot.type = 3,
hill.numbers = 1,
n.boots = 2)
## ----tidy = FALSE-------------------------------------------------------------
clonalSizeDistribution(combined.TCR,
cloneCall = "aa",
method= "ward.D2")
## ----tidy = FALSE-------------------------------------------------------------
clonalOverlap(combined.TCR,
cloneCall = "strict",
method = "morisita")
clonalOverlap(combined.TCR,
cloneCall = "strict",
method = "raw")
## ----tidy = FALSE-------------------------------------------------------------
scRep_example <- get(data("scRep_example"))
#Making a Single-Cell Experiment object
sce <- Seurat::as.SingleCellExperiment(scRep_example)
## ----tidy = FALSE-------------------------------------------------------------
sce <- combineExpression(combined.TCR,
sce,
cloneCall="gene",
group.by = "sample",
proportion = TRUE)
#Define color palette
colorblind_vector <- hcl.colors(n=7, palette = "inferno", fixup = TRUE)
plotUMAP(sce, colour_by = "cloneSize") +
scale_color_manual(values=rev(colorblind_vector[c(1,3,5,7)]))
## ----tidy = FALSE-------------------------------------------------------------
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))
Seurat::DimPlot(scRep_example, group.by = "cloneSize") +
scale_color_manual(values=rev(colorblind_vector[c(1,3,4,5,7)]))
## ----eval=FALSE, tidy = FALSE-------------------------------------------------
# #This is an example of the process, which will not be evaluated during knit
# TCR <- combineTCR(...)
# BCR <- combineBCR(...)
# list.receptors <- c(TCR, BCR)
#
#
# seurat <- combineExpression(list.receptors,
# seurat,
# cloneCall="gene",
# proportion = TRUE)
## ----tidy = FALSE-------------------------------------------------------------
#Adding patient information
scRep_example$Patient <- substr(scRep_example$orig.ident, 1,3)
clonalOverlay(scRep_example,
reduction = "umap",
cutpoint = 1,
bins = 10,
facet.by = "Patient") +
guides(color = "none")
## ----tidy = FALSE-------------------------------------------------------------
#ggraph needs to be loaded due to issues with ggplot
library(ggraph)
#No Identity filter
clonalNetwork(scRep_example,
reduction = "umap",
group.by = "seurat_clusters",
filter.clones = NULL,
filter.identity = NULL,
cloneCall = "aa")
## ----tidy = FALSE-------------------------------------------------------------
#Examining Cluster 3 only
clonalNetwork(scRep_example,
reduction = "umap",
group.by = "seurat_clusters",
filter.identity = 3,
cloneCall = "aa")
## ----tidy = FALSE-------------------------------------------------------------
shared.clones <- clonalNetwork(scRep_example,
reduction = "umap",
group.by = "seurat_clusters",
cloneCall = "aa",
exportClones = TRUE)
head(shared.clones)
## ----tidy = FALSE-------------------------------------------------------------
scRep_example <- highlightClones(scRep_example,
cloneCall= "aa",
sequence = c("CAERGSGGSYIPTF_CASSDPSGRQGPRWDTQYF",
"CARKVRDSSYKLIF_CASSDSGYNEQFF"))
Seurat::DimPlot(scRep_example, group.by = "highlight") +
ggplot2::theme(plot.title = element_blank())
## ----tidy = FALSE-------------------------------------------------------------
clonalOccupy(scRep_example,
x.axis = "seurat_clusters")
clonalOccupy(scRep_example,
x.axis = "ident",
proportion = TRUE,
label = FALSE)
## ----tidy = FALSE-------------------------------------------------------------
#Adding type information
scRep_example$Type <- substr(scRep_example$orig.ident, 4,4)
alluvialClones(scRep_example,
cloneCall = "aa",
y.axes = c("Patient", "ident", "Type"),
color = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF")) +
scale_fill_manual(values = c("grey", colorblind_vector[3]))
alluvialClones(scRep_example,
cloneCall = "gene",
y.axes = c("Patient", "ident", "Type"),
color = "ident")
## ----tidy = FALSE-------------------------------------------------------------
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)
## -----------------------------------------------------------------------------
subset <- subset(scRep_example, Type == "L")
circles <- getCirclize(subset, group.by = "ident", proportion = TRUE)
grid.cols <- scales::hue_pal()(length(unique(subset@active.ident)))
names(grid.cols) <- levels(subset@active.ident)
chordDiagram(circles,
self.link = 1,
grid.col = grid.cols,
directional = 1,
direction.type = "arrows",
link.arr.type = "big.arrow")
## ----tidy = FALSE-------------------------------------------------------------
StartracDiversity(scRep_example,
type = "Type",
group.by = "Patient")
## ----message = FALSE, tidy = FALSE--------------------------------------------
clonalBias(scRep_example,
cloneCall = "aa",
split.by = "Patient",
group.by = "seurat_clusters",
n.boots = 10,
min.expand =0)
## ----tidy = FALSE-------------------------------------------------------------
sub_combined <- clonalCluster(combined.TCR[[2]],
chain = "TRA",
sequence = "aa",
threshold = 0.85,
group.by = NULL)
head(sub_combined[,c(1,2,13)])
## ----tidy = FALSE-------------------------------------------------------------
#Define color palette
colorblind_vector <- hcl.colors(n=7, palette = "inferno", fixup = TRUE)
scRep_example <- clonalCluster(scRep_example,
chain = "TRA",
sequence = "aa",
threshold = 0.85,
group.by = "Patient")
Seurat::DimPlot(scRep_example, group.by = "TRA_cluster") +
scale_color_manual(values = hcl.colors(n=length(unique(scRep_example@meta.data[,"TRA_cluster"])), "inferno")) +
Seurat::NoLegend() +
theme(plot.title = element_blank())
## ----tidy = FALSE-------------------------------------------------------------
#Clustering Patient 19 samples
igraph.object <- clonalCluster(combined.TCR[c(5,6)],
chain = "TRB",
sequence = "aa",
group.by = "sample",
threshold = 0.85,
exportGraph = TRUE)
#Setting color scheme
col_legend <- factor(igraph::V(igraph.object)$group)
col_samples <- hcl.colors(3,"inferno")[as.numeric(col_legend)]
color.legend <- factor(unique(igraph::V(igraph.object)$group))
#Plotting
plot(
igraph.object,
vertex.size = sqrt(igraph::V(igraph.object)$size),
vertex.label = NA,
edge.arrow.size = .25,
vertex.color = col_samples
)
legend("topleft", legend = levels(color.legend), pch = 16, col = unique(col_samples), bty = "n")
## -----------------------------------------------------------------------------
sessionInfo()