Comparison of software packages for detecting differential expression in RNA-seq studies

Abstract

RNA-sequencing (RNA-seq) has rapidly become a popular tool to characterize transcriptomes. A fundamental research problem in many RNA-seq studies is the identification of reliable molecular markers that show differential expression between distinct sample groups. Together with the growing popularity of RNA-seq, a number of data analysis methods and pipelines have already been developed for this task. Currently, however, there is no clear consensus about the best practices yet, which makes the choice of an appropriate method a daunting task especially for a basic user without a strong statistical or computational background. To assist the choice, we perform here a systematic comparison of eight widely used software packages and pipelines for detecting differential expression between sample groups in a practical research setting and provide general guidelines for choosing a robust pipeline. In general, our results demonstrate how the data analysis tool utilized can markedly affect the outcome of the data analysis, highlighting the importance of this choice.

Keywords: RNA-seq; differential expression; gene expression.

Figure 1:

Effect of normalization on the detections. (A) Overlaps of the differentially expressed genes detected in the mouse and human data using the default or TMM normalization method for all R-based packages that do not use TMM as their default normalization option (see Methods for details of the significance thresholds). With edgeR and limma the TMM normalization is the default normalization and, therefore, they are not included. Overall similarity between the rankings obtained using the default or TMM normalization method (denoted by TMM after the name of the package) in the (B) mouse and (C) human data. The dendrogram was constructed using average linkage hierarchical clustering and Spearman correlation of the gene ranks.

Figure 2:

Number and consistency of differentially expressed genes detected using eight state-of-the-art software packages in the mouse and human data (upper and lower panel, respectively). (A) Number of detections (y axis) with different numbers of replicates (x axis) for each software package. The points correspond to averages over 10 randomly sampled subsets; the error bars show the standard error of the mean. (B) Differentially expressed genes in the complete data divided into four categories on the basis of their expression levels: very lowly or not expressed genes, lowly, medium and highly expressed genes. The different software packages were ordered on the basis of their total number of detections in the mouse data. (C) Precision of the detections (y axis) when increasing the number of replicates (x axis) in terms of genes identified as differentially expressed genes in the complete data using all the samples available in the mouse and human data (upper and lower panel, respectively). Only statistically significant genes were considered with each method (see Methods for details of the significance thresholds). The points correspond to averages over 10 randomly sampled subsets; the error bars show the standard error of the mean. A colour version of this figure is available at BIB online: http://bib.oxfordjournals.org.

Figure 3:

False discoveries on the basis of mock comparisons in the (A) mouse and (B) human data. In each mock comparison, differentially expressed genes were identified between two artificially constructed sample subsets from a single sample group, in which no significant detections are expected. To compare between the different software packages, we divided the number of mock detections with the average number of detections in the actual comparisons with the same number of replicates. Only statistically significant genes were considered with each method (see Methods for details of the significance thresholds). The points correspond to averages over 10 randomly sampled subsets; the error bars show the standard error of the mean. A colour version of this figure is available at BIB online: http://bib.oxfordjournals.org.

Figure 4:

Similarity between the methods in the mouse and human data (upper and lower panel, respectively). (A) Overall similarity between the methods based on Spearman correlation of gene ranks. The dendrograms were constructed using average linkage hierarchical clustering. (B) Overlap of significant detections between the methods (see Methods for details of the significance thresholds). The proportion of common detections was calculated for each pair of methods, resulting in an asymmetric matrix of percentages. A colour version of this figure is available at BIB online: http://bib.oxfordjournals.org.

Figure 5:

Runtimes of the different methods to identify differentially expressed genes in the human data. Time in seconds on log scale (y axis) is shown as a function of the number of replicates (x axis). The analyses were run on a computer cluster node with two Intel XEON Hexa-Core processors and 96 GB of memory. Cuffdiff 2 was run on all the available 12 cores; the R/Bioconductor packages were run on a single core. A colour version of this figure is available at BIB online: http://bib.oxfordjournals.org.