BACKGROUND
rna sequencing  allows an entire transcriptome to be surveyed at single-base resolution whilst concurrently profiling gene expression levels on a genome scale
<cit> .
rna-seq is an attractive approach as it profiles the transcriptome directly through sequencing and therefore does not require prior knowledge of the transcriptome under consideration.
an example of the use of rna-seq as a high-resolution exploratory tool is the discovery of thousands of additional novel coding and noncoding genes, transcripts and isoforms of known genes despite the prior extensive annotation of the mouse
<cit>  and human genomes
<cit> .
arguably, the most popular use of rna-seq is profiling of gene expression or transcript abundance between samples or differential expression .
the efficiency, resolution and cost advantages of using rna-seq as a tool for profiling de has prompted many biologists to abandon microarrays in favour of rna-seq
<cit> .
despite the advantages of using rna-seq for de analysis, there are several sources of sequencing bias and systematic noise that need to be considered when using this approach.
clearly, rna-seq analysis is vulnerable to the general biases and errors inherent in the next-generation sequencing  technology upon which it is based.
these errors and biases include: sequencing errors , biases in sequence quality, nucleotide composition and error rates relative to the base position in the read
<cit> , variability in sequence depth across the transcriptome due to preferential sites of fragmentation, variable primer and transcript nucleotide composition effects
<cit>  and finally, differences in the coverage and composition of raw sequence data generated from technical and biological replicate samples
<cit> .
recently, there have been several investigations
<cit>  into the biases that affect the accuracy with which rna-seq represents the absolute abundance of a given transcript as measured by high precision approaches such as taqman rt-pcr
<cit> .
it has been shown that these abundance measures are prone to biases correlated with the nucleotide composition
<cit>  and length of the transcript
<cit> .
several within and between sample correction and normalisation procedures have recently been developed to address these biases either as nucleotide composition effects
<cit>  or various combinations of nucleotide, length or library preparation biases
<cit> .
these approaches all yield improvements in the correspondence of rna-seq read counts with expression estimates gained by other experimental approaches.
despite the known biases, rna-seq continues to be widely and successfully used to profile relative transcript abundances across samples to identify differentially expressed transcripts
<cit> .
the profile of a given transcript across a biological population would be hoped to be less prone to nucleotide composition and length biases as these variables remain constant.
nevertheless, to accurately detect de across samples it is necessary to understand the sources of variation across technical and biological replication and where possible respond to these with an appropriate experimental design and statistically robust analysis
<cit> .
to date, there has been little discussion in the literature of efficient experimental designs for the detection of de and a lack of consensus about a standard and comprehensive approach to counter the many sources of noise and biases present in rna-seq has meant that some of the biological community remain sceptical about its reliability and unsure of how to design cost-efficient rna-seq experiments .
good experimental design and appropriate analysis is integral to maximising the power of any ngs study.
with regard to rna-seq, important experimental design decisions include the choice of sequencing depth and number of technical and/or biological replicates to use.
for researchers with a fixed budget, often a critical design question is whether to increase the sequencing depth at the cost of reduced sample numbers or to increase the sample size with limited sequencing depth for each sample
<cit> .
sequencing depth
sequencing depth is usually referenced to be the expected mean coverage at all loci over the target sequence, in the case of rna-seq experiments assuming all transcripts having similar levels of expression.
without the benefit of extensive previous rna-seq studies, it is difficult in most cases to estimate prior to data generation the optimal sequencing depth or amount of sequencing data required to adequately power the detection of de in the transcriptome of interest.
pragmatically, rna-seq sequencing depth is typically chosen based on an estimation of total transcriptome length  and the expected dynamic range of transcript abundances.
given the dynamic nature of the transcriptome, the suitability of these estimates could vary substantially across organisms, tissues, time points and biological contexts.
wang et al.
<cit>  found a significant increase in correlation between gene transcripts observed and number of sequence reads generated when increasing sequencing depth from  <dig>  to  <dig> million reads after which the gains plateau –  <dig> million reads detected about 80% of the annotated chicken transcripts.
despite the expectation of continuous sequencing depth increases in the near future, Łabaj et al.
<cit>  argue that most of the additional reads will align to the subset of already extensively sampled transcripts.
as a result, transcripts with low to moderate expression levels will remain difficult to quantify with good precision using current rna-seq protocols even at higher read depths.
greater sequencing depth will also increase sensitivity to detect smaller changes in relative expression, however this does not guarantee that these changes have functional impact in the biological system under study as opposed to tolerated fluctuations in transcript abundance
<cit> .
ideally, an efficient experimental design will be informed by an understanding of when increasing sequencing depth begins to provide rapidly diminishing returns with regard to transcript detection and de testing.
replication
replication is vital for robust statistical inference of de.
in the context of rna sequencing, multiple nested levels of technical replication exist depending upon whether it is the sequence data generation, library preparation or rna extraction technical processes that are being replicated from the same biological sample.
several published studies have incorporated technical replicates into their rna-seq experimental designs
<cit> .
the degree of technical variation present in these datasets appears to vary and the main source of technical variation appears to be library preparation
<cit> .
biological replication measures variation within the target population and simultaneously can counteract random technical variation as part of independent sample preparation
<cit> .
it has been shown that power to detect de improves when the number of biological replicates n is increased from n =  <dig> to n = 5
<cit> , however, to date few studies have incorporated extensive biological replication and extensive testing of the effects of replication on power is needed.
more recently with the increasing utility and availability of multiplex experimental designs, the incorporation of biological replicates with decreased sequencing depth is becoming a much more attractive and cost-effective strategy.
the relative merits of sacrificing sequencing depth for increased replication has not been rigorously explored.
efficient experimental design
multiplexing is an increasingly popular approach that allows the sequencing of multiple samples in a single sequencing lane or reaction and consequently the reduction in sequencing costs per sample
<cit> .
multiplexing uses indexing tags, “barcodes” or short  stretches of sequence that are ligated to the start of sample sequence fragments during the library preparation step.
barcodes are distinct between sample libraries and allow pooling for sequencing followed by allocation of reads back to individual samples after sequencing by analysis of the sequenced barcode.
multiplex barcode designs are routinely available with up to  <dig> samples in the same lane, recently up to  <dig> yeast dna samples were profiled in single lane
<cit> .
novel methods are continuing to emerge for low-cost strategies to multiplex rna-seq samples
<cit> .
with the dramatic increases in sequencing yields being achieved with current chemistries and new platforms, multiplexing is becoming the method of choice to increase sample throughput.
these designs have direct impacts on sequencing depth generated that need to be considered in the power of the experimental design.
also, when multiplex strategies are used, biologists need to be mindful of potential systematic variations between sequencing lanes.
these variations can be addressed through randomisation or blocking designs to distribute samples across lanes, see
<cit>  for a discussion of barcoding bias in multiplex sequencing, and
<cit>  for an alternative to barcoding.
in a comparison between microarray and ngs technologies in synthetic pools of small rna, willenbrock et al.
<cit>  found that multiplexing resulted in decreased sensitivity due to a reduction of sequencing depth and a loss of reproducibility; however the authors did not investigate power for detection of de in their study.
approach
improving detection of de requires not only an appropriate experimental design but also a suitably powered analysis approach.
several algorithms have recently been developed specifically to appropriately handle expected technical and biological variation arising from rna-seq experiments.
a non-exhaustive list of these algorithms is: edger
<cit> , deseq
<cit> , nbpseq
<cit> , bbseq
<cit> , fdm
<cit> , rsem
<cit> , noiseq
<cit> , myrna
<cit> , cuffdiff
<cit> .
a thorough comparison of these packages’ performance with datasets of different properties falls beyond the scope of this study, however before considering issues relating to power and experimental design, it is important to investigate whether packages for de analysis give the correct type i error rate under the null hypothesis of no de.
to do this evaluation we considered three popular packages for de analysis of rna-sequencing data.
these packages are based on a negative binomial distribution model of read counts
<cit>  and include edger
<cit> , deseq
<cit>  and nbpseq
<cit> .
to quantify the effects of different sequencing depths and replication choices we compared a range of realistic experimental designs for their ability to robustly detect de.
using simulated data with known de transcripts allowed us to estimate the false positive rate  and true positive rate  of de calls.
the changes of these rates were used to compare the detection power yielded by each choice of number of biological replicates and sequencing depth.
in the methods section, we outline the definitions used for fpr and tpr as well as explaining the method used for the construction of the synthetic data; which includes induced differential expression, simulates the variations that biological replicates introduce and simulates loss of sequencing depth.
in our study, we test a wide range of real-world experimental design scenarios for performance under the null hypothesis and in the presence of de.
in these scenarios both the numbers of biological replicates n and the sequencing depth are varied.
this provides a comprehensive quantitative comparison of different experimental design strategies and is particularly informative for those accessing modern multiplex approaches.
RESULTS
comparisons of statistical methods: edger, deseq, and nbpseq using simulated data under the null
to test the performance of each package under the null hypothesis, we simulated sets of n “control” and n “treatment” lanes of counts in accordance with the procedure described in the methods section, for a range of values of n and with no de between treatments.
for each value of n and for each package the simulation and testing were repeated  <dig> times.
figure
<dig> shows the percentage of transcripts reported as differentially expressed at the 1% significance level by each of the three packages for a range of values of n. the height of each bar is the median value obtained from  <dig> repetitions of the synthetic data generation with its associated 90% confidence intervals.
under the null hypothesis, the percentage reported is the false positive rate  defined by eq.
<dig>  and should match the significance level of α = 1% if the package is performing correctly.
also shown are fprs for high-count transcripts  and low-count transcripts .
figure
<dig> shows an example of the p-value distribution obtained for one experiment at n =  <dig> biological replicates.
ideally, p-values should have a uniform distribution in the interval  <cit>  if the package is performing correctly.
immediately noticeable in the p-value histogram is a sharp spike in the right hand bin for low count transcripts, which is observed to be present in general for all values of n and all packages.
this is a known artifact of calculating p-values for discrete random variables using the method described in
<cit>  and summarised in our methods subsection ‘under the null hypothesis’: when count sums in both conditions are equal the computed p-value is exactly  <dig>  the situation is most likely to occur for transcripts with extremely low counts, in which case it is difficult to draw meaningful conclusions regarding de via any method.
the behaviour at the left hand end of the histogram, which drives the fprs plotted in figure
<dig>  varies considerably between packages and numbers of biological replicates.
it is affected mainly by the method used for estimating a dispersion parameter ϕi for each transcript i .
the package edger performs well for large numbers of biological replicates , for which squeezing of the dispersion estimate towards the common dispersion is minimal, and a tagwise estimate is appropriate.
for small numbers of biological replicates, because the dispersion cannot be estimated accurately on a per-transcript basis, information is borrowed from the complete set of transcripts to squeeze the estimate towards a common dispersion estimate.
for the high-count transcripts in particular, the squeezing causes the dispersion of the most highly dispersed transcripts to be underestimated, causing too many transcripts to be deemed differentially expressed, leading to an inflated fpr.
in an effort to be conservative, deseq chooses as its estimate of dispersion the maximum of a per-transcript estimate and the functional form eq.
<dig> which is fitted to the per-transcript estimates for all transcripts.
our results indicate that the method performs well for the high-count transcripts when the number of biological replicates is small , but is otherwise over-conservative.
this is generally to be preferred to an inflated fpr, as one has more evidence that what is called de is truly de.
the package nbpseq imposes the functional relationship eq.
<dig>  which appears to be too restrictive for a number of relatively highly dispersed transcripts.
for those transcripts the dispersion parameter is underestimated, leading to an overestimate of significance and hence an inflated fpr irrespective of the number of biological replicates.
based on these results we selected deseq  and edger  for use in subsequent experimental design testing.
throughout these tests, results obtained using deseq and edger are mostly compatible with each other.
however, our comparison revealed a slightly inflated fpr from edger while deseq behaves more conservatively throughout.
therefore, in the following section we will focus on the results obtained using deseq while the analogous results obtained with edger are presented in the additional file
1: figure s <dig>
comparison of statistical methods: deseq and edger using simulated data with 15% de transcripts
to test the performance of packages in the presence of an alternate hypothesis, we simulated sets of n “control” and n “treatment” lanes of counts with 15% of the transcripts either up- or down-regulated according to the procedure described in the methods section.
all results presented from this point on are derived from deseq.
detection of de as a function of number of biological replicates n
with an increase in replication we saw a steady increase in the percentage de calls by deseq , increasing from  <dig> % to  <dig> % as n increased from  <dig> to  <dig> .
as n increased, the fpr, defined by eq.
<dig> at a significance level of α = 1%, remained below  <dig> % for all values of n, and the tpr, defined by eq.
<dig> with α = 1%, increased substantially from  <dig> % to  <dig> % .
effects of biological replication on power to detect de using deseq.
fpr and tpr are defined in eqs.
<dig> &  <dig> respectively at 1%.
“call rate” is the total number of reported positives / the total number of transcripts.
these values are also represented in figure
<dig> at 100% sequencing depth.
detection of de as a function of sequencing depth
figure
<dig> represents the combined results of decreasing sequencing depth for all values of n. it can be seen that as sequencing depth decreases the tpr generated by deseq decreases monotonically across all n while the fpr remains below  <dig> %.
1n×100% sequencing depth shows the increase of tpr as more biological replicates n are used despite the loss power due to the sequencing depth reduction required by the multiplexing of lanes.
this trend remains true even for the n =  <dig> and n =  <dig> cases.
table
<dig> shows the fpr for all biological replicates n and a subset of sequencing depths: 25%, 50%, 75% and 100%, the fpr remains below  <dig> % at all sequencing depths.
table
<dig> shows the tpr reported by deseq for the same subset of sequencing depths, here the tpr increases strongly as sequencing depth increases for any number of biological replicates n.
effects of sequencing depth on fpr values for a subset of our tested depths = 25%, 50%, 75% & 100%.
effects of sequencing depth on tpr values for a subset of our tested depths = 25%, 50%, 75% & 100%.
detection of de across multiplex experimental design strategies
we simulated various scenarios of multiplexing n-control samples vs. n-treatment samples into two sequencing lanes – each control and treatment sample at a sequencing depth
=1n×100%.
in figures
<dig> and
<dig>  a solid grey line connecting every value of n at its corresponding sequencing depth provides a summary of the performance of these multiplexing scenarios.
we call this trend the “multiplex line” and it provides an insight into the results obtained by increasing the number of biological replicates used into a fixed number of sequencing lanes, in this case  <dig> sequencing lanes.
1n×100% sequencing depth.
the multiplex line in figure
<dig> shows a clear increase in tpr as replication increases despite the loss of detection power that decreasing sequencing depth induces.
it can also be seen that the fpr remains below  <dig> % for all multiplex scenarios tested .
note that for completeness we also added multiplex scenarios for n =  <dig> &n =  <dig>  whose results follow the trends well.
the multiplex line strongly favours adding more biological replicates despite the inherent loss of sequencing depth as shown by its dramatic positive slope for the tpr while maintaining a roughly constant, low fpr.
fold-changes as indicators of de
it is common practice among biologists to use fold-change, rather than p-values, as an indicator of de.
figure
<dig> shows results analogous to those of figure
<dig> when the criterion of fold-change ≥  <dig>  is used to detect de: as replication n increases, both tpr and fpr decrease because more biological replicates have the effect of averaging out differences between control and treatment lanes.
note that, as sequencing depth decreases, the fpr increases owing to the growing number of transcripts with very low numbers of counts , in which case the poisson shot noise of the sequencer can easily induce a spurious doubling or halving of counts.
this effect is ameliorated by adding  <dig> count to all transcripts prior to de analysis – doing so, does not affect the calculation of p-values .
discussion
comparisons of de algorithms: edger, deseq and nbpseq
our comparison of these three de detection algorithms under the null hypothesis revealed different performances  when different numbers of biological replicates n, are used.
deseq consistently performed more conservatively across the different n biological replicates scenarios.
deseq’s performance was closest to the expected significance level when only using high-count  transcripts while for only low-count  transcripts over-conservative behaviour is shown.
edger overestimates de detection for small values of n while its performance improves as n increases.
edger’s level of detection is constant over n when only low-count transcripts are used while overestimation increases when only high-count transcripts are used.
nbpseq overestimated detection across n for the three scenarios .
this comparison led us to use both deseq and edger throughout our replication and sequencing depth simulations.
we ultimately chose deseq’s resultsa as this package behaved slightly more conservatively and appeared less sensitive to changes in replication .
in a study by tarazona et al.
<cit> , it is argued that negative-binomial based de analysis packages like deseq and edger are highly sensitive to sequence depth increases and are therefore unable to control the fpr as sequencing depth increases.
tarazona et al.
propose a non-parametric algorithm  to calculate de based on a noise distribution created with fold-changes and the absolute differences between the transcript’s control and experiment lane counts.
however, kvam and liu
<cit>  argue that due to the small number of replicates typically used for rna-seq experiments, non-parametric methods do not offer enough detection power and suggest that current statistical methods to detect de genes based on parametric models for rna-seq data  remain a more adequate approach.
in our study we also find that both deseq and edger tend to slightly increase the fpr as sequencing depth increases – as higher depths induce deseq and edger to assign smaller p-values to transcripts with small fold-changesb.
however in no instance do we obtain a fpr larger than 1% for deseq  – .
it is worth mentioning that the latest updates to deseq  and edger  – released after the studies
<cit>  and
<cit>  reduce the number of false positive calls by about 50% .
effects of replication for detection of de
to quantify the effects of replication in rna-seq de experiments, we tested n-control vs. n-treatment biological replicates  while maintaining sequencing depth constant.
we find that as n increases both algorithms increase the call rate and tpr while the fpr remains unchanged .
our results clearly support the simple message that more biological replicates are not only desirable but needed to improve the quality and reliability of de detection using rna-seq, however, due to the costs associated with rna-seq, many experiments are likely to need to use multiplex designs to achieve this level of replication.
this study is concerned with the simulation of overdispersion effects due to biological variability and it is implied that overdispersion due to technical variability is nested within this estimation .
it is worth mentioning that, while biological variability is important, the contribution to overdispersion by technical variation is not negligible, and disagreements between estimates of expression can occur at all levels of coverage
<cit> .
ideally, rna-seq experimental design with biological replication should also aim to block sources of technical variation, such as between lane variations, to constrain the dispersion of rna-seq experiments.
effects of sequencing depth for detection of de
to quantify the effects of sequencing depth in rna-seq de experiments, we simulated an extensive sequencing depth range  for every case of n-control vs. n-treatment biological replicates.
as the amount of available sequencing data is decreased, both packages decrease the call rate and tpr while the fpr remains low.
tpr decreases very slowly as sequencing depth decreases, suggesting that sequencing depth can be reduced to ∼15% without much impact on tpr.
we conclude that de analysis with rna-seq is robust to substantial loss of sequencing data as indicated by a slow decline in tpr as sequencing depth is lost accompanied by no increase in fpr.
these findings seem consistent with the results reported by bashir et al.
<cit>  who observed that lower levels of transcriptome sequencing had sufficient information to estimate the distribution of expression values arising from observed transcripts.
bashir et al.
did not directly test power to detect de, however as testing for de relies on good concordance with the expected distribution, it follows that de is reasonably robust to loss of sequencing data.
multiplexing experimental designs
to quantify the effects of varying both n and sequencing depth, we simulated multiplexing n-control vs. n-treatment lanes into two sequencing lanes.
we observed a steady increase in tpr with the increase in n, despite the corresponding decrease in sequencing data per transcript by 1/n.
similarly, for both deseq and edger the number of de calls and the tpr increases with n, as we observed previously and is unaffected by the decrease in data.
for deseq, the fpr remains roughly constant and always below  <dig> %, while for edger, the fpr decreases slowly as n increases.
our simulations strongly support that the benefits of multiplexing n-biological replicates into one sequencing lane , far outweigh the decrease of available data per sample by 1/n.
these multiplexing experimental designs improve tpr and fpr while greatly reducing the cost of the experiment.
while the detection of de appears robust to available sequence data, there remains the question of how multiplexing affects coverage of the transcriptome and detection of low abundant or rare transcripts.
this coverage issue will increasingly be counterbalanced by rapid increases in data generation capacity from a single sequencing experiment.
in a detailed study of the marioni
<cit>  human  dataset, banshir et al.
<cit> , reported that over 90% of the total observed transcripts were sampled with  <dig> million reads.
this should be considered in the context of the quickly evolving sequencing technologies like hiseq  <dig> and hiseq  <dig> which can produce up to  <dig> million reads per sequencing lane.
in an evaluation of coverage depth of the chicken transcriptome, wang et al.
<cit>  find that while  <dig> million reads allow detection of 80% of the annotated genes, an increase from  <dig> to  <dig> million reads does not have a significant effect on transcriptome coverage or reliability of mrna measurements.
that said, current estimates of transcriptome coverage and the impacts of multiplexing strategies analysed in this paper assume unbiased sampling of transcripts.
it is highly likely that the power to detect de varies across transcripts with their sequence content, isoform complexity and abundance.
fang and cui
<cit>  warn against and discuss several sequencing biases that could create the need for high sequencing coverage to accurately estimate transcript abundance and variation.
the authors mention the importance of choosing whether to increase the sequencing depth per sample or to increase the number of biological replicates when planning an experiment.
here, we quantitatively argue that given a fixed budget, the benefits of increasing the number of biological replicates outweigh the corresponding decrease of sequencing depth.
this suggestion is backed by the patterns in figures
<dig> and
<dig> in which for a given number of n-biological replicates tpr drops very slowly as depth decreases, fpr remains low when sequencing depth is decreased.
in the light of new sequencing technologies rapidly increasing the available sequencing depth per lane, the inforrmation provided by biological replicates’ variation is likely to become a priority over sequencing depth.
CONCLUSIONS
not surprisingly, our results indicate that more biological replicates are needed to improve the quality and reliability of de detection using rna-seq.
importantly however, we also find that de analysis with rna-seq is robust to substantial loss of sequencing data as indicated by a slow decline in tpr accompanied by no increase in fpr.
our simulations strongly support that multiplexing experimental designs improve tpr and fpr while greatly reducing the cost of the experiment, as the benefits of multiplexing n-biological replicates far outweigh the decrease of available data per sample by 1/n.
as many available packages for de analysis are increasingly becoming faster and easier to use, our recommendation for most rna-seq de experiments is to use  <dig> different packages for de testing.
additional file
2: figure s <dig> illustrates the detection overlap between deseq and edger for two contrasting choices of sequencing depths and n-biological replicates: n =  <dig> at 100% depth and n =  <dig> at 25% depth.
the combined use of packages based on different distribution statistics or a different set of assumptions could generate useful information about a possible bias susceptibility of a given package particular to the specific dataset of interest.
to our knowledge, this is the most up-to-date comparison of deseq and edger’s performance relative to ability to detect de in a range of experimental designs.
it directly tests the efficiency of modern multiplex experimental design strategies.
our study informs important experimental design decisions now relevant when trying to maximise an rna-seq study to reliably detect de.

 rna sequencing allows an entire transcriptome to be surveyed at single-base resolution whilst concurrently profiling gene expression levels on a genome scale. rna-seq is an attractive approach as it profiles the transcriptome directly through sequencing and therefore does not require prior knowledge of the transcriptome under consideration. an example of the use of rna-seq as a high-resolution exploratory tool is the discovery of thousands of additional novel coding and noncoding genes, transcripts and isoforms of known genes despite the prior extensive annotation of the mouse and human genomes. arguably, the most popular use of rna-seq is profiling of gene expression or transcript abundance between samples or differential expression. the efficiency, resolution and cost advantages of using rna-seq as a tool for profiling de has prompted many biologists to abandon microarrays in favour of rna-seq. despite the advantages of using rna-seq for de analysis, there are several sources of sequencing bias and systematic noise that need to be considered when using this approach. clearly, rna-seq analysis is vulnerable to the general biases and errors inherent in the next-generation sequencing technology upon which it is based. these errors and biases include: sequencing errors, biases in sequence quality, nucleotide composition and error rates relative to the base position in the read, variability in sequence depth across the transcriptome due to preferential sites of fragmentation, variable primer and transcript nucleotide composition effects and finally, differences in the coverage and composition of raw sequence data generated from technical and biological replicate samples. recently, there have been several investigations into the biases that affect the accuracy with which rna-seq represents the absolute abundance of a given transcript as measured by high precision approaches such as taqman rt-pcr. it has been shown that these abundance measures are prone to biases correlated with the nucleotide composition and length of the transcript. several within and between sample correction and normalisation procedures have recently been developed to address these biases either as nucleotide composition effects or various combinations of nucleotide, length or library preparation biases. these approaches all yield improvements in the correspondence of rna-seq read counts with expression estimates gained by other experimental approaches. despite the known biases, rna-seq continues to be widely and successfully used to profile relative transcript abundances across samples to identify differentially expressed transcripts. the profile of a given transcript across a biological population would be hoped to be less prone to nucleotide composition and length biases as these variables remain constant. nevertheless, to accurately detect de across samples it is necessary to understand the sources of variation across technical and biological replication and where possible respond to these with an appropriate experimental design and statistically robust analysis. to date, there has been little discussion in the literature of efficient experimental designs for the detection of de and a lack of consensus about a standard and comprehensive approach to counter the many sources of noise and biases present in rna-seq has meant that some of the biological community remain sceptical about its reliability and unsure of how to design cost-efficient rna-seq experiments. good experimental design and appropriate analysis is integral to maximising the power of any ngs study. with regard to rna-seq, important experimental design decisions include the choice of sequencing depth and number of technical and/or biological replicates to use. for researchers with a fixed budget, often a critical design question is whether to increase the sequencing depth at the cost of reduced sample numbers or to increase the sample size with limited sequencing depth for each sample. sequencing depth sequencing depth is usually referenced to be the expected mean coverage at all loci over the target sequence, in the case of rna-seq experiments assuming all transcripts having similar levels of expression. without the benefit of extensive previous rna-seq studies, it is difficult in most cases to estimate prior to data generation the optimal sequencing depth or amount of sequencing data required to adequately power the detection of de in the transcriptome of interest. pragmatically, rna-seq sequencing depth is typically chosen based on an estimation of total transcriptome length and the expected dynamic range of transcript abundances. given the dynamic nature of the transcriptome, the suitability of these estimates could vary substantially across organisms, tissues, time points and biological contexts. wang et al. found a significant increase in correlation between gene transcripts observed and number of sequence reads generated when increasing sequencing depth from to million reads after which the gains plateau – million reads detected about 80% of the annotated chicken transcripts. despite the expectation of continuous sequencing depth increases in the near future, Łabaj et al. argue that most of the additional reads will align to the subset of already extensively sampled transcripts. as a result, transcripts with low to moderate expression levels will remain difficult to quantify with good precision using current rna-seq protocols even at higher read depths. greater sequencing depth will also increase sensitivity to detect smaller changes in relative expression, however this does not guarantee that these changes have functional impact in the biological system under study as opposed to tolerated fluctuations in transcript abundance. ideally, an efficient experimental design will be informed by an understanding of when increasing sequencing depth begins to provide rapidly diminishing returns with regard to transcript detection and de testing. replication replication is vital for robust statistical inference of de. in the context of rna sequencing, multiple nested levels of technical replication exist depending upon whether it is the sequence data generation, library preparation or rna extraction technical processes that are being replicated from the same biological sample. several published studies have incorporated technical replicates into their rna-seq experimental designs. the degree of technical variation present in these datasets appears to vary and the main source of technical variation appears to be library preparation. biological replication measures variation within the target population and simultaneously can counteract random technical variation as part of independent sample preparation. it has been shown that power to detect de improves when the number of biological replicates n is increased from n = to n = 5, however, to date few studies have incorporated extensive biological replication and extensive testing of the effects of replication on power is needed. more recently with the increasing utility and availability of multiplex experimental designs, the incorporation of biological replicates with decreased sequencing depth is becoming a much more attractive and cost-effective strategy. the relative merits of sacrificing sequencing depth for increased replication has not been rigorously explored. efficient experimental design multiplexing is an increasingly popular approach that allows the sequencing of multiple samples in a single sequencing lane or reaction and consequently the reduction in sequencing costs per sample. multiplexing uses indexing tags, “barcodes” or short stretches of sequence that are ligated to the start of sample sequence fragments during the library preparation step. barcodes are distinct between sample libraries and allow pooling for sequencing followed by allocation of reads back to individual samples after sequencing by analysis of the sequenced barcode. multiplex barcode designs are routinely available with up to samples in the same lane, recently up to yeast dna samples were profiled in single lane. novel methods are continuing to emerge for low-cost strategies to multiplex rna-seq samples. with the dramatic increases in sequencing yields being achieved with current chemistries and new platforms, multiplexing is becoming the method of choice to increase sample throughput. these designs have direct impacts on sequencing depth generated that need to be considered in the power of the experimental design. also, when multiplex strategies are used, biologists need to be mindful of potential systematic variations between sequencing lanes. these variations can be addressed through randomisation or blocking designs to distribute samples across lanes, see for a discussion of barcoding bias in multiplex sequencing, and for an alternative to barcoding. in a comparison between microarray and ngs technologies in synthetic pools of small rna, willenbrock et al. found that multiplexing resulted in decreased sensitivity due to a reduction of sequencing depth and a loss of reproducibility; however the authors did not investigate power for detection of de in their study. approach improving detection of de requires not only an appropriate experimental design but also a suitably powered analysis approach. several algorithms have recently been developed specifically to appropriately handle expected technical and biological variation arising from rna-seq experiments. a non-exhaustive list of these algorithms is: edger, deseq, nbpseq, bbseq, fdm, rsem, noiseq, myrna, cuffdiff. a thorough comparison of these packages’ performance with datasets of different properties falls beyond the scope of this study, however before considering issues relating to power and experimental design, it is important to investigate whether packages for de analysis give the correct type i error rate under the null hypothesis of no de. to do this evaluation we considered three popular packages for de analysis of rna-sequencing data. these packages are based on a negative binomial distribution model of read counts and include edger, deseq and nbpseq. to quantify the effects of different sequencing depths and replication choices we compared a range of realistic experimental designs for their ability to robustly detect de. using simulated data with known de transcripts allowed us to estimate the false positive rate and true positive rate of de calls. the changes of these rates were used to compare the detection power yielded by each choice of number of biological replicates and sequencing depth. in the methods section, we outline the definitions used for fpr and tpr as well as explaining the method used for the construction of the synthetic data; which includes induced differential expression, simulates the variations that biological replicates introduce and simulates loss of sequencing depth. in our study, we test a wide range of real-world experimental design scenarios for performance under the null hypothesis and in the presence of de. in these scenarios both the numbers of biological replicates n and the sequencing depth are varied. this provides a comprehensive quantitative comparison of different experimental design strategies and is particularly informative for those accessing modern multiplex approaches.  comparisons of statistical methods: edger, deseq, and nbpseq using simulated data under the null to test the performance of each package under the null hypothesis, we simulated sets of n “control” and n “treatment” lanes of counts in accordance with the procedure described in the methods section, for a range of values of n and with no de between treatments. for each value of n and for each package the simulation and testing were repeated times. figure shows the percentage of transcripts reported as differentially expressed at the 1% significance level by each of the three packages for a range of values of n. the height of each bar is the median value obtained from repetitions of the synthetic data generation with its associated 90% confidence intervals. under the null hypothesis, the percentage reported is the false positive rate defined by eq. and should match the significance level of α = 1% if the package is performing correctly. also shown are fprs for high-count transcripts and low-count transcripts. figure shows an example of the p-value distribution obtained for one experiment at n = biological replicates. ideally, p-values should have a uniform distribution in the interval if the package is performing correctly. immediately noticeable in the p-value histogram is a sharp spike in the right hand bin for low count transcripts, which is observed to be present in general for all values of n and all packages. this is a known artifact of calculating p-values for discrete random variables using the method described in and summarised in our methods subsection ‘under the null hypothesis’: when count sums in both conditions are equal the computed p-value is exactly the situation is most likely to occur for transcripts with extremely low counts, in which case it is difficult to draw meaningful conclusions regarding de via any method. the behaviour at the left hand end of the histogram, which drives the fprs plotted in figure varies considerably between packages and numbers of biological replicates. it is affected mainly by the method used for estimating a dispersion parameter ϕi for each transcript i. the package edger performs well for large numbers of biological replicates, for which squeezing of the dispersion estimate towards the common dispersion is minimal, and a tagwise estimate is appropriate. for small numbers of biological replicates, because the dispersion cannot be estimated accurately on a per-transcript basis, information is borrowed from the complete set of transcripts to squeeze the estimate towards a common dispersion estimate. for the high-count transcripts in particular, the squeezing causes the dispersion of the most highly dispersed transcripts to be underestimated, causing too many transcripts to be deemed differentially expressed, leading to an inflated fpr. in an effort to be conservative, deseq chooses as its estimate of dispersion the maximum of a per-transcript estimate and the functional form eq. which is fitted to the per-transcript estimates for all transcripts. our results indicate that the method performs well for the high-count transcripts when the number of biological replicates is small, but is otherwise over-conservative. this is generally to be preferred to an inflated fpr, as one has more evidence that what is called de is truly de. the package nbpseq imposes the functional relationship eq. which appears to be too restrictive for a number of relatively highly dispersed transcripts. for those transcripts the dispersion parameter is underestimated, leading to an overestimate of significance and hence an inflated fpr irrespective of the number of biological replicates. based on these results we selected deseq and edger for use in subsequent experimental design testing. throughout these tests, results obtained using deseq and edger are mostly compatible with each other. however, our comparison revealed a slightly inflated fpr from edger while deseq behaves more conservatively throughout. therefore, in the following section we will focus on the results obtained using deseq while the analogous results obtained with edger are presented in the additional file 1: figure s comparison of statistical methods: deseq and edger using simulated data with 15% de transcripts to test the performance of packages in the presence of an alternate hypothesis, we simulated sets of n “control” and n “treatment” lanes of counts with 15% of the transcripts either up- or down-regulated according to the procedure described in the methods section. all results presented from this point on are derived from deseq. detection of de as a function of number of biological replicates n with an increase in replication we saw a steady increase in the percentage de calls by deseq, increasing from % to % as n increased from to. as n increased, the fpr, defined by eq. at a significance level of α = 1%, remained below % for all values of n, and the tpr, defined by eq. with α = 1%, increased substantially from % to %. effects of biological replication on power to detect de using deseq. fpr and tpr are defined in eqs. & respectively at 1%. “call rate” is the total number of reported positives / the total number of transcripts. these values are also represented in figure at 100% sequencing depth. detection of de as a function of sequencing depth figure represents the combined results of decreasing sequencing depth for all values of n. it can be seen that as sequencing depth decreases the tpr generated by deseq decreases monotonically across all n while the fpr remains below %. 1n×100% sequencing depth shows the increase of tpr as more biological replicates n are used despite the loss power due to the sequencing depth reduction required by the multiplexing of lanes. this trend remains true even for the n = and n = cases. table shows the fpr for all biological replicates n and a subset of sequencing depths: 25%, 50%, 75% and 100%, the fpr remains below % at all sequencing depths. table shows the tpr reported by deseq for the same subset of sequencing depths, here the tpr increases strongly as sequencing depth increases for any number of biological replicates n. effects of sequencing depth on fpr values for a subset of our tested depths = 25%, 50%, 75% & 100%. effects of sequencing depth on tpr values for a subset of our tested depths = 25%, 50%, 75% & 100%. detection of de across multiplex experimental design strategies we simulated various scenarios of multiplexing n-control samples vs. n-treatment samples into two sequencing lanes – each control and treatment sample at a sequencing depth =1n×100%. in figures and a solid grey line connecting every value of n at its corresponding sequencing depth provides a summary of the performance of these multiplexing scenarios. we call this trend the “multiplex line” and it provides an insight into the results obtained by increasing the number of biological replicates used into a fixed number of sequencing lanes, in this case sequencing lanes. 1n×100% sequencing depth. the multiplex line in figure shows a clear increase in tpr as replication increases despite the loss of detection power that decreasing sequencing depth induces. it can also be seen that the fpr remains below % for all multiplex scenarios tested. note that for completeness we also added multiplex scenarios for n = &n = whose results follow the trends well. the multiplex line strongly favours adding more biological replicates despite the inherent loss of sequencing depth as shown by its dramatic positive slope for the tpr while maintaining a roughly constant, low fpr. fold-changes as indicators of de it is common practice among biologists to use fold-change, rather than p-values, as an indicator of de. figure shows results analogous to those of figure when the criterion of fold-change ≥ is used to detect de: as replication n increases, both tpr and fpr decrease because more biological replicates have the effect of averaging out differences between control and treatment lanes. note that, as sequencing depth decreases, the fpr increases owing to the growing number of transcripts with very low numbers of counts, in which case the poisson shot noise of the sequencer can easily induce a spurious doubling or halving of counts. this effect is ameliorated by adding count to all transcripts prior to de analysis – doing so, does not affect the calculation of p-values. discussion comparisons of de algorithms: edger, deseq and nbpseq our comparison of these three de detection algorithms under the null hypothesis revealed different performances when different numbers of biological replicates n, are used. deseq consistently performed more conservatively across the different n biological replicates scenarios. deseq’s performance was closest to the expected significance level when only using high-count transcripts while for only low-count transcripts over-conservative behaviour is shown. edger overestimates de detection for small values of n while its performance improves as n increases. edger’s level of detection is constant over n when only low-count transcripts are used while overestimation increases when only high-count transcripts are used. nbpseq overestimated detection across n for the three scenarios. this comparison led us to use both deseq and edger throughout our replication and sequencing depth simulations. we ultimately chose deseq’s resultsa as this package behaved slightly more conservatively and appeared less sensitive to changes in replication. in a study by tarazona et al., it is argued that negative-binomial based de analysis packages like deseq and edger are highly sensitive to sequence depth increases and are therefore unable to control the fpr as sequencing depth increases. tarazona et al. propose a non-parametric algorithm to calculate de based on a noise distribution created with fold-changes and the absolute differences between the transcript’s control and experiment lane counts. however, kvam and liu argue that due to the small number of replicates typically used for rna-seq experiments, non-parametric methods do not offer enough detection power and suggest that current statistical methods to detect de genes based on parametric models for rna-seq data remain a more adequate approach. in our study we also find that both deseq and edger tend to slightly increase the fpr as sequencing depth increases – as higher depths induce deseq and edger to assign smaller p-values to transcripts with small fold-changesb. however in no instance do we obtain a fpr larger than 1% for deseq –. it is worth mentioning that the latest updates to deseq and edger – released after the studies and reduce the number of false positive calls by about 50%. effects of replication for detection of de to quantify the effects of replication in rna-seq de experiments, we tested n-control vs. n-treatment biological replicates while maintaining sequencing depth constant. we find that as n increases both algorithms increase the call rate and tpr while the fpr remains unchanged. our results clearly support the simple message that more biological replicates are not only desirable but needed to improve the quality and reliability of de detection using rna-seq, however, due to the costs associated with rna-seq, many experiments are likely to need to use multiplex designs to achieve this level of replication. this study is concerned with the simulation of overdispersion effects due to biological variability and it is implied that overdispersion due to technical variability is nested within this estimation. it is worth mentioning that, while biological variability is important, the contribution to overdispersion by technical variation is not negligible, and disagreements between estimates of expression can occur at all levels of coverage. ideally, rna-seq experimental design with biological replication should also aim to block sources of technical variation, such as between lane variations, to constrain the dispersion of rna-seq experiments. effects of sequencing depth for detection of de to quantify the effects of sequencing depth in rna-seq de experiments, we simulated an extensive sequencing depth range for every case of n-control vs. n-treatment biological replicates. as the amount of available sequencing data is decreased, both packages decrease the call rate and tpr while the fpr remains low. tpr decreases very slowly as sequencing depth decreases, suggesting that sequencing depth can be reduced to ∼15% without much impact on tpr. we conclude that de analysis with rna-seq is robust to substantial loss of sequencing data as indicated by a slow decline in tpr as sequencing depth is lost accompanied by no increase in fpr. these findings seem consistent with the results reported by bashir et al. who observed that lower levels of transcriptome sequencing had sufficient information to estimate the distribution of expression values arising from observed transcripts. bashir et al. did not directly test power to detect de, however as testing for de relies on good concordance with the expected distribution, it follows that de is reasonably robust to loss of sequencing data. multiplexing experimental designs to quantify the effects of varying both n and sequencing depth, we simulated multiplexing n-control vs. n-treatment lanes into two sequencing lanes. we observed a steady increase in tpr with the increase in n, despite the corresponding decrease in sequencing data per transcript by 1/n. similarly, for both deseq and edger the number of de calls and the tpr increases with n, as we observed previously and is unaffected by the decrease in data. for deseq, the fpr remains roughly constant and always below %, while for edger, the fpr decreases slowly as n increases. our simulations strongly support that the benefits of multiplexing n-biological replicates into one sequencing lane, far outweigh the decrease of available data per sample by 1/n. these multiplexing experimental designs improve tpr and fpr while greatly reducing the cost of the experiment. while the detection of de appears robust to available sequence data, there remains the question of how multiplexing affects coverage of the transcriptome and detection of low abundant or rare transcripts. this coverage issue will increasingly be counterbalanced by rapid increases in data generation capacity from a single sequencing experiment. in a detailed study of the marioni human dataset, banshir et al., reported that over 90% of the total observed transcripts were sampled with million reads. this should be considered in the context of the quickly evolving sequencing technologies like hiseq and hiseq which can produce up to million reads per sequencing lane. in an evaluation of coverage depth of the chicken transcriptome, wang et al. find that while million reads allow detection of 80% of the annotated genes, an increase from to million reads does not have a significant effect on transcriptome coverage or reliability of mrna measurements. that said, current estimates of transcriptome coverage and the impacts of multiplexing strategies analysed in this paper assume unbiased sampling of transcripts. it is highly likely that the power to detect de varies across transcripts with their sequence content, isoform complexity and abundance. fang and cui warn against and discuss several sequencing biases that could create the need for high sequencing coverage to accurately estimate transcript abundance and variation. the authors mention the importance of choosing whether to increase the sequencing depth per sample or to increase the number of biological replicates when planning an experiment. here, we quantitatively argue that given a fixed budget, the benefits of increasing the number of biological replicates outweigh the corresponding decrease of sequencing depth. this suggestion is backed by the patterns in figures and in which for a given number of n-biological replicates tpr drops very slowly as depth decreases, fpr remains low when sequencing depth is decreased. in the light of new sequencing technologies rapidly increasing the available sequencing depth per lane, the inforrmation provided by biological replicates’ variation is likely to become a priority over sequencing depth.  not surprisingly, our results indicate that more biological replicates are needed to improve the quality and reliability of de detection using rna-seq. importantly however, we also find that de analysis with rna-seq is robust to substantial loss of sequencing data as indicated by a slow decline in tpr accompanied by no increase in fpr. our simulations strongly support that multiplexing experimental designs improve tpr and fpr while greatly reducing the cost of the experiment, as the benefits of multiplexing n-biological replicates far outweigh the decrease of available data per sample by 1/n. as many available packages for de analysis are increasingly becoming faster and easier to use, our recommendation for most rna-seq de experiments is to use different packages for de testing. additional file 2: figure s illustrates the detection overlap between deseq and edger for two contrasting choices of sequencing depths and n-biological replicates: n = at 100% depth and n = at 25% depth. the combined use of packages based on different distribution statistics or a different set of assumptions could generate useful information about a possible bias susceptibility of a given package particular to the specific dataset of interest. to our knowledge, this is the most up-to-date comparison of deseq and edger’s performance relative to ability to detect de in a range of experimental designs. it directly tests the efficiency of modern multiplex experimental design strategies. our study informs important experimental design decisions now relevant when trying to maximise an rna-seq study to reliably detect de. 
