In jhsiao999/ashbun: Methods and evaluation tools for single cell RNA-seq count data

knitr::opts_chunk$set(echo = TRUE)

Background and goals

This document was created for testing the ashbun package on four UMI datasets. The goal is to run the codes without error.

---

Analysis

Load packages

library(Biobase)
library(ashbun)
library(singleCellRNASeqMouseZeiselBrain)
#library(singleCellRNASeqHumanTungiPSC)
#library(singleCellRNASeqMouseKleinESC)

Test for one simulated single cell data.

eset <- get(data("MouseZeiselBrain"))
counts <- exprs(eset)

samples_any_gene <- which( colSums(counts) != 0)
counts <- counts[,samples_any_gene]

genes_to_include <- which(apply(counts,1,sum)>0)
counts <- counts[genes_to_include, ]
dim(counts)


#---- generat simulated datasets
simdata_list <- simulationWrapper(counts, Nsim = 1,
                                  Ngenes = 1000,
                                  Nsamples = 50,
                                  sample_method = "all_genes",
                                  pi0 = .9,
                                  beta_args = args.big_normal(betapi = 1,
                                                              betamu = 0, betasd = .8))

simdata <- simdata_list[[1]]

library(doParallel)
cl <- makeCluster(8)
registerDoParallel(cl)

# ---- gather evaluation results
output  <- query.evaluation(counts = simdata$counts,
                            condition = simdata$condition,
                            is_nullgene = simdata$is_nullgene,
                            methodsNormalize = c("LIB", "TMM", "census","scran"),
                            methodsMeanExpression = c("BPSC", "MAST", "ROTS", "SCDE",
                                                      "DESeq2", "edgeR", "limmaVoom"))
results.one.data <- saveRDS(output, file = "tests.results.one.data.rds")

results.one.data <- readRDS(file = "tests.results.one.data.rds")
library(ggplot2)
qplot(x = FPR, y = TPR, data = results.one.data$roc, 
      colour = methodsMeanExpression, 
      facets = ~ methodsNormalize,
      geom = "path")

Run on 2 simulated data.

simdata_list <- simulationWrapper(counts, Nsim = 2,
                                  Ngenes = 1000,
                                  Nsamples = 50,
                                  sample_method = "all_genes",
                                  pi0 = .9,
                                  beta_args = args.big_normal(betapi = 1,
                                                              betamu = 0, betasd = .8))

output <- vector("list", 2)
for (index in 1:length(output)) {
  output[[index]]  <- query.evaluation(counts = simdata_list[[index]]$counts,
                            condition = simdata_list[[index]]$condition,
                            is_nullgene = simdata_list[[index]]$is_nullgene,
                            methodsNormalize = c("LIB", "TMM", "census","scran"),
                            methodsMeanExpression = c("BPSC", "MAST", "ROTS", "SCDE",
                                                      "DESeq2", "edgeR", "limmaVoom"),
                            nsim = index)
}
results.multiple.data <- saveRDS(output, file = "tests.results.multiple.data.rds")

results.multiple.data <- readRDS(file = "tests.results.multiple.data.rds")

# summarize ROC results
roc_combined <- do.call(rbind, lapply(results.multiple.data, "[[", 2))
roc_average <- do.call(rbind, lapply(1:nlevels(roc_combined$methodsNormalize),
                                     function(index_methodsNormalize) {
    normMethod <- levels(roc_combined$methodsNormalize)[index_methodsNormalize]
    foo2 <- 
      do.call(rbind, lapply(1:nlevels(roc_combined$methodsMeanExpression), 
                            function(index_methodsMeanExpression) {
              meanMethod <- levels(roc_combined$methodsMeanExpression)[index_methodsMeanExpression]
              roc_combined_subset <- subset(roc_combined, 
                  subset = methodsMeanExpression %in%  meanMethod & methodsNormalize %in% normMethod )
              foo <- getROC.average(roc_combined_subset) 
              foo$methodsMeanExpression <- meanMethod 
              return(foo) 
              }) )
    foo2$methodsNormalize <- normMethod
    return(foo2)
  }) )

library(ggplot2)
qplot(data = roc_average, 
      x = FPR, y = TPR, 
      colour = methodsMeanExpression, 
      facets = ~ methodsNormalize,
      geom = "path")

---

Session information

sessionInfo()

jhsiao999/ashbun documentation built on May 8, 2019, 11:17 p.m.