inst/doc/Feature-Parallel.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)

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

## ----case-study-data----------------------------------------------------------
library("limma")
library("edgeR")
library("DESeq2")
library("tximport")

data(soneson2016)

mycounts <- round(txi$counts)
mycoldat <- data.frame(condition = factor(rep(c(1, 2), each = 3)))
rownames(mycoldat) <- colnames(mycounts)
mydat <- list(coldat = mycoldat, cntdat = mycounts,
              status = truthdat$status, lfc = truthdat$logFC)

bd <- BenchDesign(data = mydat)

## ----case-study-methods-------------------------------------------------------
deseq2_run <- function(countData, colData, design, contrast) {
    dds <- DESeqDataSetFromMatrix(countData, colData = colData, design = design)
    dds <- DESeq(dds)
    results(dds, contrast = contrast)
}
deseq2_pv <- function(x) { x$pvalue }

edgeR_run <- function(countData, group, design) {
    y <- DGEList(countData, group = group)
    y <- calcNormFactors(y)
    des <- model.matrix(design)
    y <- estimateDisp(y, des)
    fit <- glmFit(y, des)
    glmLRT(fit, coef=2)
}
edgeR_pv <- function(x) { x$table$PValue }

voom_run <- function(countData, group, design) {
    y <- DGEList(countData, group = group)
    y <- calcNormFactors(y)
    des <- model.matrix(design)
    y <- voom(y, des)
    eBayes(lmFit(y, des))
}
voom_pv <- function(x) { x$p.value[, 2] }

bd <- bd %>%
    addMethod(label = "deseq2", func = deseq2_run, post = deseq2_pv,
              params = rlang::quos(countData = cntdat,
                                   colData = coldat, 
                                   design = ~condition,
                                   contrast = c("condition", "2", "1"))) %>%
    addMethod(label = "edgeR", func = edgeR_run, post = edgeR_pv,
              params = rlang::quos(countData = cntdat,
                                   group = coldat$condition,
                                   design = ~coldat$condition)) %>%
    addMethod(label = "voom", func = voom_run, post = voom_pv, 
              params = rlang::quos(countData = cntdat,
                                   group = coldat$condition,
                                   design = ~coldat$condition))

## -----------------------------------------------------------------------------
bpparam()
sbp <- buildBench(bd, parallel = TRUE,
                  BPPARAM = BiocParallel::SerialParam())
sbp

## -----------------------------------------------------------------------------
sb <- buildBench(bd, truthCols = "status")

## -----------------------------------------------------------------------------
all(assay(sbp) == assay(sb), na.rm = TRUE)

## -----------------------------------------------------------------------------
ndat <- length(mydat$status)
spIndexes <- split(seq_len(ndat), rep( 1:3, length.out = ndat))

datasetList <- lapply(spIndexes, function(x) {
    list(coldat = mydat$coldat, 
        cntdat = mydat$cntdat[x, ], 
        status = mydat$status[x], 
        lfc = mydat$lfc[x])
})

names(datasetList) <- c("dataset1", "dataset2", "dataset3")

## -----------------------------------------------------------------------------
sbL <- bplapply(datasetList, function(x) { buildBench(bd, data = x) },
                BPPARAM = MulticoreParam(3))
sbL