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