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

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

## --------------------------------------------------------------------------
data(tdat)

## --------------------------------------------------------------------------
head(tdat)

## --------------------------------------------------------------------------
adj_bonf <- p.adjust(p = tdat$pval, method = "bonferroni")

adj_bh <- p.adjust(p = tdat$pval, method = "BH")

qv <- qvalue::qvalue(p = tdat$pval)
adj_qv <- qv$qvalues

## --------------------------------------------------------------------------
adj <- cbind.data.frame(adj_bonf, adj_bh, adj_qv)
head(adj)

## --------------------------------------------------------------------------
b <- BenchDesign(tdat)

## --------------------------------------------------------------------------
b <- addMethod(bd = b, label = "bonf", func = p.adjust,
               params = rlang::quos(p = pval, method = "bonferroni"))

## --------------------------------------------------------------------------
b <- b %>% 
    addMethod(label = "BH",
              func = p.adjust,
              params = rlang::quos(p = pval, method = "BH")) %>%
    addMethod(label = "qv",
              func = qvalue::qvalue,
              params = rlang::quos(p = pval),
              post = function(x) { x$qvalues })

## --------------------------------------------------------------------------
printMethods(b)

## --------------------------------------------------------------------------
sb <- buildBench(b, truthCols = "H")

## --------------------------------------------------------------------------
head(assay(sb))

## --------------------------------------------------------------------------
colData(sb)

## --------------------------------------------------------------------------
rowData(sb)

## --------------------------------------------------------------------------
sbSub <- sb[,1:2]
colData(sbSub)

## ----addPerformanceMetric--------------------------------------------------
sb <- addPerformanceMetric(
  object = sb,
  assay = "H",
  evalMetric = "TPR",
  evalFunction = function(query, truth, alpha = 0.1) {
    goodHits <- sum((query < alpha) & truth == 1)
    goodHits / sum(truth == 1)
    }
)

performanceMetrics(sb)[["H"]]


## --------------------------------------------------------------------------
resWide <- estimatePerformanceMetrics(sb, alpha = c(0.05, 0.1, 0.2))
resWide

## ----elWide----------------------------------------------------------------
elementMetadata(resWide)

## --------------------------------------------------------------------------
sb <- estimatePerformanceMetrics(sb, 
                                 alpha = c(0.05, 0.1, 0.2), 
                                 addColData = TRUE)
colData(sb)
elementMetadata(colData(sb))

## --------------------------------------------------------------------------
estimatePerformanceMetrics(sb, 
                           alpha = c(0.05, 0.1, 0.2), 
                           tidy = TRUE)

## --------------------------------------------------------------------------
head(tidyUpMetrics(sb))

## --------------------------------------------------------------------------
tidyUpMetrics(sb) %>%
  dplyr:::filter(label == "bonf", alpha == 0.1, performanceMetric == "TPR") %>%
  dplyr:::select(value)

## --------------------------------------------------------------------------
library("limma")
library("edgeR")
library("DESeq2")
library("tximport")

## ----loadingSoneson--------------------------------------------------------
data("soneson2016")
head( txi$counts )
head( truthdat )

## --------------------------------------------------------------------------
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(mydat)

## --------------------------------------------------------------------------
deseq2_pvals <- function(countData, colData, design, contrast) {
    dds <- DESeqDataSetFromMatrix(countData,
                                  colData = colData,
                                  design = design)
    dds <- DESeq(dds)
    res <- results(dds, contrast = contrast)
    res$pvalue
}

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

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

## --------------------------------------------------------------------------
bd <- bd %>%
    addMethod(label = "deseq2",
              func = deseq2_pvals,
              params = rlang::quos(countData = cntdat,
                                   colData = coldat, 
                                   design = ~condition,
                                   contrast = c("condition", "2", "1"))
              ) %>%
    addMethod(label = "edgeR",
              func = edgeR_pvals,
              params = rlang::quos(countData = cntdat,
                                   group = coldat$condition,
                                   design = ~coldat$condition)
              ) %>%
    addMethod(label = "voom",
              func = voom_pvals,
              params = rlang::quos(countData = cntdat,
                                   group = coldat$condition,
                                   design = ~coldat$condition)
              )

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

## --------------------------------------------------------------------------
sb

## ----availableMetrics------------------------------------------------------
availableMetrics()

## --------------------------------------------------------------------------
sb <- addPerformanceMetric( sb, 
                            evalMetric=c("rejections", "TPR", "TNR", "FDR", "FNR"),
                            assay="status" )
names(performanceMetrics(sb)[["status"]])

