## ----results = "hide"--------------------------------------------------------- library(LymphoSeq) TCRB.path <- system.file("extdata", "TCRB_sequencing", package = "LymphoSeq") TCRB.list <- readImmunoSeq(path = TCRB.path) ## ----comment = ""------------------------------------------------------------- names(TCRB.list) ## ----comment = ""------------------------------------------------------------- lapply(TCRB.list, dim) ## ----comment = ""------------------------------------------------------------- CMV <- TCRB.list[grep("CMV", names(TCRB.list))] names(CMV) TRB_Unsorted_0 <- TCRB.list[["TRB_Unsorted_0"]] head(TRB_Unsorted_0) ## ----comment = ""------------------------------------------------------------- TCRB.metadata <- read.csv(system.file("extdata", "TCRB_metadata.csv", package = "LymphoSeq")) TCRB.metadata selected <- as.character(TCRB.metadata[TCRB.metadata$phenotype == "Unsorted" & TCRB.metadata$day > 300, "samples"]) TCRB.list.selected <- TCRB.list[selected] names(TCRB.list.selected) ## ----results = "hide"--------------------------------------------------------- productive.TRB.aa <- productiveSeq(file.list = TCRB.list, aggregate = "aminoAcid", prevalence = FALSE) ## ----results = "hide"--------------------------------------------------------- productive.TRB.nt <- productiveSeq(file.list = TCRB.list, aggregate = "nucleotide", prevalence = FALSE) ## ----comment = ""------------------------------------------------------------- head(TCRB.list[["TRB_Unsorted_0"]]) ## ----comment = ""------------------------------------------------------------- head(productive.TRB.aa[["TRB_Unsorted_0"]]) ## ----comment = ""------------------------------------------------------------- head(productive.TRB.nt[["TRB_Unsorted_0"]]) ## ----comment = ""------------------------------------------------------------- clonality(file.list = TCRB.list) ## ----results = "hide"--------------------------------------------------------- IGH.path <- system.file("extdata", "IGH_sequencing", package = "LymphoSeq") IGH.list <- readImmunoSeq(path = IGH.path) ## ----comment = ""------------------------------------------------------------- clonalRelatedness(list = IGH.list, editDistance = 10) ## ----results = "hide"--------------------------------------------------------- productive.IGH.nt <- productiveSeq(file.list = IGH.list, aggregate = "nucleotide") ## ----fig.width = 7, fig.height = 8, comment = ""------------------------------ phyloTree(list = productive.IGH.nt, sample = "IGH_MVQ92552A_BL", type = "nucleotide", layout = "rectangular") ## ----comment = ""------------------------------------------------------------- alignSeq(list = productive.IGH.nt, sample = "IGH_MVQ92552A_BL", type = "aminoAcid", method = "ClustalW", output = "consule") ## ----comment = ""------------------------------------------------------------- searchSeq(list = productive.TRB.aa, sequence = "CASSPVSNEQFF", type = "aminoAcid", match = "global", editDistance = 0) ## ----comment = ""------------------------------------------------------------- published <- searchPublished(list = productive.TRB.aa) head(published) ## ----fig.width = 7, fig.height = 7, comment = ""------------------------------ lorenzCurve(samples = names(productive.TRB.aa), list = productive.TRB.aa) ## ----fig.width = 7, fig.height = 5, comment = ""------------------------------ topSeqsPlot(list = productive.TRB.aa, top = 10) ## ----comment = ""------------------------------------------------------------- bhattacharyya.matrix <- bhattacharyyaMatrix(productive.seqs = productive.TRB.aa) bhattacharyya.matrix similarity.matrix <- similarityMatrix(productive.seqs = productive.TRB.aa) similarity.matrix ## ----fig.width = 6.5, fig.height = 5, comment = ""---------------------------- pairwisePlot(matrix = bhattacharyya.matrix) ## ----comment = ""------------------------------------------------------------- common <- commonSeqs(samples = c("TRB_Unsorted_0", "TRB_Unsorted_32"), productive.aa = productive.TRB.aa) head(common) ## ----fig.width = 4, fig.height = 4, comment = ""------------------------------ commonSeqsVenn(samples = c("TRB_Unsorted_32", "TRB_Unsorted_83"), productive.seqs = productive.TRB.aa) ## ----fig.width = 4, fig.height = 4, comment = ""------------------------------ commonSeqsVenn(samples = c("TRB_Unsorted_0", "TRB_Unsorted_32", "TRB_Unsorted_83"), productive.seqs = productive.TRB.aa) ## ----fig.width = 4, fig.height = 4, comment = ""------------------------------ commonSeqsPlot("TRB_Unsorted_32", "TRB_Unsorted_83", productive.aa = productive.TRB.aa, show = "common") ## ----fig.width = 7, fig.height = 5, comment = ""------------------------------ commonSeqsBar(productive.aa = productive.TRB.aa, samples = c("TRB_CD4_949", "TRB_CD8_949", "TRB_Unsorted_949", "TRB_Unsorted_1320"), color.sample = "TRB_CD8_949", labels = "no") ## ----comment = "", warning = FALSE-------------------------------------------- differentialAbundance(list = productive.TRB.aa, sample1 = "TRB_Unsorted_949", sample2 = "TRB_Unsorted_1320", type = "aminoAcid", q = 0.01) ## ----comment = ""------------------------------------------------------------- unique.seqs <- uniqueSeqs(productive.aa = productive.TRB.aa) head(unique.seqs) sequence.matrix <- seqMatrix(productive.aa = productive.TRB.aa, sequences = unique.seqs$aminoAcid) head(sequence.matrix) ## ----comment = ""------------------------------------------------------------- top.freq <- topFreq(productive.aa = productive.TRB.aa, percent = 0.1) head(top.freq) ## ----comment = ""------------------------------------------------------------- top.freq <- topFreq(productive.aa = productive.TRB.aa, percent = 0) top.freq.matrix <- merge(top.freq, sequence.matrix) head(top.freq.matrix) ## ----fig.width = 7, fig.height = 5-------------------------------------------- cloneTrack(sequence.matrix = sequence.matrix, productive.aa = productive.TRB.aa, map = c("TRB_CD4_949", "TRB_CD8_949"), label = c("CD4", "CD8"), track = "CASSPPTGERDTQYF", unassigned = FALSE) ## ----comment = ""------------------------------------------------------------- vGenes <- geneFreq(productive.nt = productive.TRB.nt, locus = "V", family = TRUE) head(vGenes) ## ----fig.width = 4, fig.height = 4-------------------------------------------- top.seqs <- topSeqs(productive.seqs = productive.TRB.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.TRB.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 == "TRB_Unsorted_83", "familyName"], freq = vGenes[vGenes$samples == "TRB_Unsorted_83", "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.TRB.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.TRB.aa, sequence = "CASSESAGSTGELFF") cleansed <- removeSeq(file.list = productive.TRB.aa, sequence = "CASSESAGSTGELFF") searchSeq(list = cleansed, sequence = "CASSESAGSTGELFF") ## ----------------------------------------------------------------------------- TRB_949_Merged <- mergeFiles(samples = c("TRB_CD4_949", "TRB_CD8_949"), file.list = TCRB.list) ## ----------------------------------------------------------------------------- sessionInfo()