How well do RNA-Seq differential gene expression tools perform in a eukaryote with a complex transcriptome?

RNA-seq experiments are usually carried out in three or fewer replicates. In order to work well with so few samples, Differential Gene Expression (DGE) tools typically assume the form of the underlying distribution of gene expression. A recent highly replicated study revealed that RNA-seq gene expression measurements in yeast are best represented as being drawn from an underlying negative binomial distribution. In this paper, the statistical properties of gene expression in the higher eukaryote Arabidopsis thaliana are shown to be essentially identical to those from yeast despite the large increase in the size and complexity of the transcriptome: Gene expression measurements from this model plant species are consistent with being drawn from an underlying negative binomial or log-normal distribution and the false positive rate performance of nine widely used DGE tools is not strongly affected by the additional size and complexity of the A. thaliana transcriptome. For RNA-seq data, we therefore recommend the use of DGE tools that are based on the negative binomial distribution.


Introduction
Short read RNA sequencing (RNA-seq) has become the method of choice for transcriptomewide quantification of gene expression and the analysis of differential gene expression between experimental conditions (Mortazavi et al. 2008;Nagalakshmi et al. 2010). RNA-seq data analysis typically involves aligning short sequence fragments (reads) to a reference genome or transcriptome, counting the resulting alignments that fall within an annotated feature region, then identifying any significant differences between two conditions. More than a dozen computational tools have been developed to identify Differential Expression (DE) from RNA-seq data and each makes assumptions about the variability of the RNA-seq expression measurements across replicates (Anders and Huber 2010;Hardcastle and Kelly 2010;Robinson et al. 2010;Wang et al. 2010;Tarazona et al. 2011;Li et al. 2012;Lund et al. 2012;Trapnell et al. 2012;Li and Tibshirani 2013;Leng et al. 2013;Frazee et al. 2014;Law et al. 2014;Love et al. 2014;Moulos and Hatzis 2014;Ritchie et al. 2015). Based on these assumptions, the tools calculate the probability that two sets of measurements come from the same statistical distribution, thus determining whether a genuine shift in expression is a more likely explanation for the observed values than random chance. Assuming an incorrect distribution can lead to poor False Discovery Rate (FDR) control and inaccurate True Positive (TP) identification in the DGE calls. Such errors will propagate downstream into the biological interpretation of the DE results. Most commonly, these tools are used to identify DE for genes (Differential Gene Expression -DGE), however they are increasingly being used to identify DE of other genic regions (i.e. exons, spliced transcripts, etc.) (Wood et al. 2013;Frazee et al. 2014;Gaidatzis et al. 2015).
However, these studies were carried out using either simulated data or biological data that was originally designed for a different purpose. Although a few of these studies have explored high biological replication by leveraging publicly available datasets on individuals within a species (Seyednasrollah et al. 2015;Guo et al. 2013;Soneson and Delorenzi 2013;Burden et al. 2014;Bottomly et al. 2011), most have a limited level of replication. Recently, a study was performed in yeast specifically to test both the underlying statistical properties of RNA-seq data across biological and technical replicates and the influence of replication on DGE results (Gierliński et al. 2015;Schurch et al. 2016). With 48 biological replicates per condition, the distribution of read counts per gene was found to be most consistent with a negative binomial (NB) distribution, and the FDR of many DGE tools was demonstrated to be low at all replication levels. However, the number of replicates was shown to affect the tools' sensitivity for identifying differential expression, particularly at low fold-changes. The authors recommended that future RNA-seq experiments that focus on identifying DE have at least six replicates per condition in order to reliably detect genes with fold-changes greater than a factor of 2, rising to at least twelve replicates to reliably detect the majority of differentially expressed genes regardless of fold-change.
In this paper, RNA-seq data from 17 wild-type biological replicates of Arabidopsis thaliana were used to explore read count measurements across replicates and the FDR of DGE tools.
Although A. thaliana has a relatively small genome, its transcriptome is similar in scale and complexity to that of model mammal species (Arabidopsis Genome Initiative 2000;Carvalho et al. 2013;Krishnakumar et al. 2014) and its genome is extensively annotated. Accordingly, conclusions from the high replicate RNA-seq study presented here should provide useful guidance for work in other complex eukaryotes as well.  Table 1) and high level of multi-mapping reads (Supplementary   Table 2). Consequently, replicate 11 was excluded from the subsequent analyses.

Consistency among replicates
Additionally, a low but uniform level of read coverage across the genome was observed in the replicates from ExpB (replicates 8-14), explaining the marginally lower correlation between the replicates of this experiment and the other replicates ( Figure 1, right panel). Figure 2 shows that the distribution of expression measurements across replicates is consistent with being drawn from a negative binomial distribution for 100% of the genes, and is consistent with a log-normal distribution for 98% of genes. In contrast, 23% of genes have expression distributions that are not consistent with a normal distribution and more than 70% are inconsistent with a Poisson distribution.

