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