benchmark_tenx.Rmd

---
title: "Time benchmark of the different methods"
author: "Fanny Perraudeau"
date: "`r Sys.Date()`"
output: 
html_document: 
fig_height: 7
fig_width: 7
toc: yes
code_folding: hide
toc_float: yes
---

```{r options, echo=FALSE, results="hide",mesasge=FALSE, error=FALSE, include=FALSE, autodep=TRUE}
knitr::opts_chunk$set(fig.align="center", cache=TRUE, error=FALSE, message=FALSE, warning=TRUE)
```

We want to time benchmark the different DE methods.

# Functions

```{r sourceEdgeR}
runEdgeR <- function(e) {
  library(edgeR)
  condition = pData(e)$condition
  design <- model.matrix(~ condition)
  dgel <- DGEList(exprs(e))
  dgel <- edgeR::calcNormFactors(dgel)
  dgel=estimateDisp(dgel,design)
  edger.fit <- glmFit(dgel, design)
  edger.lrt <- glmLRT(edger.fit, coef="condition2")
  pvals <- edger.lrt$table$PValue
  padj <- p.adjust(pvals,method="BH")
  padj[is.na(padj)] <- 1
  logfc <- edger.lrt$table$logFC
  list(pvals=pvals, padj=padj, logfc = logfc)
}

```

```{r deseq2}
runDESeq2_poscounts <- function(e) {
  library(DESeq2)
  dds <- DESeqDataSetFromMatrix(exprs(e), colData=DataFrame(pData(e)),
                                design=~ condition)
  dds <- estimateSizeFactors(dds,type="poscounts")
  dds <- estimateDispersions(dds, minmu = 1e-3)
  dds <- nbinomWaldTest(dds, betaPrior=TRUE, modelMatrixType="standard", minmu = 1e-3)
  res <- results(dds, name="condition_2_vs_1")
  pvals <- res$pvalue
  padj <- res$padj
  logfc <- res$log2FoldChange
  list(pvals=pvals, padj=padj, logfc = logfc)
}
```

```{r zinbwaveEdgeR}
runZinbwaveEdgeR <- function(e){
  library(zinbwave)
  library(edgeR)
  condition = pData(e)$condition
  design <- model.matrix(~ condition)
  # compute zinbwave weights
  zinb <- zinbFit(exprs(e), X = design, epsilon = 1e12)
  weights <- computeObservationalWeights(zinb, exprs(e))
  d <- DGEList(exprs(e))
  d <- edgeR::calcNormFactors(d)
  d$weights <- weights
  d=estimateDisp(d, design)
  fit=glmFit(d,design)
  lrt=glmWeightedF(fit,coef=2, independentFiltering = TRUE)
  pvals = lrt$table$PValue
  list(pvals = pvals, padj = lrt$table$padjFilter,
       logfc = lrt$table$logFC)
}
```

```{r zinbwavedeseq2}
runZinbwaveDESeq2 <- function(e){
  library(zinbwave)
  library(DESeq2)
  condition = pData(e)$condition
  design <- model.matrix(~ condition)
  # compute zinbwave weights
  zinb <- zinbFit(exprs(e), X = design, epsilon = 1e12)
  weights <- computeObservationalWeights(zinb, exprs(e))
  dse <- DESeqDataSetFromMatrix(exprs(e), colData=DataFrame(pData(e)),
                                design=~ condition)
  dse <- estimateSizeFactors(dse, type="poscounts")
  dimnames(weights) = NULL
  assays(dse)[["weights"]] = weights
  dse <- estimateDispersions(dse, minmu = 1e-3)
  dse = nbinomWaldTest(dse, betaPrior=TRUE, useT=TRUE, df=rowSums(weights)-2, minmu = 1e-3)
  resultsNames(dse)
  res = results(dse, name="condition2")
  list(pvals = res$pvalue, padj = res$padj, logfc = res$log2FoldChange)
}
```

```{r}
runMAST <- function(e){
  require(MAST)
  counts <- exprs(e)
  tpm <- counts*1e6/colSums(counts)
  tpm <- log2(tpm+1)
  sca <- FromMatrix(tpm,
                    cData=data.frame(group=pData(e)$condition))
  
  # Adaptive thresholding from MAST vignette
  freq_expressed <- 0 #do not filter genes so that computation time is comparable for all the methods (same number of genes)
  thres <- thresholdSCRNACountMatrix(assay(sca), nbins = 300, conditions = pData(e)$condition)
  #par(mfrow=c(5,4))
  #plot(thres)
  assays(sca) <- list(thresh=thres$counts_threshold, tpm=assay(sca))
  expressed_genes <- freq(sca) > freq_expressed
  sca <- sca[expressed_genes,]

  ngeneson <- apply(exprs(e),2,function(x) mean(x>0))
  CD <- colData(sca)
  CD$ngeneson <- ngeneson
  CD$cngeneson <- CD$ngeneson-mean(ngeneson)
  colData(sca) <- CD
  ## differential expression
  fit <- zlm(~ cngeneson + group, sca = sca,
             method = "bayesglm", ebayes = TRUE)
  L <- matrix(0, nrow = ncol(coef(fit, "D")))
  rownames(L) <- colnames(coef(fit, "D"))
  L["group2",] <- 1
  lrFit <- lrTest(fit, hypothesis = L)
  pval <- lrFit[, 'hurdle', 'Pr(>Chisq)']
  padj <- p.adjust(pval, method = "BH")
  list(pvals = pval, padj = padj, logfc = NA)
}
```

