## ----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