---
title: "Comparison Of DEgage Against Other Methods"
output: 
    rmarkdown::html_document:
        number_sections: true
        toc: true
        toc_depth: 2
        theme: cerulean
       
        
vignette: >
  %\VignetteIndexEntry{my-vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = FALSE,
  comment = "#>"
)
```

# Introduction

---

In the original DEGage publication [(1)](#ref1), DEGage was compared against multiple popular differential expression analysis tools and tests. This vignette details the general workflow used for these comparisions with a dummy example on simulated data. 

First, all of the appropriate libraries must be loaded in: 
```{r message = FALSE}
library(DEGage)
library(DESeq2)
library(edgeR)
library(DEsingle)
library(scDD)
library(monocle3)
library(MASS)
library(dplyr)
library(SingleCellExperiment)
library(ggplot2)
library(ggpubr)
```
# Count Simulation

---

For the purposes of this vignette, the packages will be compared on a small data set generated with `DEGage_Simulation()`. This function simulates counts along a negative binomial simulation and simulates differential expression by introducing a log2 fold change difference between the means of two sets of counts. The means, dispersions, and dropout proportions can be generated automatically by the simulation framework; however, in this case, we will control the dispersions, dropout proportions, and log2 fold change differences.
```{r}
ngenes = 800
ndegs = 200
ngroup1 = 20
ngroup2 = 20
lfcs = runif(ndegs, min = 2.5, max = 7.5)
dispersions = runif(ndegs + ngenes, min = 0.1, max = 10)
means = abs(rnorm(ndegs + ngenes, mean = 50, sd = 20))

counts <- DEGage_Simulation(ngenes = ngenes,
                            ndegs = ndegs,
                            ngroup1 = ngroup1,
                            ngroup2 = ngroup2,
                            lfc = lfcs,
                            dispersions = dispersions,
                            means = means ,
                            max.prop.zeros = 0.4,
                            min.prop.zeros = 0.2)
head(counts)
```
# Package Usage

---

In this section, 
## DEGage
DEGage detects pairwise differential expression by using the difference of two negative binomial (DOTNB) distribution. It uses a generalized linear model (glm) to estimate the r and p parameters of the count distribution for two conditions independently. Then, it tests for differential expression using the closed form of the DOTNB cumulative density function. DEGage only has two required inputs. The first is a `data.frame` containing counts, where columns represent cells and rows represent genes. The second is is a `factor` containing a vector of counts that indicates which cell belongs to which group. The entire DEGage pipeline is encapsulated within a single command, `DEGage()`. 
```{r message = FALSE }
groups <- factor(c(rep(1, ngroup1), rep(2, ngroup2)))
start.time <- Sys.time()
degage.results <- DEGage(counts, groups, perm.pval = 0.05)
degage.runtime <- difftime(Sys.time(), start.time, units = "secs")
```
  
DEGage outputs the estimated parameters for both gene wise distributions, pvalues,and FDRs for each gene:
```{r}
head(degage.results)
```
## DESeq2
DESeq2 [(2)](#ref2) fits negative binomial parameters for each gene using a glm, then, it uses empirical Bayes shrinkage to moderate gene wise dispersions towards a curve fit between the means and dispersions of all genes in a data set. Then, it uses a Wald test to determine significant differential expression. For a more detailed account of DESeq2 usage, please see [this Bioconductor tutorial](https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)
```{r message = FALSE, warning = FALSE}
dim.matrix <- matrix(c(rep("Group1", ngroup1), rep("Group2", ngroup2)),
                     nrow = ngroup1+ngroup2,
                     dimnames = list(colnames(counts), 'Group'))
start.time <- Sys.time()
DESeqobj <- DESeqDataSetFromMatrix(countData = counts+1,
                                   colData = dim.matrix,
                                   design = ~Group)
DESeq2.output <-DESeq(DESeqobj)
DESeq2.runtime <- difftime(Sys.time(), start.time, units = "secs")
DESeq2.results <- as.data.frame(DESeq2::results(DESeq2.output))
```
  
DESeq2 outputs shrunken log2 fold changes, alogside Wald Test p values and adjusted p values:
```{r}
head(DESeq2.results)
```
## edgeR
edgeR [(3)](#ref3) uses maxiumum likelihood to fit negative binomial parameters, then uses a similar empirical Bayes procedure to DESeq2 to shrink dispersions towards a consensus estimate. It then uses a test analagous to a Fishers Exact test to determine differential expression. For extended edgeR documentation, see [the edgeR user guide](https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf)
```{r message = FALSE}
dim.matrix <- matrix(c(rep("Group1", ngroup1), rep("Group2", ngroup2)),
                     nrow = ngroup1+ngroup2,
                     dimnames = list(colnames(counts), 'Group'))
