Getting the most out of RNA-seq data analysis

Criteria for differential expression

For edgeR, DESeq, DESeq2 and Z-test, we used a joint filtering criteria (Li, 2012) based on fold change (ϕ) and p-value (p) to call DEG. Let y = − log₁₀p and x = log₂ϕ. Thus, each gene is associated with a paired score (x, y) after differential expression analysis. Following (Xiao et al., 2014), we required p < 0.01 and ϕ ≥ 2 to call for up-regulated genes, and p < 0.01 and ϕ ≤ 1/2 to call for down-regulated genes. The product of y > 2 and |x| ≥ 1 yields the inequality y > 2/|x|. Thus, genes that fell in the region defined by y > 2/x were differentially up-regulated, and those in the region of y > − 2/x were differentially down-regulated. The union of the sets of differentially up and down-regulated genes constituted the set of DEG candidates.

To handle analysis of unreplicated experiments in edgeR, we set the biological coefficient of variation (BCV) parameter as 0.4 for the Cheung data set (see details in ‘Benchmarking’), and 0.1 for the Bottomly data set, following recommendations in Chen et al. (2015). The exact test option was used to compute p-values.

For NOISeq, we used the recommended criteria for calling DEG as described in the NOISeq documentation—q = 0.9 for unreplicated experiments, and q = 0.95 for experiments with biological replicates. For ASC, genes that had log2 fold change above 1 or −1, and posterior probability 99% or more were declared to be differentially expressed. For GFOLD, we used the default significant cut-off of 0.01. A gene with GFOLD value of 1 or larger was considered differentially up-regulated, and differentially down-regulated if GFOLD value was −1 or smaller. Except GFOLD, which is written in the C/C+ + language and requires the Linux platform, the other methods were executed in R version 3.1.3 (R Core Team, 2015).

Data sets

To set up our benchmarking exercise, we needed two RNA-seq data sets whereby variation in their phenotype classes produced mild and strong biological effect sizes in the tissue of interest, respectively. We further required the RNA-seq data sets to have fairly large replicate sizes to enable the simulation of different replicate size scenarios. To this end, we identified two suitable RNA-seq data sets in the Recount database (Frazee, Langmead & Leek, 2011). The latter contains unnormalized RNA-seq count data sets from 18 major studies that have been assembled from raw reads using the Myrna (Langmead, Hansen & Leek, 2010) pipeline.

The Bottomly data set (Bottomly et al., 2011) consists of gene expression data (22 million Illumina reads per sample, read length of ∼30 bases) obtained from the brain striatum tissues of two mice strains: C57BL/6J (n = 10) and DBA/2J (n = 11). Both mice strains are known to show large, strain-specific variation in neurological response when subjected to opiate drug treatment (Korostynski et al., 2006; Korostynski et al., 2007; Grice et al., 2007).

The Cheung data set (Cheung et al., 2010) consists of gene expression data (40 million Illumina reads per sample, read length of 50 bases) from immortalized human B-cells of 24 males and 17 females. Sex hormones are known to modulate B cell function (Klein, 2000; Verthelyi, 2001). For example, estrogen modulates B cell apoptosis and activation (Grimaldi et al., 2002), while testosterone suppresses immunoglobulin production by B cells (Kanda, Tsuchida & Tamaki, 1996). In the absence of antigenic challenge, however, it seems reasonable to expect only a modest number of DEG in male and female B cells.

After removal of transcripts with zero counts in all samples, the Bottomly count table contained 13,932 transcripts, down from an initial 36,536 transcripts, whereas the Cheung count table contained 12,410 transcripts, down from 52,580. Prior to analysis, the count data were normalized using DESeq normalization (Anders & Huber, 2010), which has been shown to be robust to library size and composition variation (Dillies et al., 2013). However, raw counts were used for DESeq2 analysis since the method explicitly requires such type of data as input.

Method for constructing a reliable reference DEG set

The construction of a reliable reference DEG set from which performance metrics for each method is evaluated is a non-trivial problem, if one eschews a simulation-based approach. To avoid circular reasoning, this reference set needs external validation from independent evidence such as confirmation from qPCR results.

