BACKGROUND
transcriptome studies usually aim at understanding inducible biological functions through an analysis of differential gene expression .
since relatively recently, the variation in gene expression is commonly studied through rna sequencing or rna-seq, based on next generation sequencing  technologies.
in these study approaches, de is usually inferred from comparison of two different treatments, developmental stages, or different tissues.
a key step in these analyses is the reliable identification of significant de.
most current statistical approaches employ a probabilistic model, such as the negative binomial  , poisson  <cit> , the generalized poisson  model  <cit> , and use information on gene expression variation in the data to account for ambiguity caused by sample size, biological and technical biases, overall levels of expression and the presence of outliers.
de inference is usually based on the null hypothesis that the means of read counts among conditions are the same or follow the same distribution.
these tests neglect the magnitude of encountered differences and might report statistically highly significant de with arbitrarily small fold change, at least if the number of sequencing counts is large enough  <cit> .
however, small fold changes may represent artifacts and often cannot be validated experimentally .
thus, they might not be worth further investigation.
a currently common solution is sought by combining the statistical indication  with a specified minimum fold change  <cit> .
this approach has the possible problem of a high number of identified candidate genes with low count numbers  and its dependence on an arbitrarily chosen fold-change cut-off value.
an alternative approach has so far only been established for chip-seq data and relies on an analysis of count differences between test and reference conditions  <cit> .
in this case, the statistics are based on a measure that considers the magnitude of count differences and the level of expression variation across replicates with the effect that genes with only minor expression levels and only small fold change are selected against.
in consideration of such potential advantages, such an approach may prove useful for reliable de identification in rna-seq data.
here, we introduce absseq , which employs an nb distribution to model count differences between conditions.
it permits testing the magnitude of observed count differences taking into consideration background expression level variation.
in particular, absseq accounts for heterogeneous dispersions in expression level across genes by adding expected values  to reads count according to the smoothed mean-variance relationship  <cit> , which thus adjusts parameters in the nb distribution .
in addition, absseq imposes a penalty on the dispersion estimation, it uses a new outlier detection strategy, and it also inroduces a procedure for shrinkage of fold change to disfavor identification of candidate genes with abnormal high dispersions and extremely low expression.
using real and simulated datasets, we demonstrate that our method is highly efficient in reducing the false discovery rate  and thus in identifying truly differentially expressed genes in rna-seq data.
it therefore shows an at least similar performance than several frequently used, alternative approaches like those implemented in the software packages deseq  <cit> , deseq <dig>  <cit> , edger  <cit>  , limma  <cit>  , bayseq  <cit> , and ebseq  <cit> .
implementation
absseq has been implemented in the software package absseq for the cross-platform environment r  <cit> .
absseq is released under the gpl- <dig> license as part of the bioconductor project  <cit>  at url: http://bioconductor.org/packages/devel/bioc/html/absseq.html.
RESULTS
we firstly introduce our approach with the help of the modencodefly and abrf datasets .
thereafter, performance of our method is compared with that of several previously developed and currently popular methods , including one, ebseq, which allows to evaluate de at both transcript and also gene level  <cit> .
we exlcude cuffdiff <dig>  <cit>  from our assessment because it was previously compared with the other available approaches and generally found to produce higher rates of false positives without an increase in sensitivity  <cit> .
method evaluation is based on two types of data sets.
on the one hand, we use simulated data, for which data structure can be efficiently controlled and which have been widely used to evalute methods of differential expression analysis .
we use the same strategy and identical simulated data sets as soneson et al.
<cit>  and compare method performance according to two criteria:  the ability to control type i error rates; and  the ability to rank truly de genes ahead of non-de ones.
on the other hand, we also evaluate our approach with the help of real data sets, as described in more detail below.
control of type i error rate
minimizing the type i error rate  or false positive rate is a primary goal of differential expression analysis  <cit> .
type i error is often introduced by under-estimation of disperson in rna-seq data and occurs at genes with very low or high counts  <cit> .
we thus compare the ability of the alternative approaches to control type i error rates, using two real data sets and also the simulated data sets from soneson et al.
<cit> .
de genes are defined by a p-value cutoff of  <dig>  for each method except bayseq and ebseq, which are excluded from this comparison since they report de by posterior probalilities instead of a p-value.
the simulated datasets are assumed to lack de genes, facilitating computation of the type i error rate by dividing the number of de genes identified by each method with the total number of genes.
figure  <dig> summarizes the results from the modencodefly data set  and two different simulation settings , including data sets of various replicate sample sizes and, in each case, ten independent repetitions .
additional file  <dig> shows the results for the abrf data set.fig.
<dig> method-dependent variation in type i error.
type i error rates for absseq and five alternative methods using the modencodefly real data set  and two simulation settings: negative binomial , and nb with random outliers .
b points show the absolute log fold change  distribution of false positives against the expression level .
c each boxplot summarizes the type i error rates across  <dig> independent simulated data sets.
asterisk indicates a statistically significant difference in type i error between absseq and any of the other methods.
n indicates the number of rna replicates considered in each case .
under all conditions, absseq reduced the type i error rate
the first comparison is based on a real data set for the fruitfly drosophila melanogaster, the modencodefly data set  <cit> , which characterizes the developmental transcriptome across  <dig> distinct stages  with technical replicates ranging from  <dig> to  <dig>  we randomly select  <dig> replicates for each condition and separate them into two groups, which should thus only be characterized by stochastic variations but not true de.
the results of our analysis is summarized in fig.
1a.
at the p-value cutoff of  <dig> , absseq identifies an average of  <dig> de genes and thus significantly fewer de genes than all alternative methods .
deseq <dig> also performs well on this real data set, while the highest type i error rate is obtained for limma .
next, we examined the distribution of false positives along absolute log <dig> fold change and expression level  using the data from fig.
1a.
as shown in fig.
1b, false positives with low logcpm  tend to have a high fold change , and vice versa.
this skewed distribution is very similar to the quadratic mean-variance relationship  <cit> , suggesting there might be a general under-estimation of variance or dispersion for these methods.
in contrast, absseq and deseq <dig> both shrink the fold change according to variance.
as a consequnce, they both exhibit a pronnounced reduction of false positives at low expression .
morever, absseq also reduces false positives at high expression level , which likely have a very low smoothed dispersion  <cit>  and are often inferred to be highly statistically significant but show only very small fold change.
as the modencodefly data set only allows us to consider two replicates per group and condition, the resulting statistical power may be limited.
therefore, we repeated this assessment using another real data set, the abrf data set  <cit>  , which is based on an rna-seq analysis of the same two samples across three independent laboraties and thus comprises for each sample three replicates that should only show variation caused by differences among the laboratories such as library preparation methods or sample processing procedures .
the analysis of the abfr data set confirms the previous results.
it demonstrates that absseq produces the smallest number of false positives and especially reduces fold change for the genes with generally low expression.
overall, the results from the two real data sets suggest that absseq has the ability to handle very small expression changes by considering the magnitude of absolute differences and penalizing the estimated dispersion .
our results also suggest that the alternative methods should allow enhanced reduction of the type i error rate if combined with additional filtering approaches, such as usage of a fold-change cut-off as discussed in  <cit>  and also further below.
in addition to the two real data sets, we also compare the ability of the alternative approaches to control type i error rates on simulated data .
generally, all methods are able to control type i error rate under  <dig>  when applied on the nb distributed data  but exhibit high diversity on the nb distributed data with randomly introduced outliers .
as already highlighted in  <cit> , deseq has excellent power to control type i error rates on nb <dig>  the performance of voom is relatively unaffected by sample size and outliers, implying advantages of log-transformation on dealing with high value outliers.
in contrast, edger does not control type i error rates efficiently when applied on data with outliers.
since both deseq <dig> and edger-robust integrate strategies to handle outliers, they expectedly reduce the type i error rate on s <dig>  <dig>  especially when compared to the earlier program versions .
absseq performs best in both cases , but slightly decreases its performance with increasing sample size .
taken together, absseq is able to efficiently control type i error rates for the real and simulated data sets  and it also reduces type i error rate at both low and high expression levels.
in addition, outliers impact the ability of controlling type i error rate for most methods except absseq and voom, which might be caused by shrinkage of the observed dispersion  or replacing the observed with a smoothed dispersion .
in contrast, absseq uses the observed dispersion directly, apparently enhancing control of type i errors to a rate of below  <dig> .
discrimination of de versus non-de genes in simulation studies
an ideal de inference method should be more sensitive to de than non-de genes, that is, it should be able to discriminate true de genes against non-de ones.
here, we evaluate the discriminative power of absseq and other selected methods in terms of the true and false positive rates and also the area under receiver operating characteristic  curve , using again the simulated data and general approach of soneson et al.
<cit> .
the auc was shown repeatedly to be informative as a measure of the overall discriminative performance of a method .
in particular, for our comparison, we extract a set of genes from the simulated data set using a given p-value or posterior probability  threshold.
thereafter, the obtained genes are divided into a truly positive group and a truly negative group according to pre-defined de genes in the simulated data.
this information then allows us to calculate the true positive and the false positive rate for all possible thresholds, construct roc curves and compute aucs using the roc package in bioconductor  <cit> .
for all simulations, we choose 10 % of the  <dig>  genes as de and symmetrically divide them into up- and down-regulated genes .
we summarize the results using boxplots for four different simulation settings, including data sets with various replicate sample sizes and, in each case, ten independent repetitions .fig.
<dig> auc comparison on simulated data.
area under the curve  for absseq and five alternative methods under two simulation settings: negative binomial  and nb with random outliers .
each boxplot summarizes the aucs across  <dig> independently simulated data sets.
asterisk indicates a statistically significant difference in auc between absseq and any of the other methods.
n indicates the number of considered rna-seq replicates, from  <dig> to  <dig>  under all conditions, absseq is highly effective in correctly identifying differentially expressed genes
when applied on the data set simulated using the nb distribution , absseq always performs at least as good as the alternative methods at the considered replicate sample sizes .
ebseq performs worse than the other approaches when applied on data with small sample size .
the performance of absseq and the other methods are generally improved as the sample size increases, revealing a positive power of sample size on identifying true des.
overall, these results suggest that our nb model fits the over-dispersion data at least as well as the nb model implemented in other methods.
we next test the influence of outliers, which we introduce into the nb distributed data using a similar approach as above  and which may show abnormally high counts, resulting in high fold changes and also false postives.
for these simulated data sets, absseq shows an advantage  at all replicate sizes, especially for the r <dig> data set  whose auc is even greater than  <dig>  at n =  <dig> .
this result indicates that absseq outlier detection is efficient.
interestingly, performance of the alternative methods also shows substantial variability.
for example, voom generally performs better at large sample size , but similar at small sample size with other methods; deseq <dig> performs better at n =  <dig> due to outlier detection but worse at n =  <dig> and n =  <dig> ; bayseq shows little improvement in performance as the sample size increases for the r <dig> data set; ebseq shows lowest auc at n =  <dig> and improves performance at large sample sizes ; edger-robust  shows an improved ability to handle outliers at small sample size  compared to edger.
overall, absseq is at least as good as alternative methods in discriminating between de and non-de genes, it is highly robust towards outliers at all sample sizes, while increasing the sample size improves the discriminative performance for all methods.
the high performance of absseq on the outlier data sets supports the efficiency of the implemented approach based on moderated median absolute deviation  in outlier detection even at small sample size .
together with the results in fig.
<dig>  our model on count differences seems to perform at least as good as other models using nb distributed data.
differential expression analysis on qrt-pcr validated real data
as simulated data are by nature artificial, we further evaluate method performance on real data sets.
the first of these relate to the maqc study, for which rna-seq-identified de genes were validated by quantitative reverse transcription pcr   <cit>  based on the commercially available taqman and primepcr methodologies.
although there is no single “gold standard” for assessment of rna-seq data reliability  <cit> , qrt-pcr based methods have widely been proposed and applied as a validation tool for de results from both microarray  <cit>  and rna-seq studies  <cit> .
here, we analysed two qrt-pcr validated data sets from the maqc study: the taqman data set from maqc-i, which included an assessment of a very small fraction of the total genes  and may thus be subject to biases, and additionally the primepcr data set from seqc , which covers more than  <dig>  genes  <cit> .
these two data sets were used to derive roc curves and auc measures for the compared analysis methods.
we consider this approach to provide at least an indication of the reliability and sensitivity of the analysis approach.
we here follow the general strategy from  <cit>  and  <cit>  and divide the taqman and primepcr gene sets into a de  group and a non-de  group based on whether their absolute log fold change  is larger or smaller than a defined threshold.
we use a logfc threshold of  <dig>   to derive roc curves.
the results for both data sets are essentially identical .
while the alternative methods can detect approximately half of the taqman validated de genes without false positives, absseq is even able to identify more than 75 % of the true de genes with a false positive rate of less than  <dig>  .
absseq reaches the highest auc of  <dig>  among six methods .
for the primepcr data set, the auc for each method decreases as the number of validated genes increases .
again, absseq performs best among all seven methods, supporting its ability to discriminate efficiently between de and non-de genes.fig.
<dig> comparison of methods using validated real data sets.
a-c based on data from the maqc study; d-e based on the abrf data set.
roc analysis for  taqman and  primepcr data sets at a qrt-pcr absolute log-ratio  threshold of  <dig> .
tpr, true positive rate; fpr, false positive rate.
absseq performs better than other methods in detecting true differential expression.
a gene was considered to be not differentially regulated if its logfc was less than  <dig> .
c minimal fold changes under various ajusted p-value cutoffs for the maqc ii data set.
d number of false postives in comparisons of samples from same condition but different lab sites and  number of de genes in comparison of samples from two conditons under additional filtering and confounding factor assessment approaches.
symbols in black show results from comparison of conditions from same laboratory and colored symbols those from comparison of conditions across laboratories.
genes are counted under  <dig> situations: orginal, without filtering ; foldchange, with a value greater than  <dig>  ; aveexp, with average logcpm greater than  <dig> ; combination of foldchange and aveexp ; and svaseq tested only for deseq <dig> and voom
analysis approaches, which do not consider the magnitude of expression differences, might yield highly statistically significant de for genes with only small fold change , which may however often be the result of chance.
the number of these type of de genes is usually not reduced by using an adjusted p-value cutoff in the alternative approaches, even if the cutoff is below  <dig> e- <dig>  therefore, other cutoff criteria are required such as fold change, which has the problem that the biologically relevant cutoff point is not clear.
the absseq-based analysis instead produces high correlation between the minimal fold change and the inferred adjusted p-value, indicating that the p-value alone will select against de genes with small fold change.
additional cutoff criteria therefore do not seem to be necessary for reliable de gene identification.
influence of cut-off criteria and confounding factor analysis procedures
we next investigate the influence of additional cut-off criteria on de detection with the help of the abrf data set, which is based on rna-seq data generated for the same sample in three different laboratories.
we apply the considered methods on this data set, which only contains variation caused by differences among the conisdered laboratories, such as biases during library preparation  <cit> , but not true de, thus allowing us to assess the efficacy of the methods to reduce the number of false positives .
in spite of varying numbers of detected de genes, absseq reports lowest number of false positives among all methods, irrespective of any additional filtering approach .
bayseq and ebseq also produce small numbers of false positives then compared to the remaining methods excluding absseq.
for all methods, the number of false positives reduces dramatically when filtered by foldchange  but less so when filtered by expression level .
this finding strongly suggests that a high foldchange cut-off increases power to control the false positive rate, yet with the problem that the choice of cut-off value will usually be arbitrary.
in addition, high specificity  might lead to low senstitivity .
to evaluate the ability of absseq to detect true positives, we apply absseq and alternative methods on the abrf data set whereby in this case we focus on the comparison of the two considered conditions  either within the considered laboratories  or across the laboratories .
the results are shown in fig.
3e .
all seven methods report similar numbers of de genes, especially after fold-change filtering.
this result indicates that absseq retains similar sensitivity than that shown by the alternative approaches.
confounding variation can originate from library preparation or other kinds of batch effects.
to remove its influence on de detection, it can be modeled and thus integrated into the statistical analysis  <cit> , as implemented in svaseq  <cit> .
to illustrate the possible influence of such variation, we applied svaseq togethor with deseq <dig> and voom.
svaseq together with voom is able to remove more than 50 % false positves for the abrf data set , in consistency with the previous application of the svaseq approach on data from the seqc study  <cit> .
however, when svaseq is combined with deseq <dig> it leads to only a small decrease in the number of false positives .
this result may suggest that the performance of svaseq depends on the de detection method itself and/or the linear model used in such methods.
moreover, the application of svaseq does not decrease sensitivity when combined with voom and only to a small extent when combined with deseq <dig> , suggesting that svaseq mainly improves removal of false positives but does not bias detection of true de.
in general, the usage of such confounding factor assessment procedures, including svaseq and also peer  <cit>  can help improve de detection.
yet, at the moment, its combination with the various de analysis methods is not straightforward, because both svaseq and peer produce non-integer values, whereas several of the current de analysis methods  rely on integer count data.
it thus represents a promising challenge to further develop these procedures as integrated modules of the common de detection methods.
assessment of statistical power via signal to noise ratio
to evaluate the statistical power of each method in measuring the magnitude of de in dependence of its variance, we repeated above comparison using genes that are exclusively expressed in only one condition of the maqc-ii data set following the approach from  <cit> .
the magnitude of de of genes expressed in only one condition is ideally shown as a signal-to-noise  ratio , which should be monotonically correlated with the p-value  <cit> .
a poor correlation between snr ratio and p-value might lead to reduced sensitivity  by assigning a large p-value to small snr ratio .
the monotonic dependency between predictor  and response  is inferred through an isotonic regression on  <dig> paired variables .
results are shown in fig.
<dig>  all methods but deseq and edger exhibit the desired monotonic behavior between snr ratio and adjusted p-value, in consistency with previous results from  <cit> .
two empirical bayes based approaches: bayseq and ebseq yield quite similar correlations between snr ratio and posterior probabilities.
in addition, voom assigns a more significant adjusted p-value for one specific gene with high snr ratio but low expression  whereas alternative methods produce adjusted p-values of around  <dig>  , suggesting voom is more sensitive to de at low expression level.
since deseq <dig> and voom test de on log fold change, we postulate that the closer correlation between snr ratio and adjusted p-value of absseq is due to modelling directly the magnitude of de difference.
overall, these results suggest that absseq seems to model the magnitude of count difference with higher accuracy, which might help de inference by reducing false positives.fig.
<dig> correlation between signal-to-noise ratio and p-value with true de present in only one condition.
evaluation is based on a total of  <dig> genes that are exclusively expressed in one condition in the maqc-ii data set.
gray points indicate genes with adjusted p-value value ≥  <dig> .
the data point highlighted by the green elipse refers to the gene with high signal-to-noise ratio but low expression.
the correlation is inferred using isotonic regression
differential expression analysis of real data with unbalanced designs
another real data set  is taken from  <cit> , consisting of  <dig> highly dispersed cdna samples from  <dig> females and  <dig> males.
de genes are inferred from male–female comparisons.
following  <cit> , a sensitivity analysis is predicted to find an over-representation of inferred de genes from the sex chromosomes.
indeed, the top ten de genes always include genes from sex chromosomes .
all methods except absseq and voom identify de genes beyond sex chromosomes.
this may indicate that absseq and voom retain higher specificity than the remaining methods and that alternative methods may not well model variance introduced by unequal smaple sizes.
ebseq produces the lowest number of de genes from sex chromosomes but the highest number from autosomes, confirming the previously observed lack of power of this method for the analysis of data with such high dispersion and uneven sample sizes  <cit> .
given the unequal sample size in this data set, the similar performance of absseq to that of alternative methods also suggests that our model is able to handle unequal sample sizes and high dispersion.
in particular, in absseq, we attempted to compensate for unequal sample size by adding expected reads counts to the smaller group until sample sizes are equal.
we always take the mean reads count of the small group as expected count, in order to minimize possible biases in subsequent variance estimations .
this compensation step is likely crucial for unbalanced data designs, especially in case of even larger differences than in our test data set.
in the future, it may be worth exploring in more detail alternative compensation procedures.table  <dig> number of de genes from sex chromosomes detected by each method in the hapmap-ceu data set at fdr-ajusted p-value of  <dig>
moderating fold change
fold change often serves as a more informative indicator for biologists to identify des.
it is also utilized in gene ranking to select candidates for further investigation and visualization .
however, the fold change neglects variance across samples and might not necessarily be informative, especially for genes with low counts .
to overcome this problem, deseq <dig> introduces an empirical bayes shrinkage for fold change estimation, which moderates the log fold change according to gene-specific dispersion  <cit> .
fold change can also be represented as a function of absolute count differences , suggesting a potential moderation of fold change via counts difference .
therefore, we introduce a fold change shrinkage procedure according to count differences and dispersion.
figure  <dig> shows how it works using the bottomly data set  <cit> .
genes with small counts tend to have high raw fold changes , which constrains reliable gene ranking by fold change at dynamic expression level.
shrinking fold change by adding pseudocounts according to expression level  removes this trend .fig.
<dig> moderation of log <dig> fold change.
a raw data  of the bottomfly study.
b the same data corrected by expression level.
c the same data corrected by expression level and gene-specific dispersion.
de genes  are shown in red.
non-de genes with high log <dig> fold change are marked by green elipses
however, this shrinkage approach neglects the gene-specific dispersion and thus shows no effects on non-de genes with high dispersion as well as high expression .
after taking account of gene-specific dispersion , the fold changes approximately reflect de genes  and produce nearly evenly distributed fold change values, apparently improving gene ranking and visualization.
notably, shrinkage with gene-specific dispersion only influences a small part of genes with high dispersion.
unlike the approach in deseq <dig>  our shrinkage method on gene-specific dispersion is based on p-values and therefore does not change the number of inferred significant de genes.
in practice, users can obtain all three types of fold change values  in absseq.
CONCLUSIONS
here we introduce a new method for differiential expression analysis of rna-seq count data, absseq.
distinct from other current methods, absseq infers de genes through the absolute differences in gene expression and assumes the differences to be influenced by two sources of variation: that found for average gene expression levels and that found for the magnitude of differential expression.
our approach employs a nb distribution to model these two parts and, as a consequence, it is able to detect de genes more effectively than existing methods, as demonstrated by our analysis of both real and simulated data.
in particular, absseq shows an advantage in discriminating de genes against non-de ones, it applies an efficient outlier detection approach and is thus robust against outliers.
morever, absseq inferred p-values correlate with the magnitude of count differences, thus producing a linear relationship between snr ratio and p-value.
as a result, it reduces type i error rates at both very low and high expression level and it also leads to a smaller number of highly significant de genes with small fold change.
in addition, absseq introduces a procedure to shrink fold change according to the smoothed dispersion across expression level and observed dispersion , which permits fold change comparisons across genes and thus might favor downstream analyses, such as gene set enrichment analysis by ranking  <cit> , clustering and visualization  or candidate selection.
a potential improvement of our approach in the future may be to adapt it to allow usage of more complex models which consider multiple conditions, its combination with additional normalization procedures, such as those implemented in peer and svaseq  <cit> , which can further help to filter out unwanted variation, and also adjustment of our approach to allow for analysis of de at the transcript level .
in summary, based on our analysis, we conclude that absseq represents a highly efficient approach for identification of significant des across a wide range of conditions and may help efficient downstram analysis of des.

 transcriptome studies usually aim at understanding inducible biological functions through an analysis of differential gene expression. since relatively recently, the variation in gene expression is commonly studied through rna sequencing or rna-seq, based on next generation sequencing technologies. in these study approaches, de is usually inferred from comparison of two different treatments, developmental stages, or different tissues. a key step in these analyses is the reliable identification of significant de. most current statistical approaches employ a probabilistic model, such as the negative binomial, poisson, the generalized poisson model, and use information on gene expression variation in the data to account for ambiguity caused by sample size, biological and technical biases, overall levels of expression and the presence of outliers. de inference is usually based on the null hypothesis that the means of read counts among conditions are the same or follow the same distribution. these tests neglect the magnitude of encountered differences and might report statistically highly significant de with arbitrarily small fold change, at least if the number of sequencing counts is large enough. however, small fold changes may represent artifacts and often cannot be validated experimentally. thus, they might not be worth further investigation. a currently common solution is sought by combining the statistical indication with a specified minimum fold change. this approach has the possible problem of a high number of identified candidate genes with low count numbers and its dependence on an arbitrarily chosen fold-change cut-off value. an alternative approach has so far only been established for chip-seq data and relies on an analysis of count differences between test and reference conditions. in this case, the statistics are based on a measure that considers the magnitude of count differences and the level of expression variation across replicates with the effect that genes with only minor expression levels and only small fold change are selected against. in consideration of such potential advantages, such an approach may prove useful for reliable de identification in rna-seq data. here, we introduce absseq, which employs an nb distribution to model count differences between conditions. it permits testing the magnitude of observed count differences taking into consideration background expression level variation. in particular, absseq accounts for heterogeneous dispersions in expression level across genes by adding expected values to reads count according to the smoothed mean-variance relationship, which thus adjusts parameters in the nb distribution. in addition, absseq imposes a penalty on the dispersion estimation, it uses a new outlier detection strategy, and it also inroduces a procedure for shrinkage of fold change to disfavor identification of candidate genes with abnormal high dispersions and extremely low expression. using real and simulated datasets, we demonstrate that our method is highly efficient in reducing the false discovery rate and thus in identifying truly differentially expressed genes in rna-seq data. it therefore shows an at least similar performance than several frequently used, alternative approaches like those implemented in the software packages deseq, deseq, edger, limma, bayseq, and ebseq. implementation absseq has been implemented in the software package absseq for the cross-platform environment r. absseq is released under the gpl- license as part of the bioconductor project at url: http://bioconductor.org/packages/devel/bioc/html/absseq.html.  we firstly introduce our approach with the help of the modencodefly and abrf datasets. thereafter, performance of our method is compared with that of several previously developed and currently popular methods, including one, ebseq, which allows to evaluate de at both transcript and also gene level. we exlcude cuffdiff from our assessment because it was previously compared with the other available approaches and generally found to produce higher rates of false positives without an increase in sensitivity. method evaluation is based on two types of data sets. on the one hand, we use simulated data, for which data structure can be efficiently controlled and which have been widely used to evalute methods of differential expression analysis. we use the same strategy and identical simulated data sets as soneson et al. and compare method performance according to two criteria: the ability to control type i error rates; and the ability to rank truly de genes ahead of non-de ones. on the other hand, we also evaluate our approach with the help of real data sets, as described in more detail below. control of type i error rate minimizing the type i error rate or false positive rate is a primary goal of differential expression analysis. type i error is often introduced by under-estimation of disperson in rna-seq data and occurs at genes with very low or high counts. we thus compare the ability of the alternative approaches to control type i error rates, using two real data sets and also the simulated data sets from soneson et al.. de genes are defined by a p-value cutoff of for each method except bayseq and ebseq, which are excluded from this comparison since they report de by posterior probalilities instead of a p-value. the simulated datasets are assumed to lack de genes, facilitating computation of the type i error rate by dividing the number of de genes identified by each method with the total number of genes. figure  summarizes the results from the modencodefly data set and two different simulation settings, including data sets of various replicate sample sizes and, in each case, ten independent repetitions. additional file shows the results for the abrf data set.fig. method-dependent variation in type i error. type i error rates for absseq and five alternative methods using the modencodefly real data set and two simulation settings: negative binomial, and nb with random outliers. b points show the absolute log fold change distribution of false positives against the expression level. c each boxplot summarizes the type i error rates across independent simulated data sets. asterisk indicates a statistically significant difference in type i error between absseq and any of the other methods. n indicates the number of rna replicates considered in each case. under all conditions, absseq reduced the type i error rate the first comparison is based on a real data set for the fruitfly drosophila melanogaster, the modencodefly data set, which characterizes the developmental transcriptome across distinct stages with technical replicates ranging from to we randomly select replicates for each condition and separate them into two groups, which should thus only be characterized by stochastic variations but not true de. the results of our analysis is summarized in fig. 1a. at the p-value cutoff of, absseq identifies an average of de genes and thus significantly fewer de genes than all alternative methods. deseq also performs well on this real data set, while the highest type i error rate is obtained for limma. next, we examined the distribution of false positives along absolute log fold change and expression level using the data from fig. 1a. as shown in fig. 1b, false positives with low logcpm tend to have a high fold change, and vice versa. this skewed distribution is very similar to the quadratic mean-variance relationship, suggesting there might be a general under-estimation of variance or dispersion for these methods. in contrast, absseq and deseq both shrink the fold change according to variance. as a consequnce, they both exhibit a pronnounced reduction of false positives at low expression. morever, absseq also reduces false positives at high expression level, which likely have a very low smoothed dispersion and are often inferred to be highly statistically significant but show only very small fold change. as the modencodefly data set only allows us to consider two replicates per group and condition, the resulting statistical power may be limited. therefore, we repeated this assessment using another real data set, the abrf data set, which is based on an rna-seq analysis of the same two samples across three independent laboraties and thus comprises for each sample three replicates that should only show variation caused by differences among the laboratories such as library preparation methods or sample processing procedures. the analysis of the abfr data set confirms the previous results. it demonstrates that absseq produces the smallest number of false positives and especially reduces fold change for the genes with generally low expression. overall, the results from the two real data sets suggest that absseq has the ability to handle very small expression changes by considering the magnitude of absolute differences and penalizing the estimated dispersion. our results also suggest that the alternative methods should allow enhanced reduction of the type i error rate if combined with additional filtering approaches, such as usage of a fold-change cut-off as discussed in and also further below. in addition to the two real data sets, we also compare the ability of the alternative approaches to control type i error rates on simulated data. generally, all methods are able to control type i error rate under when applied on the nb distributed data but exhibit high diversity on the nb distributed data with randomly introduced outliers. as already highlighted in, deseq has excellent power to control type i error rates on nb the performance of voom is relatively unaffected by sample size and outliers, implying advantages of log-transformation on dealing with high value outliers. in contrast, edger does not control type i error rates efficiently when applied on data with outliers. since both deseq and edger-robust integrate strategies to handle outliers, they expectedly reduce the type i error rate on s especially when compared to the earlier program versions. absseq performs best in both cases, but slightly decreases its performance with increasing sample size. taken together, absseq is able to efficiently control type i error rates for the real and simulated data sets and it also reduces type i error rate at both low and high expression levels. in addition, outliers impact the ability of controlling type i error rate for most methods except absseq and voom, which might be caused by shrinkage of the observed dispersion or replacing the observed with a smoothed dispersion. in contrast, absseq uses the observed dispersion directly, apparently enhancing control of type i errors to a rate of below. discrimination of de versus non-de genes in simulation studies an ideal de inference method should be more sensitive to de than non-de genes, that is, it should be able to discriminate true de genes against non-de ones. here, we evaluate the discriminative power of absseq and other selected methods in terms of the true and false positive rates and also the area under receiver operating characteristic curve, using again the simulated data and general approach of soneson et al.. the auc was shown repeatedly to be informative as a measure of the overall discriminative performance of a method. in particular, for our comparison, we extract a set of genes from the simulated data set using a given p-value or posterior probability threshold. thereafter, the obtained genes are divided into a truly positive group and a truly negative group according to pre-defined de genes in the simulated data. this information then allows us to calculate the true positive and the false positive rate for all possible thresholds, construct roc curves and compute aucs using the roc package in bioconductor. for all simulations, we choose 10 % of the genes as de and symmetrically divide them into up- and down-regulated genes. we summarize the results using boxplots for four different simulation settings, including data sets with various replicate sample sizes and, in each case, ten independent repetitions.fig. auc comparison on simulated data. area under the curve for absseq and five alternative methods under two simulation settings: negative binomial and nb with random outliers. each boxplot summarizes the aucs across independently simulated data sets. asterisk indicates a statistically significant difference in auc between absseq and any of the other methods. n indicates the number of considered rna-seq replicates, from to under all conditions, absseq is highly effective in correctly identifying differentially expressed genes when applied on the data set simulated using the nb distribution, absseq always performs at least as good as the alternative methods at the considered replicate sample sizes. ebseq performs worse than the other approaches when applied on data with small sample size. the performance of absseq and the other methods are generally improved as the sample size increases, revealing a positive power of sample size on identifying true des. overall, these results suggest that our nb model fits the over-dispersion data at least as well as the nb model implemented in other methods. we next test the influence of outliers, which we introduce into the nb distributed data using a similar approach as above and which may show abnormally high counts, resulting in high fold changes and also false postives. for these simulated data sets, absseq shows an advantage at all replicate sizes, especially for the r data set whose auc is even greater than at n = . this result indicates that absseq outlier detection is efficient. interestingly, performance of the alternative methods also shows substantial variability. for example, voom generally performs better at large sample size, but similar at small sample size with other methods; deseq performs better at n =  due to outlier detection but worse at n =  and n = ; bayseq shows little improvement in performance as the sample size increases for the r data set; ebseq shows lowest auc at n =  and improves performance at large sample sizes; edger-robust shows an improved ability to handle outliers at small sample size compared to edger. overall, absseq is at least as good as alternative methods in discriminating between de and non-de genes, it is highly robust towards outliers at all sample sizes, while increasing the sample size improves the discriminative performance for all methods. the high performance of absseq on the outlier data sets supports the efficiency of the implemented approach based on moderated median absolute deviation in outlier detection even at small sample size. together with the results in fig.  our model on count differences seems to perform at least as good as other models using nb distributed data. differential expression analysis on qrt-pcr validated real data as simulated data are by nature artificial, we further evaluate method performance on real data sets. the first of these relate to the maqc study, for which rna-seq-identified de genes were validated by quantitative reverse transcription pcr based on the commercially available taqman and primepcr methodologies. although there is no single “gold standard” for assessment of rna-seq data reliability, qrt-pcr based methods have widely been proposed and applied as a validation tool for de results from both microarray and rna-seq studies. here, we analysed two qrt-pcr validated data sets from the maqc study: the taqman data set from maqc-i, which included an assessment of a very small fraction of the total genes and may thus be subject to biases, and additionally the primepcr data set from seqc, which covers more than genes. these two data sets were used to derive roc curves and auc measures for the compared analysis methods. we consider this approach to provide at least an indication of the reliability and sensitivity of the analysis approach. we here follow the general strategy from and and divide the taqman and primepcr gene sets into a de group and a non-de group based on whether their absolute log fold change is larger or smaller than a defined threshold. we use a logfc threshold of to derive roc curves. the results for both data sets are essentially identical. while the alternative methods can detect approximately half of the taqman validated de genes without false positives, absseq is even able to identify more than 75 % of the true de genes with a false positive rate of less than. absseq reaches the highest auc of among six methods. for the primepcr data set, the auc for each method decreases as the number of validated genes increases. again, absseq performs best among all seven methods, supporting its ability to discriminate efficiently between de and non-de genes.fig. comparison of methods using validated real data sets. a-c based on data from the maqc study; d-e based on the abrf data set. roc analysis for taqman and primepcr data sets at a qrt-pcr absolute log-ratio threshold of. tpr, true positive rate; fpr, false positive rate. absseq performs better than other methods in detecting true differential expression. a gene was considered to be not differentially regulated if its logfc was less than. c minimal fold changes under various ajusted p-value cutoffs for the maqc ii data set. d number of false postives in comparisons of samples from same condition but different lab sites and number of de genes in comparison of samples from two conditons under additional filtering and confounding factor assessment approaches. symbols in black show results from comparison of conditions from same laboratory and colored symbols those from comparison of conditions across laboratories. genes are counted under situations: orginal, without filtering; foldchange, with a value greater than; aveexp, with average logcpm greater than; combination of foldchange and aveexp; and svaseq tested only for deseq and voom analysis approaches, which do not consider the magnitude of expression differences, might yield highly statistically significant de for genes with only small fold change, which may however often be the result of chance. the number of these type of de genes is usually not reduced by using an adjusted p-value cutoff in the alternative approaches, even if the cutoff is below e- therefore, other cutoff criteria are required such as fold change, which has the problem that the biologically relevant cutoff point is not clear. the absseq-based analysis instead produces high correlation between the minimal fold change and the inferred adjusted p-value, indicating that the p-value alone will select against de genes with small fold change. additional cutoff criteria therefore do not seem to be necessary for reliable de gene identification. influence of cut-off criteria and confounding factor analysis procedures we next investigate the influence of additional cut-off criteria on de detection with the help of the abrf data set, which is based on rna-seq data generated for the same sample in three different laboratories. we apply the considered methods on this data set, which only contains variation caused by differences among the conisdered laboratories, such as biases during library preparation, but not true de, thus allowing us to assess the efficacy of the methods to reduce the number of false positives. in spite of varying numbers of detected de genes, absseq reports lowest number of false positives among all methods, irrespective of any additional filtering approach. bayseq and ebseq also produce small numbers of false positives then compared to the remaining methods excluding absseq. for all methods, the number of false positives reduces dramatically when filtered by foldchange but less so when filtered by expression level. this finding strongly suggests that a high foldchange cut-off increases power to control the false positive rate, yet with the problem that the choice of cut-off value will usually be arbitrary. in addition, high specificity might lead to low senstitivity. to evaluate the ability of absseq to detect true positives, we apply absseq and alternative methods on the abrf data set whereby in this case we focus on the comparison of the two considered conditions either within the considered laboratories or across the laboratories. the results are shown in fig. 3e. all seven methods report similar numbers of de genes, especially after fold-change filtering. this result indicates that absseq retains similar sensitivity than that shown by the alternative approaches. confounding variation can originate from library preparation or other kinds of batch effects. to remove its influence on de detection, it can be modeled and thus integrated into the statistical analysis, as implemented in svaseq. to illustrate the possible influence of such variation, we applied svaseq togethor with deseq and voom. svaseq together with voom is able to remove more than 50 % false positves for the abrf data set, in consistency with the previous application of the svaseq approach on data from the seqc study. however, when svaseq is combined with deseq it leads to only a small decrease in the number of false positives. this result may suggest that the performance of svaseq depends on the de detection method itself and/or the linear model used in such methods. moreover, the application of svaseq does not decrease sensitivity when combined with voom and only to a small extent when combined with deseq, suggesting that svaseq mainly improves removal of false positives but does not bias detection of true de. in general, the usage of such confounding factor assessment procedures, including svaseq and also peer can help improve de detection. yet, at the moment, its combination with the various de analysis methods is not straightforward, because both svaseq and peer produce non-integer values, whereas several of the current de analysis methods rely on integer count data. it thus represents a promising challenge to further develop these procedures as integrated modules of the common de detection methods. assessment of statistical power via signal to noise ratio to evaluate the statistical power of each method in measuring the magnitude of de in dependence of its variance, we repeated above comparison using genes that are exclusively expressed in only one condition of the maqc-ii data set following the approach from. the magnitude of de of genes expressed in only one condition is ideally shown as a signal-to-noise ratio, which should be monotonically correlated with the p-value. a poor correlation between snr ratio and p-value might lead to reduced sensitivity by assigning a large p-value to small snr ratio. the monotonic dependency between predictor and response is inferred through an isotonic regression on paired variables. results are shown in fig.  all methods but deseq and edger exhibit the desired monotonic behavior between snr ratio and adjusted p-value, in consistency with previous results from. two empirical bayes based approaches: bayseq and ebseq yield quite similar correlations between snr ratio and posterior probabilities. in addition, voom assigns a more significant adjusted p-value for one specific gene with high snr ratio but low expression whereas alternative methods produce adjusted p-values of around, suggesting voom is more sensitive to de at low expression level. since deseq and voom test de on log fold change, we postulate that the closer correlation between snr ratio and adjusted p-value of absseq is due to modelling directly the magnitude of de difference. overall, these results suggest that absseq seems to model the magnitude of count difference with higher accuracy, which might help de inference by reducing false positives.fig. correlation between signal-to-noise ratio and p-value with true de present in only one condition. evaluation is based on a total of genes that are exclusively expressed in one condition in the maqc-ii data set. gray points indicate genes with adjusted p-value value ≥ . the data point highlighted by the green elipse refers to the gene with high signal-to-noise ratio but low expression. the correlation is inferred using isotonic regression differential expression analysis of real data with unbalanced designs another real data set is taken from, consisting of highly dispersed cdna samples from females and males. de genes are inferred from male–female comparisons. following, a sensitivity analysis is predicted to find an over-representation of inferred de genes from the sex chromosomes. indeed, the top ten de genes always include genes from sex chromosomes. all methods except absseq and voom identify de genes beyond sex chromosomes. this may indicate that absseq and voom retain higher specificity than the remaining methods and that alternative methods may not well model variance introduced by unequal smaple sizes. ebseq produces the lowest number of de genes from sex chromosomes but the highest number from autosomes, confirming the previously observed lack of power of this method for the analysis of data with such high dispersion and uneven sample sizes. given the unequal sample size in this data set, the similar performance of absseq to that of alternative methods also suggests that our model is able to handle unequal sample sizes and high dispersion. in particular, in absseq, we attempted to compensate for unequal sample size by adding expected reads counts to the smaller group until sample sizes are equal. we always take the mean reads count of the small group as expected count, in order to minimize possible biases in subsequent variance estimations. this compensation step is likely crucial for unbalanced data designs, especially in case of even larger differences than in our test data set. in the future, it may be worth exploring in more detail alternative compensation procedures.table number of de genes from sex chromosomes detected by each method in the hapmap-ceu data set at fdr-ajusted p-value of moderating fold change fold change often serves as a more informative indicator for biologists to identify des. it is also utilized in gene ranking to select candidates for further investigation and visualization. however, the fold change neglects variance across samples and might not necessarily be informative, especially for genes with low counts. to overcome this problem, deseq introduces an empirical bayes shrinkage for fold change estimation, which moderates the log fold change according to gene-specific dispersion. fold change can also be represented as a function of absolute count differences, suggesting a potential moderation of fold change via counts difference. therefore, we introduce a fold change shrinkage procedure according to count differences and dispersion. figure  shows how it works using the bottomly data set. genes with small counts tend to have high raw fold changes, which constrains reliable gene ranking by fold change at dynamic expression level. shrinking fold change by adding pseudocounts according to expression level removes this trend.fig. moderation of log fold change. a raw data of the bottomfly study. b the same data corrected by expression level. c the same data corrected by expression level and gene-specific dispersion. de genes are shown in red. non-de genes with high log fold change are marked by green elipses however, this shrinkage approach neglects the gene-specific dispersion and thus shows no effects on non-de genes with high dispersion as well as high expression. after taking account of gene-specific dispersion, the fold changes approximately reflect de genes and produce nearly evenly distributed fold change values, apparently improving gene ranking and visualization. notably, shrinkage with gene-specific dispersion only influences a small part of genes with high dispersion. unlike the approach in deseq our shrinkage method on gene-specific dispersion is based on p-values and therefore does not change the number of inferred significant de genes. in practice, users can obtain all three types of fold change values in absseq.  here we introduce a new method for differiential expression analysis of rna-seq count data, absseq. distinct from other current methods, absseq infers de genes through the absolute differences in gene expression and assumes the differences to be influenced by two sources of variation: that found for average gene expression levels and that found for the magnitude of differential expression. our approach employs a nb distribution to model these two parts and, as a consequence, it is able to detect de genes more effectively than existing methods, as demonstrated by our analysis of both real and simulated data. in particular, absseq shows an advantage in discriminating de genes against non-de ones, it applies an efficient outlier detection approach and is thus robust against outliers. morever, absseq inferred p-values correlate with the magnitude of count differences, thus producing a linear relationship between snr ratio and p-value. as a result, it reduces type i error rates at both very low and high expression level and it also leads to a smaller number of highly significant de genes with small fold change. in addition, absseq introduces a procedure to shrink fold change according to the smoothed dispersion across expression level and observed dispersion, which permits fold change comparisons across genes and thus might favor downstream analyses, such as gene set enrichment analysis by ranking, clustering and visualization or candidate selection. a potential improvement of our approach in the future may be to adapt it to allow usage of more complex models which consider multiple conditions, its combination with additional normalization procedures, such as those implemented in peer and svaseq, which can further help to filter out unwanted variation, and also adjustment of our approach to allow for analysis of de at the transcript level. in summary, based on our analysis, we conclude that absseq represents a highly efficient approach for identification of significant des across a wide range of conditions and may help efficient downstram analysis of des. 
