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

## -----------------------------------------------------------------------------
suppressMessages(library(scRepertoire))

## -----------------------------------------------------------------------------
data("contig_list") #the data built into scRepertoire

head(contig_list[[1]])

## ----eval=FALSE---------------------------------------------------------------
#  for (i in seq_along(contig_list)) {
#      contig_list[[i]] <- stripBarcode(contig_list[[i]], column = 1, connector = "_", num_connects = 3)
#  }

## -----------------------------------------------------------------------------
combined <- combineTCR(contig_list, 
                        samples = c("PY", "PY", "PX", "PX", "PZ","PZ"), 
                        ID = c("P", "T", "P", "T", "P", "T"), cells ="T-AB")

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

## -----------------------------------------------------------------------------
subset <- subsetContig(combined, name = "sample", 
                        variables = c("PX", "PY"))

## -----------------------------------------------------------------------------
quantContig(combined, cloneCall="gene+nt", scale = TRUE)

## -----------------------------------------------------------------------------
quantContig_output <- quantContig(combined, cloneCall="gene+nt", 
                                    scale = TRUE, exportTable = TRUE)
quantContig_output

## -----------------------------------------------------------------------------
quantContig(combined, cloneCall="gene", group = "ID", scale = TRUE)

## -----------------------------------------------------------------------------
abundanceContig(combined, cloneCall = "gene", scale = FALSE)
abundanceContig(combined, cloneCall = "gene", group = "ID", scale = FALSE)

## -----------------------------------------------------------------------------
abundanceContig(combined, group = "ID", scale = TRUE)

## -----------------------------------------------------------------------------
lengthContig(combined, cloneCall="aa", chains = "combined") 

## -----------------------------------------------------------------------------
lengthContig(combined, cloneCall="nt", chains = "single") 

## -----------------------------------------------------------------------------
compareClonotypes(combined, numbers = 10, samples = c("PX_P", "PX_T"), 
                    cloneCall="aa", graph = "alluvial")

## -----------------------------------------------------------------------------
clonalHomeostasis(combined, cloneCall = "gene")
clonalHomeostasis(combined, cloneCall = "aa")

## -----------------------------------------------------------------------------
clonalProportion(combined, cloneCall = "gene") 
clonalProportion(combined, cloneCall = "nt") 

## -----------------------------------------------------------------------------
clonalOverlap(combined, cloneCall = "gene+nt", method = "morisita")

## -----------------------------------------------------------------------------
clonesizeDistribution(combined, cloneCall = "gene+nt", 
                        method="ward.D2")

## -----------------------------------------------------------------------------
clonalDiversity(combined, cloneCall = "gene", group = "samples", n.boots = 100)
clonalDiversity(combined, cloneCall = "gene", group = "ID")

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

## -----------------------------------------------------------------------------
sub_combined <- clusterTCR(combined[[1]], chain = "TCRA", sequence = "aa", threshold = 0.85)

sub_combined[[1]]$TCRA_cluster[1:20]

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

## -----------------------------------------------------------------------------
table(screp_example$seurat_clusters)

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

## -----------------------------------------------------------------------------
colorblind_vector <- colorRampPalette(c("#FF4B20", "#FFB433", "#C6FDEC", "#7AC5FF", "#0348A6"))

DimPlot(screp_example, group.by = "Type") + NoLegend() +
    scale_color_manual(values=colorblind_vector(2))

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

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

plotUMAP(sce, colour_by = "cloneType")

## -----------------------------------------------------------------------------
clonalOverlay(screp_example, reduction = "umap", 
              freq.cutpoint = 30, bins = 10, facet = "Patient") + 
                guides(color = FALSE)

## -----------------------------------------------------------------------------
screp_example <- highlightClonotypes(screp_example, cloneCall= "aa", 
                sequence = c("CAVNGGSQGNLIF_CSAEREDTDTQYF", "NA_CATSATLRVVAEKLFF"))
Seurat::DimPlot(screp_example, group.by = "highlight")

## -----------------------------------------------------------------------------
occupiedscRepertoire(screp_example, x.axis = "cluster")

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

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

## -----------------------------------------------------------------------------
StartracDiversity(screp_example, type = "Type", sample = "Patient", by = "overall")

## -----------------------------------------------------------------------------
combined2 <- expression2List(screp_example, group = "cluster")
combined3 <- expression2List(sce, group = "cluster")

## -----------------------------------------------------------------------------
clonalDiversity(combined2, cloneCall = "nt")
clonalDiversity(combined3, cloneCall = "nt")

## -----------------------------------------------------------------------------
clonalHomeostasis(combined2, cloneCall = "nt")
clonalHomeostasis(combined3, cloneCall = "nt")

## -----------------------------------------------------------------------------
clonalProportion(combined2, cloneCall = "nt")
clonalProportion(combined3, cloneCall = "nt")

## -----------------------------------------------------------------------------
clonalOverlap(combined2, cloneCall="aa", method="overlap")
clonalOverlap(combined3, cloneCall="aa", method="overlap")

## -----------------------------------------------------------------------------
sessionInfo()