## ---- results = "hide"---------------------------------------------------
library(LymphoSeq)

file.path <- system.file("extdata", "TCRB_sequencing", package = "LymphoSeq")

file.list <- readImmunoSeq(path = file.path, 
                           columns = c("aminoAcid", "nucleotide", "count", "count (templates)", 
                                       "count (reads)", "frequencyCount", "frequencyCount (%)",
                                       "estimatedNumberGenomes", "vFamilyName", "dFamilyName", 
                                       "jFamilyName"), 
                           recursive = FALSE)

## ---- comment = ""-------------------------------------------------------
names(file.list)

## ---- comment = ""-------------------------------------------------------
sapply(file.list, dim)

## ---- comment = ""-------------------------------------------------------
CMV <- file.list[grep("CMV", names(file.list))]
names(CMV)
TCRB_Day0_Unsorted <- file.list[["TCRB_Day0_Unsorted"]]
head(TCRB_Day0_Unsorted)

## ---- comment = ""-------------------------------------------------------
metadata <- read.csv(system.file("extdata", "metadata.csv", package = "LymphoSeq"))
metadata
selected <- as.character(metadata[metadata$phenotype == "Unsorted" & 
                                 metadata$day > 300, "samples"])
file.list.selected <- file.list[selected]
names(file.list.selected)

## ---- results = "hide"---------------------------------------------------
productive.aa <- productiveSeq(file.list = file.list, aggregate = "aminoAcid", 
                               prevalence = FALSE)

## ---- results = "hide"---------------------------------------------------
productive.nt <- productiveSeq(file.list = file.list, aggregate = "nucleotide", 
                               prevalence = FALSE)

## ---- comment = ""-------------------------------------------------------
head(file.list[["TCRB_Day0_Unsorted"]])

## ---- comment = ""-------------------------------------------------------
head(productive.aa[["TCRB_Day0_Unsorted"]])

## ---- comment = ""-------------------------------------------------------
head(productive.nt[["TCRB_Day0_Unsorted"]])

## ---- comment = ""-------------------------------------------------------
clonality(file.list = file.list)

## ---- comment = ""-------------------------------------------------------
searchSeq(list = productive.aa, sequence = "CASSPVSNEQFF", type = "aminoAcid", 
          match = "global", editDistance = 0)

## ---- comment = ""-------------------------------------------------------
published <- searchPublished(list = productive.aa)
head(published)

## ---- fig.width = 7, fig.height = 7, comment = ""------------------------
lorenzCurve(samples = names(productive.aa), list = productive.aa)

## ---- fig.width = 7, fig.height = 5, comment = ""------------------------
topSeqsPlot(list = productive.aa, top = 10)

## ---- comment = ""-------------------------------------------------------
bhattacharyya.matrix <- bhattacharyyaMatrix(productive.seqs = productive.aa)
bhattacharyya.matrix[,1:2]
similarity.matrix <- similarityMatrix(productive.seqs = productive.aa)
similarity.matrix[,1:2]

## ---- fig.width = 6.5, fig.height = 5, comment = ""----------------------
pairwisePlot(matrix = bhattacharyya.matrix)

## ---- comment = ""-------------------------------------------------------
common <- commonSeqs(samples = c("TCRB_Day0_Unsorted", "TCRB_Day32_Unsorted"), 
                    productive.aa = productive.aa)
head(common)

## ---- fig.width = 4, fig.height = 4, comment = ""------------------------
commonSeqsVenn(samples = c("TCRB_Day32_Unsorted", "TCRB_Day83_Unsorted"), 
               productive.seqs = productive.aa)

## ---- fig.width = 4, fig.height = 4, comment = ""------------------------
commonSeqsVenn(samples = c("TCRB_Day0_Unsorted", "TCRB_Day32_Unsorted", "TCRB_Day83_Unsorted"), 
               productive.seqs = productive.aa)