## ----echo=FALSE------------------------------------------------------------
assay(sb)[,"deseq2"][is.na(assay(sb)[, "deseq2"])] <- 1

## --------------------------------------------------------------------------
estimatePerformanceMetrics(
  sb, 
  alpha = c(0.01, 0.05, 0.1, 0.2), 
  tidy = TRUE) %>%
  dplyr:::select(label, value, performanceMetric, alpha) %>%
  tail()

## ---- fig.width=4.5, fig.height=4------------------------------------------
plotMethodsOverlap( sb, assay="status", alpha=0.1, order.by="freq")

## ---- fig.width=5, fig.height=4--------------------------------------------
SummarizedBenchmark::plotROC(sb, assay="status")

## --------------------------------------------------------------------------
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 <- BenchDesign(mydat) %>%
    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 = bd, truthCols = c(pv = "status", lfc = "lfc"))
sb

## --------------------------------------------------------------------------
head(assay(sb, "pv"))
head(assay(sb, "lfc"))

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

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

## --------------------------------------------------------------------------
BenchDesign(mydat) %>%
    addMethod(label = "deseq2",
              func = deseq2_pvals,
              meta = list(reason = "recommended by friend X"),
              params = rlang::quos(countData = cntdat,
                                   colData = coldat, 
                                   design = ~condition,
                                   contrast = c("condition", "2", "1"))
              ) %>%
    addMethod(label = "edgeR",
              func = edgeR_pvals,
              meta = list(reason = "recommended by friend Y"), 
              params = rlang::quos(countData = cntdat,
                                   group = coldat$condition,
                                   design = ~coldat$condition)
              ) %>%
    buildBench() %>%
    colData()

## --------------------------------------------------------------------------
BenchDesign(mydat) %>%
    addMethod(label = "deseq2",
              func = deseq2_pvals,
              meta = list(pkg_func = rlang::quo(DESeq2::DESeq)),
              params = rlang::quos(countData = cntdat,
                                   colData = coldat, 
                                   design = ~condition,
                                   contrast = c("condition", "2", "1"))
              ) %>%
    addMethod(label = "edgeR",
              func = edgeR_pvals,
              meta = list(pkg_name = "edgeR",
                          pkg_vers = as.character(packageVersion("edgeR"))), 
              params = rlang::quos(countData = cntdat,
                                   group = coldat$condition,
                                   design = ~coldat$condition)
              ) %>%
    buildBench() %>%
    colData()

## --------------------------------------------------------------------------
printMethod(bd, "deseq2")

## --------------------------------------------------------------------------
modifyMethod(bd, label = "deseq2", 
             params = rlang::quos(contrast = c("condition", "1", "2"),
                                  bd.meta = list(note = "modified post hoc"))
             ) %>%
    printMethod("deseq2")

## --------------------------------------------------------------------------
modifyMethod(bd, label = "deseq2", 
             params = rlang::quos(contrast = c("condition", "1", "2"),
                                  bd.meta = list(note = "modified post hoc")),
             .overwrite = TRUE) %>%
    printMethod("deseq2")

## --------------------------------------------------------------------------
bde <- expandMethod(bd, label = "deseq2", 
                    onlyone = "contrast",
                    params = rlang::quos(deseq2_v1 = c("condition", "1", "2"),
                                         deseq2_v2 = c("condition", "2", "2"))
                    )
printMethod(bde, "deseq2_v1")
printMethod(bde, "deseq2_v2")

## --------------------------------------------------------------------------
bde <- expandMethod(bd, label = "deseq2", 
                    params = list(deseq2_v1 = rlang::quos(contrast = c("condition", "1", "2"),
                                                          meta = list(note = "filp order")),
                                  deseq2_v2 = rlang::quos(contrast = c("condition", "2", "2"),
                                                          meta = list(note = "nonsensical order"))
                                  )
                    )
printMethod(bde, "deseq2_v1")
printMethod(bde, "deseq2_v2")

## --------------------------------------------------------------------------
dropMethod(bd, "deseq2") %>%
    printMethods()

## --------------------------------------------------------------------------
bdnull <- BenchDesign()
bdnull

## --------------------------------------------------------------------------
bdnull <- bdnull %>%
    addMethod(label = "bonf",
              func = p.adjust,
              params = rlang::quos(p = pval,
                                   method = "bonferroni")) %>%
    addMethod(label = "BH",
              func = p.adjust,
              params = rlang::quos(p = pval,
                                   method = "BH"))

## --------------------------------------------------------------------------
buildBench(bdnull, data = tdat)