## ---- echo = FALSE------------------------------------------------------------ library(knitr) # For a smaller vignette # This can be deleted when trying to follow along opts_chunk$set(dpi = 58) ## ---- eval = FALSE------------------------------------------------------------ # if(!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # BiocManager::install("tradeSeq") ## ---- warning=F, message=F---------------------------------------------------- library(tradeSeq) library(RColorBrewer) library(SingleCellExperiment) library(slingshot) # For reproducibility RNGversion("3.5.0") palette(brewer.pal(8, "Dark2")) data(countMatrix, package = "tradeSeq") counts <- countMatrix rm(countMatrix) data(crv, package = "tradeSeq") data(celltype, package = "tradeSeq") ## ---- warning=FALSE, out.width="50%", fig.asp=.6------------------------------ set.seed(200) rd <- reducedDims(crv) cl <- kmeans(rd, centers = 7)$cluster plot(rd, col = brewer.pal(9, "Set1")[cl], pch = 16, asp = 1, cex = 2/3) lin <- getLineages(rd, clusterLabels = cl, start.clus = 1) crv <- getCurves(lin) ## ---- out.width="50%", fig.asp=.6--------------------------------------------- plotGeneCount(curve = crv, clusters = cl) plotGeneCount(curve = crv, clusters = celltype, title = "Colored by cell type") ## ---- eval=FALSE-------------------------------------------------------------- # icMat <- evaluateK(counts = counts, sds = crv, k=3:20, nGenes = 200, # verbose=FALSE) # # alternatively: # icMat <- evaluateK(counts = counts, k=3:20, nGenes = 200, # pseudotime = slingPseudotime(crv, na = FALSE), # cellWeights = slingCurveWeights(crv)) ## ---- echo=FALSE-------------------------------------------------------------- knitr::include_graphics("evalK_paul1Cropped.png") knitr::include_graphics("evalK_paul2Cropped.png") ## ----------------------------------------------------------------------------- library(BiocParallel) library(magrittr) # Register BiocParallel Serial Execution (no parallelization in that case) BiocParallel::register(BiocParallel::SerialParam()) counts <- as.matrix(counts) sce <- fitGAM(counts = counts, sds = crv) # This takes about 1mn to run ## ---- out.width="50%", fig.asp=.6--------------------------------------------- plotGeneCount(curve = crv, counts = counts, clusters = cl, models = sce) ## ----------------------------------------------------------------------------- assoRes <- associationTest(sce) head(assoRes) ## ----------------------------------------------------------------------------- startRes <- startVsEndTest(sce) ## ---- out.width="40%", fig.asp=1---------------------------------------------- oStart <- order(startRes$waldStat, decreasing = TRUE) sigGeneStart <- names(sce)[oStart[1]] plotSmoothers(sce, counts, gene = sigGeneStart) ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(crv, counts, gene = sigGeneStart) ## ----------------------------------------------------------------------------- customRes <- startVsEndTest(sce, pseudotimeValues = c(0.1, 0.8)) ## ----------------------------------------------------------------------------- endRes <- diffEndTest(sce) ## ---- out.width="40%", fig.asp=1---------------------------------------------- o <- order(endRes$waldStat, decreasing = TRUE) sigGene <- names(sce)[o[1]] plotSmoothers(sce, counts, sigGene) ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(crv, counts, gene = sigGene) ## ---- out.width="40%", fig.asp=1---------------------------------------------- patternRes <- patternTest(sce) oPat <- order(patternRes$waldStat, decreasing = TRUE) head(rownames(patternRes)[oPat]) plotSmoothers(sce, counts, gene = rownames(patternRes)[oPat][1]) ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(crv, counts, gene = rownames(patternRes)[oPat][1]) ## ---- out.width="50%", fig.asp=.8--------------------------------------------- library(dplyr) library(ggplot2) library(tidyr) compare <- inner_join(patternRes %>% mutate(Gene = rownames(patternRes), pattern = waldStat) %>% select(Gene, pattern), endRes %>% mutate(Gene = rownames(endRes), end = waldStat) %>% select(Gene, end), by = c("Gene" = "Gene")) %>% mutate(transientScore = (min_rank(desc(end)))^2 + (dense_rank(pattern))^2) ggplot(compare, aes(x = log(pattern), y = log(end))) + geom_point(aes(col = transientScore)) + labs(x = "patternTest Wald Statistic (log scale)", y = "diffEndTest Wald Statistic (log scale)") + scale_color_continuous(low = "yellow", high = "red") + theme_classic() ## ---- out.width="40%", fig.asp=1---------------------------------------------- topTransient <- (compare %>% arrange(desc(transientScore)))[1, "Gene"] plotSmoothers(sce, counts, gene = topTransient) ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(crv, counts, gene = topTransient) ## ----------------------------------------------------------------------------- head(compare %>% arrange(desc(transientScore)) %>% select(Gene), n = 5) ## ---- out.width="40%", fig.asp=1---------------------------------------------- plotSmoothers(sce, counts, gene = "Irf8") ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(crv, counts, gene = "Irf8") ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(curve = crv, counts = counts, clusters = cl, models = sce) earlyDERes <- earlyDETest(sce, knots = c(1, 2)) oEarly <- order(earlyDERes$waldStat, decreasing = TRUE) head(rownames(earlyDERes)[oEarly]) ## ---- out.width="40%", fig.asp=1---------------------------------------------- plotSmoothers(sce, counts, gene = rownames(earlyDERes)[oEarly][2]) ## ---- out.width="50%", fig.asp=.5--------------------------------------------- plotGeneCount(crv, counts, gene = rownames(earlyDERes)[oEarly][2]) ## ---- warning=FALSE,message=F------------------------------------------------- library(clusterExperiment) nPointsClus <- 20 clusPat <- clusterExpressionPatterns(sce, nPoints = nPointsClus, genes = rownames(counts)[1:200]) clusterLabels <- primaryCluster(clusPat$rsec) ## ---- out.width="50%", fig.asp=1---------------------------------------------- library(cowplot) cUniq <- unique(clusterLabels) cUniq <- cUniq[!cUniq == -1] # remove unclustered genes plots <- list() for (xx in cUniq[1:4]) { cId <- which(clusterLabels == xx) p <- ggplot(data = data.frame(x = 1:nPointsClus, y = rep(range(clusPat$yhatScaled[cId, ]), nPointsClus / 2)), aes(x = x, y = y)) + geom_point(alpha = 0) + labs(title = paste0("Cluster ", xx), x = "Pseudotime", y = "Normalized expression") + theme_classic() for (ii in 1:length(cId)) { geneId <- rownames(clusPat$yhatScaled)[cId[ii]] p <- p + geom_line(data = data.frame(x = rep(1:nPointsClus, 2), y = clusPat$yhatScaled[geneId, ], lineage = rep(0:1, each = nPointsClus)), aes(col = as.character(lineage), group = lineage), lwd = 1.5) } p <- p + guides(color = FALSE) + scale_color_manual(values = c("orange", "darkseagreen3"), breaks = c("0", "1")) plots[[as.character(xx)]] <- p } plots$ncol <- 2 do.call(plot_grid, plots) ## ----------------------------------------------------------------------------- gamList <- fitGAM(counts, pseudotime = slingPseudotime(crv, na = FALSE), cellWeights = slingCurveWeights(crv), nknots = 6) ## ----------------------------------------------------------------------------- summary(gamList[["Irf8"]]) ## ---- eval=FALSE-------------------------------------------------------------- # pvalLineage <- getSmootherPvalues(gamList) # statLineage <- getSmootherTestStats(gamList) ## ----------------------------------------------------------------------------- library(mgcv) control <- gam.control() control$maxit <- 1000 #set maximum number of iterations to 1K # pass to control argument of fitGAM as below: # # gamList <- fitGAM(counts = counts, # pseudotime = slingPseudotime(crv, na = FALSE), # cellWeights = slingCurveWeights(crv), # control = control) ## ---- echo = F---------------------------------------------------------------- # ggdraw() + draw_image("cheatsheet_highRes.jpeg") knitr::include_graphics("cheatsheet_highRes.jpeg") ## ----------------------------------------------------------------------------- sessionInfo()