Here, we chose voom (Law et al., 2014; Ritchie et al., 2015) as the method of choice for setting the reference DEG set. Unlike other DEG methods that primarily model mean–variance relationships in the count data using discrete distributions such as the Poisson or negative binomial distributions, voom log-transforms count data into a microarray-like data type suitable for analysis using the robust limma pipeline (Smyth, 2004; Ritchie et al., 2015). Because of this, using voom to set the reference DEG set can avoid biasing results of the called DEG due to algorithmic similarities. A gene was defined as differentially expressed using the same joint filtering criteria for edgeR, DESeq, DESeq2 and Z-test. We found the nonparametric SAMSeq (Li & Tibshirani, 2013), which has also been reported to have strong DEG mining performance, unsuitable for setting the reference DEG set as it returned different DEG sets for different random seeds and number of permutation parameters (Fig. S1).

The validity of voom as a tool for constructing reasonable in silico reference DEG sets for the Bottomly and Cheung data set requires justification. To this end, we compared its performance with other DEG call methods on an RNA-seq data set in which qPCR validation results for sufficiently large numbers of genes are available. Briefly, the Rajkumar data set (Rajkumar et al., 2015) consists of gene expression count data (26,119 genes; minimum of 10 million Illumina reads per sample, read length of ∼50 bases) from the amygdala tissues of C57BL/6NTac strain mice. There are two phenotype classes: wild type (n = 8), and heterozygotes for the Brd1 gene deletion (n = 8). A total of 115 genes were selected for qPCR validation (additional Table 5 in Rajkumar et al., 2015); differential expression was observed in 60 of them, and not in the remaining 55. Each DEG call method returns N_g differentially expressed gene candidates. We considered a method to be sound for setting the reference DEG set if it did not return too few (tens) or too many (thousands) candidates. Among methods that satisfied this criterion, the one that had relatively higher positive predictive value (PPV; the complement of the false discovery rate) would be preferred. Let N_TP be the number of true positives, and N_FP the number of false positives. Then the number of DEG that lack validation result is U = N_g − N_FP − N_TP. If U is not too large or too small, then the expected number of true positives can be estimated as ${\displaystyle N_{TP}^{\ast }=N_{TP}+\frac{N_{TP}}{N_{TP}+N_{FP}}U.}$ The expected PPV is therefore given by ${\displaystyle PP{V}^{\ast }=\frac{N_{TP}+N_{TP}^{\ast }}{N_g}.}$

Characteristics of constructed reference DEG sets

Ideally, the in silico reference DEG set called using voom for the two test data sets should be independently validated using qPCR, but evidence at such level may not always be available. Where microarray data are available for the same study, a DEG candidate can be considered reliable if it is called in both RNA-seq and microarray analyses, since fold change of DEG from the latter has been found to correlate strongly with fold change from qPCR (Wang et al., 2014). A total of 362 DEG for the Bottomly data set were thus called (Fig. 1A). About 88% (320/362) of the DEG for the Bottomly data called using voom were identical with those called in Bottomly et al. (2011) using edgeR (1,727 DEG). Approximately two fifths of them (153/362) were detected using limma applied on Affymetrix or Illumina microarray expression data (Table S1).

For the Cheung data set, gender difference was the source of phenotype class variation. We exploited this fundamental biological difference to infer the most reliable DEG from the candidates returned using voom. Only DEG which were located on the sex chromosomes, or interacted with at least one gene product from the sex chromosomes were used to construct the reference DEG set. This strategy resulted in a set of 19 DEG (Fig. 1B). Five of them were located on the Y chromosome, three on the X chromosome and the remainder had known gene-gene interactions with at least one gene located on sex chromosomes (based on BioGRID (Stark et al., 2006; Chatr-Aryamontri et al., 2015); Table S2).

Differentially expressed genes are characterized by between-phenotype variation that is significantly larger than within-phenotype variation. However, occasionally some genes may be wrongly declared as differentially expressed because some outliers within a phenotype class were sufficiently extreme to cause relatively large between-phenotype variation. To assess the quality of DEG called using voom, we used the Bland-Altman plot (Bland & Altman, 1986; Fig. S2). Among the 362 DEG called for the Bottomly data set, the majority of DEG showed good agreement of replicate variation between the two phenotype classes—about 88% (318/362) were within 2SD (standard deviation) from perfect agreement, and about 95% (345/362) were within 3SD. Similarly, among the 19 DEG called for the Cheung data set, within-phenotype variation difference was within 2SD from perfect agreement for about 74% (14/19) of DEG, and within 3SD for about 84% (16/19) of DEG. Generally, genes that showed large within-phenotype variation in both phenotypes were not called by voom.