Distribution of gene read counts across replicates
As mentioned previously, replicates 8-14 presented a low level of uniform read coverage along the genome. We believe this to be noise consistent with a small amount of DNA contamination in the affected samples. Such reads have the potential to interfere with the fitting of statistical distributions to the data, as they can make silent genes artificially appear as expressed. Indeed, Figure 3 shows that approximately 6,000 of the genes annotated in Araport11 appear to be lowly expressed in the affected replicates, but are not detected in the ten replicates from the other two experiments (replicates 1-7 and 15-17). The potential for this noise to impact on the distribution measurements was assessed by comparing the fit of the distributions to the full dataset of 16 replicates against the the fit without the noisy replicates of ExpB. As a control, the distribution fitting was also performed with the exclusion of replicates 1-7 in place of the noisy ones. In both cases, the fit of the distributions to the reduced dataset was very similar to that of the full dataset, with the exception of the normal distribution where, in both the noise exclusion and the control, the fraction of the gene data that is inconsistent with the distribution was reduced from 23% to ~10%. The lack of difference between excluding the noisy replicates and excluding the control replicates demonstrates that the apparent improvement of model fit is due to the reduction of statistical power caused by the smaller number of replicates, rather than to a cleaner signal.  performance is worse than the other tools at all the tested numbers of replicates, suggesting that it is a poor choice for calling SDE.

Discussion
In this study the statistical assumptions made by tools that identify differential gene expression from RNA-seq read count data were validated in a high-replicate experiment. The findings show that the negative binomial (NB) and log-normal distributions are both good choices as models for the cross-replicate variability of RNA-seq read counts. The study demonstrates that 6 out of the 9 DGE tools examined here control their identification of false positives well even with only 3 replicates. In contrast, tools that assume distributions other than negative binomial or log-normal showed higher FPR, at all replication levels. These results reinforce similar conclusions previously reached by studies of the yeast transcriptome (Gierliński et al. 2015;Schurch et al. 2016).
The transcriptome of A. thaliana is considerably more complex than S. cerevisiae, with almost four times the number of protein-coding genes (27,667 in A. thaliana, 7,126 in S. cerevisiae) and widespread alternative splicing and alternative polyadenylation. The similarity of the results from these two diverse organisms suggests that the conclusions of both studies are applicable to a wide range of eukaryotes.
However, the concept of gene expression in complex transcriptomes is confounded by the presence of alternative transcript isoforms, which give the organism additional means to regulate a gene's expression. This type of regulation is not necessarily reflected in changes to the total transcriptional output by a gene. Ideally, expression studies should aim to quantify the abundance of alternative isoforms individually and independently. Interestingly, the sum of independent random variables with a NB distribution itself has a NB distribution. Thus, the conclusion that a NB is a suitable model for gene expression variability across replicates is consistent with the hypothesis that the underlying variability of expression of the individual isoforms also follows the NB distribution. This suggests that tools originally intended for the study of differential gene expression may also be appropriate for studying differential transcript expression (DTE).
In summary, our analyses show that the statistical properties of gene expression are similar between simple and complex eukaryotes, and validate the model assumptions of the bestperforming DGE tools, thus justifying their use in a broad range of organisms. thaliana genome using the splice-aware aligner STAR v2.5.0a (Dobin et al. 2012). The index was built with "--sjdbOverhang 99" and the alignment was run with parameters:"--
These read counts were used without further processing to examine the FPR performance of nine DGE tools, allowing each tool to carry out its default normalization. The tools were used in the R v3.2.2 environment and installed through Bioconductor v3.2.
For the purposes of comparing the expression distribution models, consistently normalized data was required. As some of the distributions in question are discrete, normalized integer read counts were used for this purpose, which were calculated by randomly down-sampling read-pairs from each replicate to the level of the replicate with the lowest read depth. In the present study, the focus is on the collective behavior of gene expression, rather than the biological interpretation of the expression of any specific gene, so this type of normalization is appropriate here. However, it is not recommended for typical gene expression analysis studies, as some low expression signals can randomly be lost during resampling. After the normalization, each replicate consisted of ~ 7 7 ൈ 1 0 read pairs, which were then aligned to the genome and quantified as described above.
14 Performing the tests The read counts of each gene were tested against four theoretical distributions across replicates: normal, log-normal, Poisson and negative binomial (NB). For the normal and lognormal distributions the goodness of fit was determined using the test for normality from D 'Agostino et. al. (1990). For the NB distribution, the method described by Meintanis (2005) was employed and for the Poisson distribution a χ 2 test (Fisher 1950). In each case, rejection of the null hypothesis was based on a Benjamini-Hochberg corrected critical p-value of 0.05 (Benjamini and Hochberg 1995). In order to test the false positive rate (FPR) of the DGE tools, two sets of ݊ replicates were randomly selected without replacement from the pool of 16 WT replicates. Differential gene expression was then called on each of the set pairs with each of nine DGE tools (Table 1), allowing each tool to apply its default normalization. Since all sets are drawn from the same WT pool, every gene identified as significantly differentially expressed (SDE) is, by definition, a false positive. This process was repeated 100 times for each sample size in the range 2 ݊ 7 for each tool.

Availability
The raw data for the 17 WT Arabidopsis thaliana datasets is available from the European Nucleotide Archive ([database accession to be inserted here]).