absseq specifically considers the overall magnitude of expression differences, which enhances the power in detecting truly differentially expressed genes by reducing false positives at both very low and high expression level. identification of significant differential expression represents a crucial initial step in these analyses, on which most subsequent inferences of biological functions are built. yet, for identification of these subsequently analysed genes, most studies use an additional minimal threshold of differential expression that is not captured by the applied statistical procedures. here we introduce a new analysis approach, absseq, which uses a negative binomal distribution to model absolute expression differences between conditions, taking into account variations across genes and samples as well as magnitude of differences. in comparison to alternative methods, absseq shows higher performance on controling type i error rate and at least a similar ability to correctly identify differentially expressed genes.
 the recent advances in next generation sequencing technology have made the sequencing of rna an extemely popular approach for gene expression analysis. identification of significant differential expression represents a crucial initial step in these analyses, on which most subsequent inferences of biological functions are built. yet, for identification of these subsequently analysed genes, most studies use an additional minimal threshold of differential expression that is not captured by the applied statistical procedures.  here we introduce a new analysis approach, absseq, which uses a negative binomal distribution to model absolute expression differences between conditions, taking into account variations across genes and samples as well as magnitude of differences. in comparison to alternative methods, absseq shows higher performance on controling type i error rate and at least a similar ability to correctly identify differentially expressed genes.  absseq specifically considers the overall magnitude of expression differences, which enhances the power in detecting truly differentially expressed genes by reducing false positives at both very low and high expression level. in addition, absseq offers to calculate shrinkage of fold change to facilitate gene ranking and effective outlier detection. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords rna-seqtranscriptome analysisdifferential gene expressionabsseqnegative binomial distributionhttp://dx.doi.org/ /501100001659deutsche forschungsgemeinschaft schu 1415/ schu 1415/9schulenburg hinrich http://dx.doi.org/ /501100001659deutsche forschungsgemeinschaft ro 2994/3http://dx.doi.org/ /501100002347bundesministerium für bildung und forschung deep bmbf tp http://dx.doi.org/ /501100002347bundesministerium für bildung und forschung sysinflame tp4issue-copyright-statement© the author 2016 