In JohnReid/DeLorean: Estimates Pseudotimes for Single Cell Expression Data

library(rmarkdown)
render("generate-figures.Rmd")

# library(devtools)
# load_all("/home/john/Dev/DeLorean")
# suppressMessages(loadfonts())
# library(devtools)
# load_all("/home/john/Dev/DeLorean")
# library(DeLorean)
# library(matrixcalc)
library(MASS)
library(ggplot2)
library(dplyr)
library(reshape2)
library(stringr)
# library(grDevices)

pts.in.inch <- 72
golden.ratio <- (1 + sqrt(5)) / 2
page.width <- 5
page.height <- page.width / golden.ratio
theme.present <- theme_grey(base_size=7)
mrc.colors <- c(
    rgb(138, 121, 103, maxColorValue=255),
    rgb(217, 165, 41 , maxColorValue=255),
    rgb(153, 152, 40 , maxColorValue=255),
    rgb(117, 139, 121, maxColorValue=255),
    rgb(33 , 103, 126, maxColorValue=255),
    rgb(208, 114, 50 , maxColorValue=255),
    rgb(106, 59 , 119, maxColorValue=255),
    rgb(130, 47 , 90 , maxColorValue=255)
)
scale.color.mrc <- scale_color_manual(values = mrc.colors)
scale.fill.mrc <- scale_fill_manual(values = mrc.colors)

set.seed(1)
tmin <- 1
tmax <- 4
twidth <- tmax - tmin
tlims <- c(tmin, tmax)
#
# Set up cells
#
capture.times <- seq(tmin, tmax, length=twidth+1)
num.cells <- 40
sigma.tau <- .5
cells <- data.frame(
    cell = 1:num.cells,
    capture = sample(capture.times, num.cells, replace=TRUE)
) %>% mutate(pseudotime=rnorm(num.cells, mean=capture, sd=sigma.tau))
#
# Set up genes
#
num.genes <- 12
gene.levels <- str_c('Gene.', 1:num.genes)
genes <- data.frame(
    gene=factor(gene.levels, levels=gene.levels, ordered=TRUE)
)
#
# Set up covariance function
#
length.scale <- twidth / 3
cov.fn <- function(r) exp(-.5*r**2)
calc.cov <- function(input.A, input.B=input.A) {
    # Calculate distances between inputs
    r <- outer(input.A, input.B, FUN="-")
    # Apply covariance function to distances scaled by length.scale
    cov.fn(r/length.scale)
}
noise <- .2
K00 <- calc.cov(cells$pseudotime) + noise**2 * diag(num.cells)
#
# Sample genes
#
expr.mat <- mvrnorm(num.genes, rep(3, num.cells), K00)
class(expr.mat)
dimnames(expr.mat) <- list(gene=genes$gene, cell=cells$cell)
expr <- (
    melt(expr.mat, value.name="expression")
    %>% mutate(gene=factor(gene, levels=gene.levels, ordered=TRUE))
    %>% left_join(genes)
    %>% left_join(cells)
)

Make plots

scale.color.expr <- scale_color_brewer(name="capture\ntime", palette="Set1")
gp.capture.single <- (
    ggplot(
        expr %>% filter("Gene.3" == gene),
        aes(x=capture, y=expression, color=factor(capture))
    )
    + geom_point(shape=3, size=1)
    + scale.color.expr
)
ggsave('Figures/capture-single.pdf',
       gp.capture.single + theme.present + theme(legend.position="none"),
       width=page.width/2, height=page.height/2)
gp.capture <- (
    ggplot(
        expr,
        aes(x=capture, y=expression, color=factor(capture))
    )
    + geom_point(shape=3, size=1)
    + scale.color.expr
    + facet_wrap(~ gene)
)
ggsave('Figures/capture-multi.pdf', gp.capture + theme.present,
       width=page.width, height=page.height)
print(gp.capture)
gp.pseudo.single <- (
    ggplot(
        expr %>% filter("Gene.3" == gene),
        aes(x=pseudotime, y=expression, color=factor(capture))
    )
    + geom_point(shape=3, size=1)
    + scale.color.expr
)
ggsave('Figures/pseudo-single.pdf',
       gp.pseudo.single + theme.present + theme(legend.position="none"),
       width=page.width/2, height=page.height/2)
gp.pseudo <- (
    ggplot(
        expr,
        aes(x=pseudotime, y=expression, color=factor(capture))
    )
    # + geom_point(aes(x=capture))
    + geom_point(shape=3, size=1)
    + scale.color.expr
    + facet_wrap(~ gene)
)
ggsave('Figures/pseudo-multi.pdf', gp.pseudo + theme.present,
       width=page.width, height=page.height)
print(gp.pseudo)