## ----echo=FALSE, results="hide", message=FALSE--------------------------------
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
library(BiocStyle)
library(scater)
library(Seurat)

## -----------------------------------------------------------------------------
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"

## ----eval = FALSE, tidy = FALSE-----------------------------------------------
# screp_example <- quietVDJgenes(screp_example)

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

