library(devtools) load_all('../..') library(rmarkdown) render('Windram-Monocle.Rmd')
library(knitr) library(knitcitations) library(rmarkdown) # # knitr options # opts_chunk$set( fig.path = 'figures/Windram-', stop_on_error = TRUE, fig.width = 12.5, fig.height = 8) # # Citations # cleanbib() cite_options( # hyperlink = 'to.doc', hyperlink = TRUE, # style = 'html', # citation_format = 'text', citation_format = "pandoc", cite.style = "numeric", check.entries = TRUE) # hyperlink = TRUE) bib.file <- system.file("inst/scripts/DeLorean.bib", package="DeLorean") bib <- read.bibtex(bib.file) if (file.exists("config.R")) { source("config.R") } source(system.file("inst/scripts/shared.R", package="DeLorean"))
# suppressMessages(loadfonts()) library(DeLorean) # # Stylesheet # options(markdown.HTML.stylesheet = system.file("inst/Rmd/foghorn.css", package="DeLorean")) font.family <- "Verdana" font.theme <- theme_update(text=element_text(family=font.family)) theme_set(font.theme)
Load the data that we used for the DeLorean analysis and create a
CellDataSet
for use with Monocle
.
library(monocle) # vignette('monocle-vignette') # Read the input, this is generated by the Windram-DeLorean.Rmd script which # must be run first. input <- readRDS('Data/Windram-input.rds') colnames(input$expr) <- rownames(input$cell.meta) <- input$cell.meta$cell rownames(input$expr) <- rownames(input$gene.meta) <- input$gene.meta$gene cds <- new("CellDataSet", exprs=exp(input$expr), phenoData=new("AnnotatedDataFrame", data=as.data.frame(input$cell.meta)), featureData=new("AnnotatedDataFrame", data=as.data.frame(input$gene.meta))) cds@expressionFamily <- gaussianff()
# Log-transform each value in the expression matrix. L <- log(exprs(cds)) # Standardize each gene, so that they are all on the same scale, # Then melt the data with reshape2 so we can plot it easily melted_dens_df <- reshape2::melt(t(scale(t(L)))) # Plot the distribution of the standardized gene expression values. qplot(value, geom = "density", data =melted_dens_df) + stat_function(fun=dnorm, size=0.5, color='red') + # + xlab("Standardized log(FPKM)") ylab("Density")
set.seed(37) # For reproducibility cds <- reduceDimension(cds, use_irlba=F)
cds <- orderCells(cds, num_paths=1, reverse=F)
Plot the spanning tree.
gp.mst <- plot_spanning_tree(cds, color_by="capture") print(gp.mst) pdf('Windram-Monocle-order.pdf', width=fig.width, height=fig.height) print(gp.mst) dev.off() png('Windram-Monocle-order.png', width=fig.width, height=fig.height, units="in", res=300) print(gp.mst) dev.off()
Plot some genes in pseudotime.
genes.to.plot <- sample(rownames(fData(cds)), 3) p <- plot_genes_in_pseudotime(cds[genes.to.plot,], color_by="capture") print(p) pdf('Windram-Monocle-pseudotime.pdf', width=fig.width, height=fig.height) print(p) dev.off()
Plot pseudotime against true capture time.
.data <- as(phenoData(cds), "data.frame") gp.monocle.pseudo <- ( ggplot(.data, aes(x=Pseudotime, y=obstime.orig)) + geom_point() + scale_x_continuous(name="pseudotime") + scale_y_continuous(name="capture time")) print(gp.monocle.pseudo) do.call(ggsave, c(list('Windram-Monocle-compare.pdf', gp.monocle.pseudo + bioinf.config), bioinf.sizes)) # What is the correlation? with(.data, cor(obstime.orig, Pseudotime, method="spearman"))
Examine the roughness under the ordering given by the Monocle pseudotime.
monocle.order <- order(pData(cds)$Pseudotime) monocle.roughness <- mean(apply(input$held.out.expr[,monocle.order], 1, calc.roughness)) print(monocle.roughness)
date()
R version and packages used:
devtools::session_info()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.