Comparison of Bioconductor Packages edgeR and DESeq for Analyzing Differential Expression

Next-generation sequencing (NGS) technologies have greatly improved the landscape of genomic research and have, within the past decade, served to reduce dramatically the cost and increase the data output potential.  One of the specific applications that utilizes NGS methods is RNA-sequencing (RNA-Seq), which has a number of uses, one of which is the determination of differential (gene) expression (DE) among two or more RNA samples.  Although the cost of NGS technology continues to decrease and open the field of genomic research to a broader scope of scientists, the cost of large-scale experiments in which sufficient biological and technical replicates are used and for which standard, large-sample statistical methods are appropriate are still cost prohibitive to all but the most well-funded industry or government enterprises.  Thus, academic scientists, while prolifically conducting genomic research on DE, rely on experimental designs that include few biological and/or technical replicates.  Additionally, several other experimental conditions, for example, differences among sample library sizes and the conduct of multiple hypothesis tests, confound the analysis.  As a consequence, new analytic tools were needed to account for these adjusted model assumptions.  R/Bioconductor supports a multitude of statistical analysis packages designed for the analysis of DE.  These packages include edgeR, DESeq, DESeq2, NBPSeq, baySeq, and ShrinkSeq, among others—each of which differs in one or more aspects, including the way it handles normalization, parameter estimation, hypothesis testing, and correction for multiple testing.  Whereas different analysis methods are not at all unusual, a complicating factor for the analysis of DE is the fact that, given the same assumptions, the analysis packages produce different results—often significantly.  This project explores the analytic differences and performance among several R/Bioconductor packages, including edgeR and DESeq.  To conduct this analysis, a count dataset was simulated using model parameters estimated from a real RNA-seq dataset.  Count data for the treatment condition were then simulated to reflect DE (i.e., either up- or down-regulated) for a pre-specified number of genes.  The simulated dataset—for which DE status was known—was then challenged using each of the analysis packages.  The results and performance, including Type I and Type II error rates, were compared across packages.