dim.matrix <- factor(dim.matrix)

list <-DGEList(counts)

design <- model.matrix(~0+dim.matrix)
colnames(design) <- levels (dim.matrix)

start.time <- Sys.time()
AveLogCPM <-aveLogCPM(list)
list <- calcNormFactors(list)
list <- estimateDisp(list, design, Robust = TRUE)
fit <-glmQLFit(list, design, robust = TRUE)
onev.two <-makeContrasts(Group1-Group2, levels = design)
res <- glmQLFTest(fit, contrast = onev.two)
edgeR.runtime <- difftime(Sys.time(), start.time, units = "secs")

edgeR.results <- as.data.frame(topTags(res, n =nrow(counts) ,adjust.method = "fdr",p.value = 0.05))
```
  
Similarly to DESeq2, edgeR provides shrunken log2 fold change estimates, an F statistic, a p value, and an adjusted p value for each gene:
```{r}
head(edgeR.results)
```

## DESingle
DESingle [(4)](#ref4) uses maximum likelihood estimation to fit gene wise parameters along a zero inflated negative binomial (ZINB) distribution. Then, it uses a likelihood ratio test to determine differential expression. By using information provided by the third dispersion parameter for the ZINB distrbution, DESingle classifies DE genes into three categories: differential expression (DE), differential expression abundance (DEa), and general differential expression (DEg). For more information on DESingle, see [its documentation](https://bioconductor.org/packages/devel/bioc/manuals/DEsingle/man/DEsingle.pdf). 
```{r message = FALSE, warning = FALSE, results = "hide"}
groups <- factor(c(rep(1, ngroup1), rep(2, ngroup2)))
start.time <- Sys.time()
DEsingle.output <- DEsingle(counts, groups)
DESingle.runtime <- difftime(Sys.time(), start.time, units = "secs")
DEsingle.results <- DEtype(DEsingle.output, threshold = 0.05)
```
DEsingle returns gene parameters, p values, adjusted pvalues, DE classifications, and information as to whether a gene was up or down regulated:
```{r}
head(DEsingle.results)
```

## monocle3
This monocle3 [(5)](#ref5) workflow fits a generalized linear model for each gene, and then it uses a Wald Test to determine differential expression. monocle3 provides multiple distribution families to fit the counts to; however, a negative binomial glm was chosen in this case, given that `DEGage_Simulation()` simulates counts along an NB distribution. See the [monocle3 differential expression analysis tutorial](https://cole-trapnell-lab.github.io/monocle3/docs/differential/) for more information. 
```{r message = FALSE, warning = FALSE}
cell_metadata = data.frame(Group = c(rep(1,ngroup1), rep(2,ngroup2)),
                           row.names = colnames(counts))
gene_metadata = data.frame(gene_short_name = rownames(counts),
                           row.names = rownames(counts))

start.time <- Sys.time()
cds <- new_cell_data_set(data.matrix(counts),
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_metadata)
gene_fits <- fit_models(cds,
                        model_formula_str = "~Group",
                        expression_family = "negbinomial")
fit_coefs <- coefficient_table(gene_fits)
monocle3.runtime <- difftime(Sys.time(), start.time, units = "secs")

fit_coefs
```
fit_coefs contains parameters and p values for each gene. For simplicity, we will select the gene names and the p values from fit_coefs: 
```{r}
intermediate <- fit_coefs %>% filter(term == "Group")

monocle.results <-intermediate %>%
  filter (p_value < 0.05) %>%
  select(gene_short_name,p_value)
monocle.results <- monocle.results %>%
  distinct(monocle.results$gene_short_name, .keep_all = TRUE)

monocle.results <- as.data.frame(monocle.results)
rownames(monocle.results) <- monocle.results$gene_short_name
monocle.results <- monocle.results[,c(1,2)]
monocle.results$p.adj <- p.adjust(monocle.results$p_value, method = "fdr")

