## ----options, results="hide", include=FALSE, cache=FALSE, results='hide', message=FALSE---- knitr::opts_chunk$set(fig.align="center", cache=FALSE,error=FALSE, #make it stop on error fig.width=7, fig.height=7, autodep=TRUE, out.width="600px", out.height="600px", results="markup", echo=TRUE, eval=TRUE) #knitr::opts_knit$set(stop_on_error = 2L) #really make it stop #knitr::dep_auto() options(getClass.msg=FALSE) #get rid of annoying messages about cache until fixed internally in R set.seed(98883) ## for reproducibility library(bioc2016singlecell) library(slingshot) library(clusterExperiment) library(gam) ## ----datain-------------------------------------------------------------- data('full_pca') ## Examine dimensionality reduction # plot3d(pcaX, aspect = 'iso') pairs(pcaX[,1:3], asp = 1) ## ----lines_unsup--------------------------------------------------------- l1 <- get_lineages(pcaX, clus) # plot_tree(pcaX, clus, l1, threeD = TRUE) plot_tree(pcaX, clus, l1, dim = 3) ## ----lines_sup_start----------------------------------------------------- l2 <- get_lineages(pcaX, clus, start.clus = 'm10') # plot_tree(pcaX, clus, l2, threeD = TRUE) plot_tree(pcaX, clus, l2, dim = 3) ## ----lines_sup_end------------------------------------------------------- l3 <- get_lineages(pcaX, clus, start.clus = 'm10', end.clus = 'm17') # plot_tree(pcaX, clus, l3, threeD = TRUE) plot_tree(pcaX, clus, l3, dim = 3) ## ----curves-------------------------------------------------------------- crv <- get_curves(pcaX, clus, l2) # plot_curves(pcaX, clus, c, threeD = TRUE) plot_curves(pcaX, clus, crv, dim = 3) ## ----genedata------------------------------------------------------------ data('var_genes') ## ----fitgam-------------------------------------------------------------- gam.pval <- vector("list",length(crv)) for(l in 1:length(crv)){ t <- crv[[l]]$pseudotime y <- vargenes[,! is.na(t)] t <- t[! is.na(t)] gam.pval[[l]] <- apply(y,1,function(z){ d <- data.frame(z=z, t=t) tmp <- gam(z ~ lo(t), data=d) p <- summary(tmp)[4][[1]][1,5] p }) } ## ----heatmaps------------------------------------------------------------ topgenes1 <- names(sort(gam.pval[[1]], decreasing = FALSE))[1:100] heatdata1 <- vargenes[rownames(vargenes) %in% topgenes1, order(crv[[1]]$pseudotime, na.last = NA)] heatclus1 <- clus[order(crv[[1]]$pseudotime, na.last = NA)] ce1 <- clusterExperiment(heatdata1, heatclus1, transformation=identity) plotHeatmap(ce1, clusterSamplesData="orderSamplesValue") topgenes2 <- names(sort(gam.pval[[2]], decreasing = FALSE))[1:100] heatdata2 <- vargenes[rownames(vargenes) %in% topgenes2, order(crv[[2]]$pseudotime, na.last = NA)] heatclus2 <- clus[order(crv[[2]]$pseudotime, na.last = NA)] ce2 <- clusterExperiment(heatdata2, heatclus2, transformation=identity) plotHeatmap(ce2, clusterSamplesData="orderSamplesValue") ## ----session------------------------------------------------------------- sessionInfo()