## ----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 } deseq2_lfc <- function(x) { x$log2FoldChange } 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 } edgeR_lfc <- function(x) { x$table$logFC } 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] } voom_lfc <- function(x) { x$coefficients[, 2] } bd <- bd %>% addMethod(label = "deseq2", func = deseq2_run, post = list(pv = deseq2_pv, lfc = deseq2_lfc), params = rlang::quos(countData = cntdat, colData = coldat, design = ~condition, contrast = c("condition", "2", "1"))) %>% addMethod(label = "edgeR", func = edgeR_run, post = list(pv = edgeR_pv, lfc = edgeR_lfc), params = rlang::quos(countData = cntdat, group = coldat$condition, design = ~coldat$condition)) %>% addMethod(label = "voom", func = voom_run, post = list(pv = voom_pv, lfc = voom_lfc), params = rlang::quos(countData = cntdat, group = coldat$condition, design = ~coldat$condition)) sb <- buildBench(bd, truthCols = c(pv = "status", lfc = "lfc")) ## -------------------------------------------------------------------------- updateBench(sb = sb) ## -------------------------------------------------------------------------- BDData(BenchDesign(sb)) ## -------------------------------------------------------------------------- bd2 <- modifyMethod(bd, "deseq2", rlang::quos(bd.meta = list(note = "minor update"))) ## -------------------------------------------------------------------------- updateBench(sb = sb, bd = bd2) ## -------------------------------------------------------------------------- deseq2_ashr <- function(countData, colData, design, contrast) { dds <- DESeqDataSetFromMatrix(countData, colData = colData, design = design) dds <- DESeq(dds) lfcShrink(dds, contrast = contrast, type = "ashr") } bd2 <- addMethod(bd2, "deseq2_ashr", deseq2_ashr, post = list(pv = deseq2_pv, lfc = deseq2_lfc), params = rlang::quos(countData = cntdat, colData = coldat, design = ~condition, contrast = c("condition", "2", "1"))) ## -------------------------------------------------------------------------- updateBench(sb = sb, bd = bd2) ## -------------------------------------------------------------------------- sb2 <- updateBench(sb = sb, bd = bd2, dryrun = FALSE) sb2 ## -------------------------------------------------------------------------- head(assay(sb2, "lfc")) ## -------------------------------------------------------------------------- colData(sb2)[, "session.idx", drop = FALSE] ## -------------------------------------------------------------------------- lapply(metadata(sb2)$sessions, `[[`, "methods") ## -------------------------------------------------------------------------- updateBench(sb2, bd, keepAll = FALSE) ## ----simpleMetric---------------------------------------------------------- sb <- addPerformanceMetric(sb, assay = "pv", evalMetric = "rejections", evalFunction = function(query, truth) { sum(p.adjust(query, method = "BH") < 0.1, na.rm = TRUE) }) sb <- estimatePerformanceMetrics(sb, addColData = TRUE) ## ----modifyBench----------------------------------------------------------- sb2 <- updateBench(sb, bd2, dryrun = FALSE) estimatePerformanceMetrics(sb2, rerun = FALSE, addColData = TRUE)