head(monocle.results)
```
## scDD
scDD [(6)](#ref6) models genes using a Dirichlet process mixture. It classifies genes into four categories, DE (differential expression), DP (differential proportions), DM (differential modalities), and DB (combined DP + DM). A vignette detailing scDD usage is available [here](https://bioconductor.org/packages/devel/bioc/vignettes/scDD/inst/doc/scDD.pdf).
```{r warning = FALSE, message = FALSE}
sce <- SingleCellExperiment(assays = list(counts = counts))
colData(sce)$condition = c(rep("Group1", ngroup1), rep("Group2", ngroup2))
start.time <- Sys.time()
sce <- preprocess(sce, scran_norm = TRUE) #counts must be normalized
scDD.output <- scDD(sce)
scDD.runtime <- difftime(Sys.time(), start.time, units = "secs")
scDD.results <- as.data.frame(scDD::results(scDD.output))
```
scDD returns a categorization and p values for each gene:
```{r}
head(scDD.results)
```

## Wilcoxon
Although it's not a specific scRNA-seq analysis framework, Wilcoxon tests are frequently used to test for differential expression, and are the default test used in large frameworks like [Scanpy](https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering.html). The Wilcoxon rank sum test and non parametric and operates under the null hypothesis that two sets of counts follow a continuous distribution. 
```{r warning = FALSE, message = FALSE}
groups <- c(rep(1, ngroup1), rep(2, ngroup2))

calc.wilcox <- function(genewise.counts, groups){
  g1 <- as.numeric(genewise.counts[groups == 1])
  g2 <- as.numeric(genewise.counts[groups == 2])
  t <-wilcox.test(g1, g2)
  return(t$p.value)
}

start.time <- Sys.time()
wilcox.output <- apply(X = counts, FUN = calc.wilcox, MARGIN = 1,groups = groups)
wilcox.runtime <- difftime(Sys.time(), start.time, units = "secs")
wilcox.results <- data.frame(
                    pval = as.numeric(wilcox.output),
                    p.adj = p.adjust(as.numeric(wilcox.output), method = "fdr"),
                    row.names = names(wilcox.output)
                  )
head(wilcox.results)
```

## glm.nb
In addition to using a wilcoxon test, Petrany et. al also tested the performance of the `glm.nb` function from the `MASS` package. Some of the most popular DE analysis tools use a glm to estimate NB parameters, and then perform some form of post-processing on the parameters, such as the Bayes shrinkage strategy used by DESeq2. The following code runs a NB glm, and then performs a likelihood ratio test to determine differential expression. This tests the capabilites of glm's in DE testing alone, without the bells and whistles of other DE analysis tools. 
```{r warning = FALSE, message = FALSE}
groups <- c(rep(1, ngroup1), rep(2, ngroup2))

calc.glm <- function(genewise.counts, groups){
  design <- data.frame(counts = genewise.counts,
                         group = groups)
    res <- tryCatch({res<-glm.nb(counts~group, design)},
                    error = function(e){return(NA)})
    if(length(res) == 1){
      p <- NA
    }else{
      lrt.res <- anova(res)
      p <- lrt.res$`Pr(>Chi)`[2]
    }
    return(p)
}

start.time <- Sys.time()
glm.output <- apply(X = counts, 
                    FUN = calc.glm, 
                    MARGIN = 1, 
                    groups = groups)
glm.runtime <- difftime(Sys.time(), start.time, units = "secs")
glm.results <- data.frame(
                pval = as.numeric(glm.output),
                p.adj = p.adjust(as.numeric(glm.output), method = "fdr"),
                row.names = names(glm.output)
              )
head(glm.results)
```

# Method Comparison

---

Next, we will perform a brief comparison of all the methods detailed above. First, we will ensure all of the dataframes containing the results of each package have a column called `p.adj`, if they don't already. 
``` {r}
colnames(degage.results)[15] <- "p.adj"
colnames(DESeq2.results)[6] <- "p.adj"
colnames(edgeR.results)[5] <- "p.adj"
colnames(DEsingle.results)[21] <- "p.adj"
colnames(scDD.results)[11] <- "p.adj"
```
Next, we will obtain basic performace statistics for each framework: 
```{r}
get_performance <- function(results, counts, packagename){
  
  results<- results[!is.na(results$p.adj),]
  expDEs<- rownames(results[results$p.adj <= 0.05,])
  
  trueDEs <- rownames(counts[1:ndegs,])
  
  trueEEs <- rownames(counts[ndegs + 1:nrow(counts),])
  
  tp <- sum(trueDEs %in% expDEs)
  fp <- sum(trueEEs %in% expDEs)
  tn <- sum(!(trueEEs %in% expDEs))
  fn <- sum(!(trueDEs %in% expDEs))
  
  sensitivity <- tp/(tp+fn)
  specificity <- tn/(tn+fp)
  f1 <- (2*tp)/(2*tp+fp+fn)
  
  performance <- data.frame("Package" = packagename, 
                            "ndegs"= length(expDEs), 
                            "Sensitivity" = sensitivity, 
                            "Specificity" = specificity, 
                            "F1" = f1)
  return(performance)
}


