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