```{r limmavoom}
runVoom <- function(e) {
  library(limma)
  condition = pData(e)$condition
  design <- model.matrix(~ condition)
  library(edgeR) # for DGEList
  dgel <- DGEList(exprs(e))
  dgel <- edgeR::calcNormFactors(dgel)
  v <- voom(dgel,design,plot=FALSE)
  fit <- lmFit(v,design)
  fit <- eBayes(fit)
  tt <- topTable(fit,coef="condition2",n=nrow(dgel),sort.by="none")
  pvals <- tt$P.Value
  padj <- p.adjust(pvals,method="BH")
  padj[is.na(padj)] <- 1
  list(pvals=pvals, padj=padj, logfc=tt$logFC)
}
```

```{r}
runDESeq2 <- function(e, retDDS=FALSE) {
  library(DESeq2)
  dds <- DESeqDataSetFromMatrix(exprs(e), DataFrame(pData(e)), ~ condition)
  dds <- DESeq(dds,betaPrior=TRUE,quiet=TRUE, minReplicatesForReplace=Inf) # imputation step disabled
  res <- results(dds)
  beta <- res$log2FoldChange
  pvals <- res$pvalue
  padj <- res$padj
  pvals[is.na(pvals)] <- 1
  padj[is.na(padj)] <- 1
  list(pvals=pvals, padj=padj, logfc=NA)
}
```

```{r scde}
runSCDE <- function(e){
  require(scde)
  # calculate models
  counts = exprs(e)
  counts<-apply(counts,2,function(x) {storage.mode(x) <- 'integer'; x})
  rownames(counts) = rownames(fData(e))
  o.ifm <- scde.error.models(counts = counts,
                             groups = pData(e)$condition, 
                             n.cores = 1,
                             threshold.segmentation = TRUE,
                             save.crossfit.plots = FALSE, 
                             save.model.plots = FALSE, 
                             verbose = 0)
  # estimate gene expression prior
  o.prior <- scde.expression.prior(models = o.ifm,
                                   counts = exprs(e),
                                   length.out = 400,
                                   show.plot = FALSE)
  # run differential expression tests on all genes.
  ediff <- scde.expression.difference(o.ifm, exprs(e), 
                                      o.prior, 
                                      groups  =  pData(e)$condition,
                                      n.randomizations  =  100,
                                      n.cores  =  1, 
                                      verbose  =  0)
  pvals = (1 - pnorm(abs(ediff$Z))) * 2
  padj = p.adjust(pvals, method = "BH")
  list(pvals = pvals, padj = padj, logfc = NA)
}
```

# Data
```{r}
load("../../realdata/data/core.rda")
tenx = core[,colData(core)$seurat %in% 1:2]
tenx = tenx[rowSums(assay(tenx)) > 0, ]
colData(tenx)$condition = factor(colData(tenx)$seurat)
tenx = tenx[sample(nrow(tenx), 10000), ]
# convert to expressionSet
tenx <- as(tenx, "ExpressionSet")
dim(tenx)
```

# Benchmark

## Time benchmarking
```{r res}
algos <- list("limma-voom"      = runVoom,
              "MAST"            = runMAST,
              "edgeR"           = runEdgeR,
              "edgeR-zinbwave"  = runZinbwaveEdgeR,
              "DESeq2"          = runDESeq2_poscounts)

algos <- list("DESeq2-zinbwave" = runZinbwaveDESeq2,
              "SCDE"            = runSCDE)
namesAlgos <- names(algos)
names(namesAlgos) <- namesAlgos
```

```{r}
tt <- lapply(namesAlgos, function(n){
  print(n)
  timing = system.time(algos[[n]](tenx))
  save(timing, file = sprintf('benchmark_tenx_%s.rda', n))
  timing
})
```

```{r}
cpuTime_tenx <- sapply(namesAlgos, function(n){
  load(sprintf('benchmark_tenx_%s.rda', n))
  timing[[3]]
})
```

## Plot

```{r cpu}
library(ggplot2)
library(cowplot)
library(RColorBrewer)

dfTime = data.frame(time = cpuTime_tenx, method = namesAlgos)
pd <- position_dodge(0.1)
cpu = ggplot(dfTime, aes(x = factor(method), y = time)) +
  geom_point(position=pd) + 
  background_grid(major = 'xy', minor = 'xy') +
  xlab('Methods') + ylab('CPU time (sec.)') +
  coord_trans(y = "log10") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(breaks = c(3, 5, 10, 25, 50, 100, 250,
                                500, 1000, 2000))
cpu
ggsave(filename="cpu.png", plot = cpu, device = 'png')
```