## ---- fig.width = 4, fig.height = 4, comment = ""------------------------
commonSeqsPlot("TCRB_Day32_Unsorted", "TCRB_Day83_Unsorted", 
               productive.aa = productive.aa)

## ---- comment = ""-------------------------------------------------------
unique.seqs <- uniqueSeqs(productive.aa = productive.aa)
head(unique.seqs)
sequence.matrix <- seqMatrix(productive.aa = productive.aa, sequences = unique.seqs$aminoAcid)
head(sequence.matrix)[1:6]

## ---- comment = ""-------------------------------------------------------
top.freq <- topFreq(productive.aa = productive.aa, percent = 0.1)
head(top.freq)

## ---- comment = ""-------------------------------------------------------
top.freq <- topFreq(productive.aa = productive.aa, percent = 0)
top.freq.matrix <- merge(top.freq, sequence.matrix)
head(top.freq.matrix)[1:12]

## ---- fig.width = 7, fig.height = 5--------------------------------------
cloneTrack(sequence.matrix = sequence.matrix, 
           productive.aa = productive.aa, 
           map = c("TCRB_Day949_CD4", "TCRB_Day949_CD8"), 
           label = c("CD4", "CD8"), 
           track = "CASSPPTGERDTQYF", 
           unassigned = FALSE)

## ---- comment = ""-------------------------------------------------------
vGenes <- geneFreq(productive.nt = productive.nt, locus = "V", family = TRUE)
head(vGenes)

## ---- fig.width = 4, fig.height = 4--------------------------------------
top.seqs <- topSeqs(productive.seqs = productive.nt, top = 1)
chordDiagramVDJ(sample = top.seqs, 
                association = "VJ", 
                colors = c("darkred", "navyblue"))

## ---- fig.width = 4, fig.height = 4, warning = FALSE, message = FALSE----
vGenes <- geneFreq(productive.nt = productive.nt, locus = "V", family = TRUE)
library(RColorBrewer)
library(grDevices)
RedBlue <- grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(256)
library(wordcloud)
wordcloud::wordcloud(words = vGenes[vGenes$samples == "TCRB_Day83_Unsorted", "familyName"], 
                     freq = vGenes[vGenes$samples == "TCRB_Day83_Unsorted", "frequencyGene"], 
                     colors = RedBlue)

## ---- fig.width = 5, fig.height = 7, warning = FALSE, message = FALSE----
library(reshape)
vGenes <- reshape::cast(vGenes, familyName ~ samples, value = "frequencyGene", sum)
rownames(vGenes) = as.character(vGenes$familyName)
vGenes$familyName = NULL
library(pheatmap)
pheatmap::pheatmap(vGenes, color = RedBlue, scale = "row")

## ---- fig.width = 7, fig.height = 5.6, warning = FALSE, message = FALSE----
vGenes <- geneFreq(productive.nt = productive.nt, locus = "V", family = TRUE)
library(ggplot2)
multicolors <- grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(9, "Set1")))(28)
ggplot2::ggplot(vGenes, aes(x = samples, y = frequencyGene, fill = familyName)) +
  geom_bar(stat = "identity") +
  theme_minimal() + 
  scale_y_continuous(expand = c(0, 0)) + 
  guides(fill = guide_legend(ncol = 2)) +
  scale_fill_manual(values = multicolors) + 
  labs(y = "Frequency (%)", x = "", fill = "") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

## ---- comment = ""-------------------------------------------------------
searchSeq(list = productive.aa, sequence = "CASSDLIGNGKLFF")
cleansed <- removeSeq(file.list = productive.aa, sequence = "CASSDLIGNGKLFF")
searchSeq(list = cleansed, sequence = "CASSDLIGNGKLFF")

## ------------------------------------------------------------------------
TCRB_Day949_Merged <- mergeFiles(samples = c("TCRB_Day949_CD4", "TCRB_Day949_CD8"), 
                                file.list = file.list)

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