with regard to testing of various experimental designs, this work strongly suggests that greater power is gained through the use of biological replicates relative to library replicates and sequencing depth. rna sequencing has emerged as a powerful approach for the detection of differential gene expression with both high-throughput and high resolution capabilities possible depending upon the experimental design chosen. multiplex experimental designs are now readily available, these can be utilised to increase the numbers of samples or replicates profiled at the cost of decreased sequencing depth generated per sample. these datasets were used to evaluate the performance of three commonly used differential expression analysis algorithms and to quantify the changes in power with respect to true and false positive rates when simulating variations in sequencing depth, biological replication and multiplex experimental design choices. this work quantitatively explores comparisons between contemporary analysis tools and experimental design choices for the detection of differential expression using rna-seq. this study presents a detailed analysis of the power to detect differential expression in a range of scenarios including simulated null and differential expression distributions with varying numbers of biological or technical replicates, sequencing depths and analysis methods. these strategies impact on the power of the approach to accurately identify differential expression. differential and non-differential expression datasets were simulated using a combination of negative binomial and exponential distributions derived from real rna-seq data.
 rna sequencing has emerged as a powerful approach for the detection of differential gene expression with both high-throughput and high resolution capabilities possible depending upon the experimental design chosen. multiplex experimental designs are now readily available, these can be utilised to increase the numbers of samples or replicates profiled at the cost of decreased sequencing depth generated per sample. these strategies impact on the power of the approach to accurately identify differential expression. this study presents a detailed analysis of the power to detect differential expression in a range of scenarios including simulated null and differential expression distributions with varying numbers of biological or technical replicates, sequencing depths and analysis methods.  differential and non-differential expression datasets were simulated using a combination of negative binomial and exponential distributions derived from real rna-seq data. these datasets were used to evaluate the performance of three commonly used differential expression analysis algorithms and to quantify the changes in power with respect to true and false positive rates when simulating variations in sequencing depth, biological replication and multiplex experimental design choices.  this work quantitatively explores comparisons between contemporary analysis tools and experimental design choices for the detection of differential expression using rna-seq. we found that the deseq algorithm performs more conservatively than edger and nbpseq. with regard to testing of various experimental designs, this work strongly suggests that greater power is gained through the use of biological replicates relative to library replicates and sequencing depth. strikingly, sequencing depth could be reduced as low as 15% without substantial impacts on false positive or true positive rates. rna-seqdifferential expression analysissequencing depthreplicationexperimental designmultiplex 