performance <- rbind(get_performance(degage.results, counts, "DEGage"), 
                      get_performance(DESeq2.results, counts, "DESeq2"),
                      get_performance(edgeR.results, counts, "edgeR"), 
                      get_performance(DEsingle.results, counts, "DEsingle"), 
                      get_performance(monocle.results, counts, "monocle3"), 
                      get_performance(scDD.results, counts, "scDD"), 
                      get_performance(wilcox.results, counts, "Wilcoxon"), 
                      get_performance(glm.results, counts, "glm.nb"))

performance$Runtime <- as.numeric(c(degage.runtime, 
                         DESeq2.runtime, 
                         edgeR.runtime, 
                         DESingle.runtime, 
                         monocle3.runtime, 
                         scDD.runtime, 
                         wilcox.runtime, 
                         glm.runtime))
performance
```
```{r echo = FALSE}
plot_style <- function(plot){
  return(
    plot + 
      theme_classic() + 
      theme( axis.text.x = element_text(angle = 45, hjust = 1), 
             legend.position = "None",
             plot.title = element_text(hjust = 0.5)) +
      xlab("")+
      scale_y_continuous(expand = expansion(mult = c(0, 0.05)))
  )
}


Sensitivity.plot <- ggplot(performance, aes(x = reorder(Package, -Sensitivity), 
                                            y = Sensitivity, 
                                            fill = Package))+
                          geom_bar(stat= "identity", color = "black")+
                          ggtitle("Sensitivity")
Sensitivity.plot <- plot_style(Sensitivity.plot)  

Specificity.plot <- ggplot(performance, aes(x = reorder(Package, -Specificity), 
                                            y = Specificity, 
                                            fill = Package))+
                          geom_bar(stat= "identity", color = "black")+
                          ggtitle("Specificity")
Specificity.plot <- plot_style(Specificity.plot)  

F1.plot <- ggplot(performance, aes(x = reorder(Package, -F1), 
                                            y = F1, 
                                            fill = Package))+
                geom_bar(stat= "identity", color = "black")+
                ggtitle("F1") 
F1.plot <- plot_style(F1.plot)  

Runtime.plot <- ggplot(performance, aes(x = reorder(Package, Runtime), 
                                   y = Runtime, 
                                   fill = Package))+
  geom_bar(stat= "identity", color = "black")+
  ggtitle("Runtime") + 
  ylab("Runtime (s)")
Runtime.plot <- plot_style(Runtime.plot)  

ggarrange(Sensitivity.plot, Specificity.plot, F1.plot, Runtime.plot, ncol = 4)
```
   
---


As is seen above, DEGage provides one of the best trade offs between performance and run time. Note that this is not meant to be a representative comparative analysis of these methods, and is purely for illustrative purposes. For a more robust comparison, see [(1)](#ref1)

# References

---

[(1)]{#ref1} Petrany A., Chen R., Zhang S. and Chen, Y. "Theoretical Framework for the Difference of Two Negative Binomial Distributions and its Application in Comparative Analysis of Sequencing Data". Genome Research, Under Review.   
[(2)]{#ref2} Love, M.I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014). https://doi.org/10.1186/s13059-014-0550-8   
[(3)]{#ref3} Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010 Jan 1;26(1):139-40. doi: 10.1093/bioinformatics/btp616. Epub 2009 Nov 11. PMID: 19910308; PMCID: PMC2796818.    
[(4)]{#ref4} Zhun Miao, Ke Deng, Xiaowo Wang, Xuegong Zhang, DEsingle for detecting three types of differential expression in single-cell RNA-seq data, Bioinformatics, Volume 34, Issue 18, September 2018, Pages 3223–3224, https://doi.org/10.1093/bioinformatics/bty332  
[(5)]{#ref5} Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn JL. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014 Apr;32(4):381-386. doi: 10.1038/nbt.2859. Epub 2014 Mar 23. PMID: 24658644; PMCID: PMC4122333.    
[(6)]{#ref6} Korthauer, K.D., Chu, LF., Newton, M.A. et al. A statistical approach for identifying differential distributions in single-cell RNA-seq experiments. Genome Biol 17, 222 (2016). https://doi.org/10.1186/s13059-016-1077-y  














