## ----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()