benchmark_usoskin.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
  pickingSession = pData(e)$pickingSession
  design <- model.matrix(~ condition + pickingSession)
  dgel <- DGEList(exprs(e))
  dgel <- edgeR::calcNormFactors(dgel)
  dgel=estimateDisp(dgel,design)
  edger.fit <- glmFit(dgel, design)
  edger.lrt <- glmLRT(edger.fit, coef="conditionB")
  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 + pickingSession)
  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_B_vs_A")
  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
  pickingSession = pData(e)$pickingSession
  design <- model.matrix(~ condition + pickingSession)
  # 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
  pickingSession = pData(e)$pickingSession
  design <- model.matrix(~ condition + pickingSession)
  # 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 + pickingSession)
  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)
  res = results(dse, name="conditionB")
  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,
                                     pickingSession=pData(e)$pickingSession))
  
  # 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 = 10, min_per_bin = 50,
                                     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 + pickingSession, sca = sca,
             method = "bayesglm", ebayes = TRUE)
  L <- matrix(0, nrow = ncol(coef(fit, "D")))
  rownames(L) <- colnames(coef(fit, "D"))
  L["groupB",] <- 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
  pickingSession = pData(e)$pickingSession
  design <- model.matrix(~ condition + pickingSession)
  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="conditionB",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)
  rownames(counts) = fData(e)[,1]
  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,
                                      batch=factor(pData(e)$pickingSession))
  pvals = (1 - pnorm(abs(ediff[[1]]$Z))) * 2
  padj = p.adjust(pvals, method = "BH")
  list(pvals = pvals, padj = padj, logfc = NA)
}
```

# Data

```{r data}
set.seed(7828)
library(zinbwave)
path = '../../datasets/usoskin/'
load(paste0(path, 'esetUsoskin.RData'))
eset=eset[rowSums(exprs(eset)>0)>=20,]
exprs(eset) <- apply(exprs(eset),2,function(x) {storage.mode(x) <- 'integer'; x})
file = paste0(path, "subsetMatrixUsoskinFPR_randomCellTypes.txt")
subsets <- read.table(file)
cData = pData(eset)[,c('Level 3', 'Picking sessions')]
colnames(cData) = c('Level3', 'Pickingsessions')

usoskin <- eset[,as.numeric(subsets[1,])]
usoskin <- usoskin[rowSums(exprs(usoskin)) != 0, ]
usoskin <- usoskin[sample(nrow(usoskin), 10000), ]
condition=factor(rep(c("A","B"),each=45))
pickingSession=factor(rep(rep(c("Cold","rt1","rt2"),each=15),2))
pData(usoskin)$condition=condition
pData(usoskin)$pickingSession=pickingSession
dim(usoskin)
```

# Benchmark

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

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

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

## Plot

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

dfTime = data.frame(time = cpuTime_usoskin,
                    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')
```