inst/doc/SummarizedBenchmark-FullCaseStudy.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")

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(data = mydat)

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

## -----------------------------------------------------------------------------
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_lfc <- function(x) {
    x$log2FoldChange
}

edgeR_lfc <- function(x) {
    x$table$logFC
}

voom_lfc <- function(x) {
    x$coefficients[, 2]
}

## -----------------------------------------------------------------------------
bd <- BenchDesign(data = 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"))

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

## ----download-txi, eval = FALSE-----------------------------------------------
#  library(tximport)
#  library(readr)
#  
#  d <- tempdir()
#  download.file(url = paste0("https://www.ebi.ac.uk/arrayexpress/files/",
#                             "E-MTAB-4119/E-MTAB-4119.processed.3.zip"),
#                destfile = file.path(d, "samples.zip"))
#  unzip(file.path(d, "samples.zip"), exdir = d)
#  
#  fl <- list.files(d, pattern = "*_rsem.txt", full.names=TRUE)
#  names(fl) <- gsub("sample(.*)_rsem.txt", "\\1", basename(fl))
#  
#  txi <- tximport(fl, txIn = TRUE, txOut = TRUE,
#                  geneIdCol = "gene_id",
#                  txIdCol = "transcript_id",
#                  countsCol = "expected_count",
#                  lengthCol = "effective_length",
#                  abundanceCol = "TPM",
#                  countsFromAbundance = "scaledTPM",
#                  importer = function(x){ readr::read_tsv(x) })

## ----download-truthdat, eval = FALSE------------------------------------------
#  download.file(url = paste0("https://www.ebi.ac.uk/arrayexpress/files/",
#                             "E-MTAB-4119/E-MTAB-4119.processed.2.zip"),
#                destfile = file.path(d, "truth.zip"))
#  unzip(file.path(d, "truth.zip"), exdir = d)
#  
#  truthdat <- read_tsv(file.path(d, "truth_transcript.txt"))

## ----save-rda, eval = FALSE---------------------------------------------------
#  save(txi, truthdat, file = "../data/soneson2016.rda")