inst/doc/CaseStudy-SingleCellSimulation.R

In SummarizedBenchmark: Classes and methods for performing benchmark comparisons

## ----echo=FALSE, include=FALSE------------------------------------------------
knitr::opts_chunk$set(tidy = FALSE,
                      cache = TRUE,
                      dev = "png",
                      message = FALSE,
                      error = FALSE,
                      warning = TRUE)

## -----------------------------------------------------------------------------
library("SummarizedBenchmark")
library("magrittr")

## ----load-fluidigm-data-------------------------------------------------------
library("splatter")
library("scRNAseq")

fluidigm <- ReprocessedFluidigmData()
se <- fluidigm[, colData(fluidigm)[, "Coverage_Type"] == "High"]
assays(se) <- assays(se)["rsem_counts"]
assayNames(se) <- "counts"

## ----subset-data--------------------------------------------------------------
set.seed(1912)
se <- se[sample(nrow(se), 1e4), sample(ncol(se), 20)]

## ----convert-to-sce-----------------------------------------------------------
sce <- as(se, "SingleCellExperiment")

## ----construct-sim-benchdesign------------------------------------------------
bd <- BenchDesign() %>%
    addMethod(label = "splat",
              func = splatSimulate,
              params = rlang::quos(params = splatEstimate(in_data),
                                   verbose = in_verbose,
                                   seed = in_seed)) %>%
    addMethod(label = "simple",
              func = simpleSimulate,
              params = rlang::quos(params = simpleEstimate(in_data),
                                   verbose = in_verbose,
                                   seed = in_seed)) %>%
    addMethod(label = "lun",
              func = lunSimulate,
              params = rlang::quos(params = lunEstimate(in_data),
                                   verbose = in_verbose,
                                   seed = in_seed))

## ----run-sim-buildbench-------------------------------------------------------
fluidigm_dat <- list(in_data = assay(sce, "counts"),
                     in_verbose = FALSE,
                     in_seed = 19120128)
sb <- buildBench(bd, data = fluidigm_dat) 
sb

## ----check-buildbench-results-------------------------------------------------
assay(sb)
sapply(assay(sb), class)

## ----compute-sim-result-comparison--------------------------------------------
res_compare <- compareSCEs(c(ref = sce, assay(sb)[1, ]))
res_diff <- diffSCEs(c(ref = sce, assay(sb)[1, ]), ref = "ref")

## ----plot-sim-result-comparison-----------------------------------------------
res_compare$Plots$MeanVar

res_diff$Plots$MeanVar

## ----add-performance-metrics--------------------------------------------------
sb <- sb %>%
    addPerformanceMetric(
        assay = "default",
        evalMetric = "zerosPerCell",
        evalFunction = function(query, truth) {
            colMeans(assay(query[[1]], "counts") == 0)
        }) %>%
    addPerformanceMetric(
        assay = "default",
        evalMetric = "zerosPerGene",
        evalFunction = function(query, truth) {
            rowMeans(assay(query[[1]], "counts") == 0)
        } )

## ----compute-performance-metrics----------------------------------------------
sbmets <- estimatePerformanceMetrics(sb, tidy = TRUE)
sbmets  <- dplyr::select(sbmets, label, value, performanceMetric)
head(sbmets)

## ----check-performance-metrics------------------------------------------------
sbmets <- tidyr::unnest(sbmets)
head(sbmets)

## ----plot-performace-metrics--------------------------------------------------
ggplot(sbmets, aes(x = label, y = value,
                   color = label, fill = label)) +
    geom_boxplot(alpha = 1/2) +
    xlab("method") +
    scale_color_discrete("method") + 
    scale_fill_discrete("method") + 
    facet_grid(performanceMetric ~ .) +
    theme_bw()

## ----analyze-new-data---------------------------------------------------------
library(scater)
scec <- mockSCE()

scater_dat <- list(in_data = scec,
                   in_verbose = FALSE,
                   in_seed = 19120128)
buildBench(bd, data = scater_dat)