Once the DEG set had been constructed for the Bottomly and Cheung sets, it became possible to operationally define what we meant by mild or strong biological effect size. For the ith differentially expressed gene, define ${\displaystyle T_i^2=\frac{{\left({\bar{X}}_{i,1}-{\bar{X}}_{i,2}\right)}^2}{S_{i,1}^2/n_1+S_{i,2}^2/n_2},}$ where i indexes the genes, and j the phenotype classes (j = 1, 2); ${\bar{X}}_{i,j}$ and $S_{i,j}^2$ are the mean and variance of normalized read counts respectively, and n_j are the replicate sizes. Thus, T² is essentially the square of the t-statistic, which measures the magnitude of squared deviation between mean counts in two different phenotype classes relative to the latters’ variances. By definition, the median values of T² should be large in a data set that shows strong biological effect size, and vice versa. For the Bottomly data set (strong effect size), median T² was 27.6; for the Cheung data set (mild effect size), it was 4.6. Both data sets had approximately equal spread of T² values around the median, the interquartile range being 38.3 and 34.5 for the Bottomly and Cheung data sets, respectively.

Positive predictive value and sensitivity

The ASC package provided by Wu et al. (2010) failed to run for particular combinations of sample pairs. Only 24.5% (27/110) and 41.3% (124/300) pairs of samples from the Bottomly and Cheung data set respectively could analyzed using ASC. The simulation results show that optimality of a DEG call method for a given replicate size depended on whether biological effect size was mild or strong (Fig. 2). In the Cheung data set (mild biological effect size), all methods had very low (about 1%) mean positive predictive value (PPV) for unreplicated experiment (Fig. 2A and Table S3), suggesting that no meaningful biological insights were possible. However, mean PPV (±SD) increased substantially for NOISeq to 43.5 ± 31.5%, and for GFOLD to 29.6 ± 15.8%, for n = 3 (Fig. 2B). Doubling and approximately tripling replicate size to n = 6 (Fig. 2C) and n = 10 (Fig. S3) further improved mean PPV for NOISeq to 87.0 ± 16.1% and 92.2 ± 12.9%, and for GFOLD to 36.3 ± 14.9% and 52.6 ± 18.8%, respectively. In all four replicate sizes, mean PPV was low for the other methods. It did not exceed 12% for DESeq2, and was never more than 3% for edgeR, DESeq and Z-test.

A markedly different pattern of method performance was observed in the analysis of the Bottomly data set (strong biological effect size). In unreplicated experiments (Fig. 2D), mean PPV was relatively high for NOISeq (50.6 ± 20.3%), ASC (47.2 ± 25.9%) and GFOLD (31.2 ± 25.6%), compared to just about 15% in edgeR and 5% in DESeq and Z-test. DESeq2 did not perform well, with mean PPV (29.6 ± 28.4%), and an extremely low sensitivity (0.2 ± 0.6%) as a result of making too few calls. Interestingly, GFOLD attained very high mean PPV at n = 3 (94.3 ± 6.9%; Fig. 2E), with marginal change to 92.5 ± 3.3% at n = 6. However, GFOLD was also the method with the lowest sensitivity (below 10%) for these two replicate sizes, which was caused by its small DEG set size (Fig. 3).

DESeq2 struck the best balance between PPV and sensitivity as replicate size increased, but edgeR showed reasonable performance too. At n = 3 (Fig. 2E) and n = 6 (Fig. 2F), DESeq2 had mean PPV of 52.5 ± 10.8% and 62.1 ± 7.7%, with mean sensitivity of 36.0 ± 5.7% and 65.1 ± 4.5%, respectively. For edgeR, its mean PPV was 28.7 ± 4.1% and 33.9 ± 3.0%, with mean sensitivity of 59.8 ± 5.4% and 79.0 ± 4.6%, respectively. At n = 6, DESeq2 had similar sensitivity compared to its older version DESeq, and a superior mean PPV that was about four times higher. Unsurprisingly, the Z-test remained the worst performer, with mean PPV just about 6%.

The general increase in mean sensitivity for replicated experiments was consistent with the finding that the increase in statistical power for detecting DEG is primarily determined by biological replicate size, and less by sequencing depth (Liu, Zhou & White, 2014).

DEG set size

Figure 3 shows the distribution of DEG set size in the Cheung and Bottomly data sets for different replicate sizes (for details, see Table S3). Although DESeq2 could be used to call DEG for unreplicated experiments, it was shown to behave erratically in the Bottomly data set, with extremely low mean DEG set size (2.5 ± 6.8). In general, for replicated studies, methods such as DESeq2, DESeq, edgeR and Z-test made large numbers of calls that were typically one or two order of magnitudes more (depending on underlying biological effect size) compared to GFOLD or NOISeq. Consequently, it is expected that the formers’ sensitivity would increase at the expense of their PPV.