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

devtools::load_all('../..')
devtools::load_all('../../../DeLoreanData')
rmarkdown::render('Tang-DeLorean.Rmd')

library(knitr)
library(knitcitations)
library(rmarkdown)
#
# knitr options
#
opts_chunk$set(
    fig.path = 'figures/Tang-',
    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 <- read.bibtex("DeLorean.bib")
if (file.exists("config.R")) {
    source("config.R")
}
source(system.file("scripts/shared.R", package="DeLorean"))

# suppressMessages(loadfonts())
library(DeLorean)
library(DeLoreanData)
#
# Stylesheet
#
options(markdown.HTML.stylesheet = system.file(file.path('Rmd', 'foghorn.css'),
                                               package="DeLorean"))
font.family <- "Verdana"
font.theme <- theme_update(text=element_text(family=font.family))
theme_set(font.theme)

r citet(bib[["tang_tracing_2010"]]) assayed leaves at 24 time points in 2 conditions.

Data

Tang et al.'s data is available in the DeLorean R package.

library(DeLorean)
data(TangDeLorean)
dl <- de.lorean(
    tang.rna.seq,
    tang.rna.seq.gene.meta,
    tang.rna.seq.cell.meta)

Estimate hyperparameters

Examine data for empirical Bayes estimation of hyperparameters.

dl <- estimate.hyper(dl, sigma.tau=1.5)

Choose cells and genes

Select some cells at random if we have too many.

set.seed(1)
max.at.each.stage <- min(getOption("Tang.max.at.stage", nrow(dl$cell.meta)))
sample.up.to <- function(.data, size) {
    if (size < nrow(.data)) {
        sample_n(.data, size)
    } else {
        .data
    }
}
sampled.cells <- (
    dl$cell.meta
    %>% group_by(capture)
    %>% do(sample.up.to(., max.at.each.stage)))
sampled.cells

dl <- filter_cells(dl, function(cells) cells %in% sampled.cells$cell)

Select some genes at random if we have too many.

max.genes <- min(getOption("Tang.max.genes", nrow(dl$gene.meta)))
if (max.genes <= length(tang.key.genes)) {
    sampled.genes <- (dl$gene.meta
                      %>% filter(key)
                      %>% sample_n(max.genes))$gene
} else {
    sampled.genes <- sample(dl$gene.meta$gene, max.genes)
}
dl <- filter_genes(dl, function(genes) genes %in% sampled.genes)

Save expression data and meta data.

saveRDS(list(expr = dl$expr, cell.meta = dl$cell.map, gene.meta=dl$gene.map),
        file='Data/Tang-input.rds')

Format the data for Stan and fit the model.

dl <- prepare.for.stan(dl)
dl <- compile.model(dl)
dl <- find.best.tau(dl)
system.time(dl <- fit.model(dl))

Examine convergence.

dl <- examine.convergence(dl)

Analyse posterior

Examine posterior.

dl <- process.posterior(dl)
dl <- analyse.noise.levels(dl)

Profiles

Calculate expression profiles.

dl <- make.predictions(dl)

# Save DeLorean object without fit component
saveRDS({dl2 <- dl; dl2$fit <- NULL; dl2}, "Data/Tang.rds")
# dl <- readRDS("Windram.rds")

date()

R version and packages used:

sessionInfo()

JohnReid/DeLorean documentation built on Sept. 27, 2021, 5:45 a.m.