BACKGROUND
rna sequencing , which applies high-throughput dna sequencing technologies to directly sequence complementary dnas , has completely transformed the way in which transcriptomes are studied.
in particular, it permits the quantitative analysis of gene expression and transcript variant discovery, which was not made possible with the previous microarray technologies  <cit> .
rna-seq is increasingly being used to investigate a wide range of biological and medical questions, e.g., in genomics research  <cit>  and in clinic use  <cit> .
in rna-seq experiments, millions of short fragments  are sequenced from samples and aligned back to a reference genome.
the expression level of a feature  is then measured by the read count which is the number of short reads that map to the feature.
when rna-seq measurements are made for multiple samples from different biological conditions, a question of particular interest is to identify genes/features that are differentially expressed across conditions.
this is the primary aim of rna-seq differential expression analysis.
there have been a number of statistical approaches proposed for differential expression analysis of rna-seq data, and they broadly fall into two categories: parametric or nonparametric.
in  <cit> , the over-dispersed rna-seq data is transformed so that the poisson distribution can be used to model read counts.
edger  <cit> , deseq  <cit> , and sseq  <cit>  instead assume the negative binomial distribution on read counts--a flexible probability model allowing a larger variance than mean.
the differences among these three approaches lie mainly in their different ways to estimate the dispersion parameter.
ebseq  <cit>  and bayseq  <cit>  also assume the negative binomial distribution, but they were cast within an empirical bayesian framework.
all the above parametric approaches are generally very efficient when the distributional assumption holds.
however, violation of distributional assumptions or a poor estimation of parameters often leads to unreliable results.
noiseq  <cit>  is a data-adaptive nonparametric approach that uses both log fold changes and absolute expression differences as test statistics.
it is effective in controlling the false discovery rate  while being robust against the sequencing depth.
samseq  <cit>  is another nonparametric approach that utilizes a wilcoxon statistic.
it estimates the false discovery rate by a permutation plug-in procedure and thus is not sensitive to outliers in the data.
recently, an efficient algorithm based on a markov random field model, called mrfseq, was developed  <cit> .
different from previous methods, mrfseq takes advantage of the additional gene co-expression data to effectively alleviate the selection bias of differentially expressed genes against genes with low read counts.
for more discussions and comparisons of these differential expression analysis methods, we refer readers to  <cit>  and  <cit> .
in this paper, we propose a new data-driven nonparametric approach called lfcseq for differential expression analysis of rna-seq data.
basically, it is based on a similar principle to noiseq, but uses only log fold changes as the test statistic.
to conduct a statistical test for each gene, lfcseq estimates a null or noise probability distribution by contrasting log fold changes for a selected set of genes at similar expression levels.
in contrast, noiseq relies on a null distribution estimated from all genes within an experimental condition regardless of their expression levels.
however, as we shall demonstrate later, the null distribution of log fold changes varies considerably for genes at different expression levels, which makes the results from noiseq less reliable.
methods
notation
although biological experimental designs may vary greatly, rna-seq data generated for differential expression analysis can all be written into a matrix n, whose element nij is the number of reads mapped to gene i in sample j from an experimental condition a or b. without ambiguity, we also use a  to denote the set of samples j under the condition a .
that is, if j∈a, it indicates that sample j is under the experimental condition a rather than condition b. let xi be a binary random variable indicating whether gene i is differentially expressed between two conditions a and b. we have xi =  <dig> if gene i is differentially expressed  and xi =  <dig> if gene i is not differentially expressed .
typically, only a few samples are available in current rna-seq experimental data; however, there could instead have up to tens of thousands of genes under examination.
in the present study, we limit our discussions to two experimental conditions only, although our proposed approach can be extended to three or more conditions.
normalization
since different samples may have different sequencing depths, the read counts nij are not directly comparable across samples before being properly normalized  <cit> .
a simple normalization scheme is to divide the read counts by the sample library size and gene length  <cit> .
however, this total-count normalization was shown to be problematic, as the normalized read count of a gene is adversely affected by expression levels of all the other genes  <cit> .
many sophisticated normalization procedures have been proposed, including the trimmed means of m values  normalization in edger  <cit> , quantile normalization  <cit> , a 'median' normalization method in deseq  <cit>  and a goodness-of-fit method in poissonseq  <cit> .
in our experiments below, we use the goodness-of-fit method to normalize read counts.
it defines the sequencing depth for sample j as d ^j= ∑i∈snij/ ∑i∈s∑jnij, where s is a half set of genes that are least differentially expressed between two conditions as estimated by a poisson goodness-of-fit statistic  <cit> .
the normalized read count nij is subsequently computed as nij=nij/d ^j.
lfcseq
let n ¯ia and n ¯ib be the means of the normalized read counts for gene i of samples under conditions a and b, respectively.
that is, n ¯ia=1|a|Σj∈anij and n ¯ib=1|b|Σj∈bnij.
in lfcseq, we use the log fold change of mean read counts, i.e.,
li=log2n ¯ian ¯ib,
as the statistic to test differential expression.
because there are usually only a small number of samples under one condition, no read counts could be reliably identified as outliers.
therefore, we choose to use the mean instead of the median of read counts in the above definition .
however, when there is obvious evidence that the outliers of read count exist or when the number of samples is large enough, median may be a better choice than mean.
on the other hand, to avoid the division by zero, genes with zero read counts in all samples are removed from the analysis, and the zero counts are replaced by  <dig>  for the rest of genes, as in  <cit> .
we try to build a null or noise distribution for log fold changes by contrasting gene read counts within the same condition.
to this end, we first divide the samples within a same condition into two groups of almost equal size.
let a <dig> and a <dig> be the two resulting groups of samples under condition a such that a = a <dig> ∪ a <dig> and |a1|=|a| <dig>  as we did in the preceding paragraph, let n ¯ia1=1|a1|Σj∈a1nij and n ¯ia2=1|a2|Σj∈a2nij.
they are the means of the normalized read counts for gene i within each group of samples.
then, the log fold change of read counts between two groups a <dig> and a <dig> is computed as
lia1∪a2=log2n ¯ia1n ¯ia <dig>
when |a| ≤  <dig>  we may compute the log fold change value lia1∪a <dig> for all the possible partitions of a into a <dig> and a <dig>  however, when |a| > <dig>  we compute it only for  <dig> random partitions in order to reduce the computational cost.
finally, we pool all these log fold change values together, and denote the resulting collection by lia.
by applying the same procedure as above, we can obtain a collection lb of log fold change values of read counts for gene i within condition b.
given a gene i, we define its neighborhood as a set of genes with similar expression strength across conditions.
specifically, we define the neighborhood n of gene i as n={i′:|n ¯i′a∪b-n ¯ia∪b|<∈i}}, where n ¯ia∪b=1|a∪b|∑j∈a∪bnij and ∈i set to a value such that n  would contain a predefined number of genes .
then, we build a null fold change distribution li for gene i by using
li= ⋃i′∈nli′a∪li′b
note that this null distribution is gene-specific, as it takes into account only genes from the neighborhood of gene i. a special case of the above proposed approach is obtained when the neighborhood of a gene includes all the genes in a sample under investigation.
with the log fold change li of read counts of gene i between two conditions and a null fold change distribution li, we approximate the probability of gene i being not differentially expressed as the fraction of points from li that correspond to a larger absolute fold change value than |li|.
therefore, we may write this probability as
p=|{l:|l|>|li|,l∈li}||li|.
the above proposed approach to estimate the probability of a gene being non-de is motivated by a previous observation in  <cit>  that the squared coefficient of variation  decreases as gene expression levels increase.
we further found that the standard error of the null distribution li decreases considerably  as gene expression levels increase, as demonstrated in figure  <dig>  it clearly tells us that using a common null distribution to approximate the probability of genes being de or non-de, regardless of their expression levels, is not sufficient or appropriate.
therefore, we choose to group genes at similar expression levels and estimate the null fold change probability distribution based only on genes within the same group.
as shown in figure  <dig>  the estimated null distributions vary substantially across different groups of genes.
in general, the null distributions from genes of lower expression levels tend to shift to the right with heavier tails.
lfcseq is implemented in r and publicly available at http://www <dig> spms.ntu.edu.sg/~chenxin/lfcseq/.
relation to noiseq
lfcseq was developed based on a similar principle to the nonparametric approach noiseq  <cit> .
it is worth pointing out their major differences.
first, noiseq uses not only the log fold change li but also the absolute difference |di| of mean read counts as the statistics to test gene i for differential expression, where the absolute difference |di| is defined as |di|=|n ¯ia-n ¯ib|.
second, noiseq estimates the null joint probability distribution  by computing the log fold change r and absolute difference d for every pair of samples within a same condition  and for every gene .
consequently, a common null distribution is applied to all genes in noiseq to compute the probability of a gene being de.
that is,
pnoiseq=|{:|l|<|li|,|d|<|di|,∈}|||
RESULTS
datasets
we test the performance of lfcseq on two simulated and three real rna-seq datasets, and compare it with six existing parametric and nonparametric approaches, including noiseq, samseq, edger, deseq, sseq and ebseq .
simulation  <dig>  in this simulated dataset, there are a total of  <dig>  genes and their read counts are generated from a negative binomial distribution under each condition a or b,
nij~nb
where µij and σij <dig> are the mean and variance, respectively.
as in  <cit> , we further let µij = e{nij } = qia dj under condition a and µij = e{nij } = qib dj under condition b, where qia and qib represent the true expression values of gene i under condition a and b, respectively, and where dj represents the sequencing library size of sample j. for the variance, we let σij2=μij+ϕi.μij <dig>  where φi is the dispersion parameter of the negative binomial distribution.
as a typical setting, 30% of the genes are simulated to be differentially expressed, among which 70% are set to be up-regulated.
the library size factors are generated from the uniform distribution dj ~ u .
we consider three different sample sizes |a| = |b| =  <dig>   <dig> and  <dig> under each condition.
simulation  <dig>  we generate read counts for  <dig>  genes using the same procedure as above in simulation  <dig>  except that the parameter values of qia, qib , and φi are randomly sampled with replacement from the experimental bottomly's dataset  <cit> .
thus we expect this setting is more realistic than the previous one in simulation  <dig>
maqc dataset.
maqc dataset  <cit>  contains two rna sample types, stratagene's human universal reference rna  and ambion's human brain reference rna .
each sample type has seven replicates.
in this dataset,  <dig> genes have been assayed by the quantitative real-time polymerase chain reaction .
as in  <cit> , a gene is considered as differentially expressed if the log fold change ratio of its cycle threshold values exceeds  <dig> or as non-differentially expressed if this log fold change ratio is smaller than  <dig> .
as a result, we obtain  <dig> de genes and  <dig> non-de genes from the qrt-pcr gold-standard to assess the performance of the proposed approach.
griffith's dataset.
gene expression is compared between two human colorectal cell lines  <cit> , mip <dig> and mip/5-fu, of the fluorouracil -resistant and -nonresistant phenotype, respectively.
qrt-pcr measurements were made for  <dig> genes.
a two-tailed t-test was applied to identify de and non-de genes with a cutoff point  <dig> , which left  <dig> de genes and  <dig> non-de genes for performance evaluation.
sultan's dataset.
gene expression of two human cell lines, ramos b and hek 293t, were compared using rna-seq  <cit> .
in this dataset, there are two replicates for each cell line.
see additional file  <dig> for further details of these testing datasets.
evaluation criteria
we evaluate the performance of lfcseq from the following two aspects.
first, we evaluate its ability to discriminate between de genes and non-de genes by ranking genes in order of significance for differential expression between conditions.
with the gene ranking list, we plot a receiver operating characteristic  curve and compute the area under the curve  to measure the overall discriminating ability.
then, lfcseq is compared with six other approaches in terms of auc without imposing any arbitrary cutoffs.
for lfcseq, we rank genes in increasing order of the probability p .
for three parametric approaches that assumed the negative binomial distribution , genes are ranked according to their estimated nominal p-values.
for samseq, we use the false discovery rates  estimated by a permutation plug-in method and, for noiseq and ebseq, we use their estimated probabilities of genes being differentially expressed for ranking.
second, we evaluate the experimental results of lfcseq in terms of precision, sensitivity, and f-score.
these evaluation metrics are defined as follows: pre  = tp/, sen  =tp/, and fs  = <dig> × pre × sen/, where tp, fp, and fn are the number of true positives, the number of false positives and the number of false negatives, respectively.
note that the metric f-score is the harmonic mean of sensitivity and precision and thus measure the overall differential expression inference performance of a method.
in general, the higher the f-score, the better the inference performance.
in order to compute precision and sensitivity scores, all the approaches used their respective default settings to call a list of de genes.
specifically, lfcseq, noiseq, and ebseq used a probability cutoff of  <dig> ,  <dig> , and  <dig> , respectively.
samseq used a fdr cutoff of  <dig> , while edger, deseq, and sseq all used a p-value cutoff of  <dig>  after adjusted for multiple testing.
performance on simulated data
the experimental results of precision, sensitivity and f-scores are summarized in table  <dig> and in table s <dig> in additional file  <dig> for simulation  <dig> and  <dig>  respectively.
clearly, lfcseq has the best overall performance as it achieved the highest f-scores in all the tests.
noiseq achieved relatively higher precision scores than lfcseq.
however, its sensitivity scores are much lower than those of lfcseq so that its overall performance becomes inferior to lfcseq.
for example, in simulation  <dig>  the f-scores of lfcseq are  <dig> ,  <dig> , and  <dig>  for the tests with  <dig>   <dig>  and  <dig> replicates, respectively.
they are  <dig> ,  <dig> , and  <dig>  greater than the corresponding f-scores of noiseq.
samseq did not call any de genes in the experiments with  <dig> replicates.
this is not surprising considering that the power of the wilcoxon test is generally low with a few replicates.
however, it is not clear why sseq did not call any de genes either, while the other two similar parametric approaches edger and deseq performed relatively well in both precision and sensitivity.
the numbers of replicates per condition are  <dig>   <dig> and  <dig>  respectively.
the highest precision, sensitivity and f-scores achieved are highlighted in bold.
the corresponding results obtained in simulation  <dig> can be found in additional file  <dig>
in addition, we compared lfcseq with a simple hypergeometric test  when the numbers of replicates per condition are  <dig> and  <dig> in simulation  <dig>  in the simple hypergeometric test, the null distribution for gene i is built on the randomly permuted samples of gene i between conditions a and b, instead of using the neighborhood n .
from figure s <dig> and table s <dig> in additional file  <dig>  it can be seen that lfcseq performs better than shgt in the terms of both auc values and f-scores.
<dig>  performance on real data
on maqc dataset, the curves of precision, sensitivity and f-scores obtained with varying number of replicates, as well as the precision-sensitivity curves obtained with  <dig> replicates per condition, are shown in figure  <dig>  similar to the results in the previous simulation study, lfcseq achieved higher sensitivity and the comparable levels of precision with other methods.
as a result, it has the highest f-scores and hence the best overall performance in all the tests.
in comparison, noiseq provides higher precision than lfcseq, but its sensitivity scores are significantly lower than lfcseq's by up to 22%.
on the other hand, samseq achieved comparably high sensitivity scores with lfcseq, but its precision scores are always the lowest among all the tested methods.
it is interesting to note that samseq behaved differently in the simulation study, where its precision scores are instead higher than lfcseq's in most cases.
as the number of replicates increases, noiseq maintains a relatively stable level of precision while all other approaches lose some precision.
this result is in agreement with the observation in  <cit> .
in addition, the precision-sensitivity curves also clearly indicate that lfcseq is a high-performing approach for differential expression analysis of rna-seq data, as it yields the improved balance between precision and sensitivity.
on griffith's dataset, the corresponding curves of precision, sensitivity and fscores are presented in figures s4-s <dig> in additional file  <dig>  overall, we observed a similar pattern of performance to that observed on maqc dataset.
one noticeable difference is that although lfcseq still achieved the best overall performance in terms of f-score, there are only marginal improvements over the two parametric approaches edger and deseq.
recall that only  <dig> truly non-de genes were identified from griffith's limited rt-pcr data for the validation of prediction results.
such a small true negative dataset is hardly sufficient to fully characterize the performance behavior of a method.
on sultan's dataset, no gold-standard is available for performance validation.
instead of computing precision and sensitivity scores, we plotted in figure s <dig> in additional file  <dig> the fold changes of genes against their mean expression levels on the logarithmic scale.
in those scatter plots, each red dot represents a gene being called de while each black dot represents a gene being called non-de.
as we can see, lfcseq called de genes at both high and low expression ranges.
however, noiseq called few de genes at low expression ranges, which might suggest that noiseq is biased against genes with low read counts and that its sensitivity could still be very low as we observed earlier.
we also notice that sseq called a considerably less number of de genes than other approaches, which indicates that it is very conservative when calling de genes.
CONCLUSIONS
in this paper we proposed a new nonparametric approach for differential expression analysis of rna-seq data.
it relies on the statistical tests of log fold changes of gene read counts between and within biological conditions.
following the observation that the standard errors of log fold changes vary considerably with gene expression levels, we choose to create a gene-specific null probability distribution for each gene rather than a common null probability distribution for all genes.
this is done by considering the gene neighborhood, which is defined as a set of genes at similar expression levels.
as a result, the estimated probability of a gene being de depends only on the read counts of genes from its neighborhood.
our experimental results demonstrate that the proposed approach lfcseq outperforms its competitors in better ranking the truly de genes ahead of non-de genes.
it has the best overall performance as it achieved the highest f-scores in almost all our tests .
the improvements over other methods are especially remarkable when the number of replicates is small.
in such cases, those parametric methods based on negative binomial distribution, such as edger, deseq and sseq, could not estimate the distributional parameters accurately, while for the nonparametric samseq method, its wilcoxon statistic has a relatively low testing power.
in this study, we applied a pre-specified probability cutoff of  <dig>  for our approach lfcseq.
this cutoff generally works well, as shown in our experiments on both simulated data and real rna-seq data.
however, it is certainly of interest to develop a data-driven cutoff selection method for a wide applicability of the approach.
in addition, it is also interesting to formulate a framework to control the false discovery rate  <cit>  for our approach.
we will explore these in future work.
competing interests
the authors declare that they have no competing interests.
authors' contributions
bql and xc conceived the idea.
bql and xc contributed to the design of the study.
bql processed the data and conducted simulation and real dataset experiments.
bql, lfz and xc wrote the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary text and figures.
this file contains related codes to use existing approaches, information and results for simulated and real datasets.
click here for file
acknowledgements
the authors would like to thank ei-wen yang for his assistance in data analysis and prof tao jiang for his helpful comments.
this work was partially supported by the singapore national medical research council grant  and the ministry of education academic research fund .
declarations
the publication cost for this article was funded by a grant from the singapore ministry of education academic research fund .
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2014: proceedings of the 25th international conference on genome informatics : genomics.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/15/s <dig>

 rna sequencing, which applies high-throughput dna sequencing technologies to directly sequence complementary dnas, has completely transformed the way in which transcriptomes are studied. in particular, it permits the quantitative analysis of gene expression and transcript variant discovery, which was not made possible with the previous microarray technologies. rna-seq is increasingly being used to investigate a wide range of biological and medical questions, e.g., in genomics research and in clinic use. in rna-seq experiments, millions of short fragments are sequenced from samples and aligned back to a reference genome. the expression level of a feature is then measured by the read count which is the number of short reads that map to the feature. when rna-seq measurements are made for multiple samples from different biological conditions, a question of particular interest is to identify genes/features that are differentially expressed across conditions. this is the primary aim of rna-seq differential expression analysis. there have been a number of statistical approaches proposed for differential expression analysis of rna-seq data, and they broadly fall into two categories: parametric or nonparametric. in, the over-dispersed rna-seq data is transformed so that the poisson distribution can be used to model read counts. edger, deseq, and sseq instead assume the negative binomial distribution on read counts--a flexible probability model allowing a larger variance than mean. the differences among these three approaches lie mainly in their different ways to estimate the dispersion parameter. ebseq and bayseq also assume the negative binomial distribution, but they were cast within an empirical bayesian framework. all the above parametric approaches are generally very efficient when the distributional assumption holds. however, violation of distributional assumptions or a poor estimation of parameters often leads to unreliable results. noiseq is a data-adaptive nonparametric approach that uses both log fold changes and absolute expression differences as test statistics. it is effective in controlling the false discovery rate while being robust against the sequencing depth. samseq is another nonparametric approach that utilizes a wilcoxon statistic. it estimates the false discovery rate by a permutation plug-in procedure and thus is not sensitive to outliers in the data. recently, an efficient algorithm based on a markov random field model, called mrfseq, was developed. different from previous methods, mrfseq takes advantage of the additional gene co-expression data to effectively alleviate the selection bias of differentially expressed genes against genes with low read counts. for more discussions and comparisons of these differential expression analysis methods, we refer readers to and. in this paper, we propose a new data-driven nonparametric approach called lfcseq for differential expression analysis of rna-seq data. basically, it is based on a similar principle to noiseq, but uses only log fold changes as the test statistic. to conduct a statistical test for each gene, lfcseq estimates a null or noise probability distribution by contrasting log fold changes for a selected set of genes at similar expression levels. in contrast, noiseq relies on a null distribution estimated from all genes within an experimental condition regardless of their expression levels. however, as we shall demonstrate later, the null distribution of log fold changes varies considerably for genes at different expression levels, which makes the results from noiseq less reliable. methods notation although biological experimental designs may vary greatly, rna-seq data generated for differential expression analysis can all be written into a matrix n, whose element nij is the number of reads mapped to gene i in sample j from an experimental condition a or b. without ambiguity, we also use a to denote the set of samples j under the condition a. that is, if j∈a, it indicates that sample j is under the experimental condition a rather than condition b. let xi be a binary random variable indicating whether gene i is differentially expressed between two conditions a and b. we have xi = if gene i is differentially expressed and xi = if gene i is not differentially expressed. typically, only a few samples are available in current rna-seq experimental data; however, there could instead have up to tens of thousands of genes under examination. in the present study, we limit our discussions to two experimental conditions only, although our proposed approach can be extended to three or more conditions. normalization since different samples may have different sequencing depths, the read counts nij are not directly comparable across samples before being properly normalized. a simple normalization scheme is to divide the read counts by the sample library size and gene length. however, this total-count normalization was shown to be problematic, as the normalized read count of a gene is adversely affected by expression levels of all the other genes. many sophisticated normalization procedures have been proposed, including the trimmed means of m values normalization in edger, quantile normalization, a'median' normalization method in deseq and a goodness-of-fit method in poissonseq. in our experiments below, we use the goodness-of-fit method to normalize read counts. it defines the sequencing depth for sample j as d ^j= ∑i∈snij/ ∑i∈s∑jnij, where s is a half set of genes that are least differentially expressed between two conditions as estimated by a poisson goodness-of-fit statistic. the normalized read count nij is subsequently computed as nij=nij/d ^j. lfcseq let n ¯ia and n ¯ib be the means of the normalized read counts for gene i of samples under conditions a and b, respectively. that is, n ¯ia=1|a|Σj∈anij and n ¯ib=1|b|Σj∈bnij. in lfcseq, we use the log fold change of mean read counts, i.e., li=log2n ¯ian ¯ib, as the statistic to test differential expression. because there are usually only a small number of samples under one condition, no read counts could be reliably identified as outliers. therefore, we choose to use the mean instead of the median of read counts in the above definition. however, when there is obvious evidence that the outliers of read count exist or when the number of samples is large enough, median may be a better choice than mean. on the other hand, to avoid the division by zero, genes with zero read counts in all samples are removed from the analysis, and the zero counts are replaced by for the rest of genes, as in. we try to build a null or noise distribution for log fold changes by contrasting gene read counts within the same condition. to this end, we first divide the samples within a same condition into two groups of almost equal size. let a and a be the two resulting groups of samples under condition a such that a = a ∪ a and |a1|=|a| as we did in the preceding paragraph, let n ¯ia1=1|a1|Σj∈a1nij and n ¯ia2=1|a2|Σj∈a2nij. they are the means of the normalized read counts for gene i within each group of samples. then, the log fold change of read counts between two groups a and a is computed as lia1∪a2=log2n ¯ia1n ¯ia when |a| ≤ we may compute the log fold change value lia1∪a for all the possible partitions of a into a and a however, when |a| > we compute it only for random partitions in order to reduce the computational cost. finally, we pool all these log fold change values together, and denote the resulting collection by lia. by applying the same procedure as above, we can obtain a collection lb of log fold change values of read counts for gene i within condition b. given a gene i, we define its neighborhood as a set of genes with similar expression strength across conditions. specifically, we define the neighborhood n of gene i as n={i′:|n ¯i′a∪b-n ¯ia∪b|<∈i}}, where n ¯ia∪b=1|a∪b|∑j∈a∪bnij and ∈i set to a value such that n would contain a predefined number of genes. then, we build a null fold change distribution li for gene i by using li= ⋃i′∈nli′a∪li′b note that this null distribution is gene-specific, as it takes into account only genes from the neighborhood of gene i. a special case of the above proposed approach is obtained when the neighborhood of a gene includes all the genes in a sample under investigation. with the log fold change li of read counts of gene i between two conditions and a null fold change distribution li, we approximate the probability of gene i being not differentially expressed as the fraction of points from li that correspond to a larger absolute fold change value than |li|. therefore, we may write this probability as p=|{l:|l|>|li|,l∈li}||li|. the above proposed approach to estimate the probability of a gene being non-de is motivated by a previous observation in that the squared coefficient of variation decreases as gene expression levels increase. we further found that the standard error of the null distribution li decreases considerably as gene expression levels increase, as demonstrated in figure it clearly tells us that using a common null distribution to approximate the probability of genes being de or non-de, regardless of their expression levels, is not sufficient or appropriate. therefore, we choose to group genes at similar expression levels and estimate the null fold change probability distribution based only on genes within the same group. as shown in figure the estimated null distributions vary substantially across different groups of genes. in general, the null distributions from genes of lower expression levels tend to shift to the right with heavier tails. lfcseq is implemented in r and publicly available at http://www spms.ntu.edu.sg/~chenxin/lfcseq/. relation to noiseq lfcseq was developed based on a similar principle to the nonparametric approach noiseq. it is worth pointing out their major differences. first, noiseq uses not only the log fold change li but also the absolute difference |di| of mean read counts as the statistics to test gene i for differential expression, where the absolute difference |di| is defined as |di|=|n ¯ia-n ¯ib|. second, noiseq estimates the null joint probability distribution by computing the log fold change r and absolute difference d for every pair of samples within a same condition and for every gene. consequently, a common null distribution is applied to all genes in noiseq to compute the probability of a gene being de. that is, pnoiseq=|{:|l|<|li|,|d|<|di|,∈}|||  datasets we test the performance of lfcseq on two simulated and three real rna-seq datasets, and compare it with six existing parametric and nonparametric approaches, including noiseq, samseq, edger, deseq, sseq and ebseq. simulation in this simulated dataset, there are a total of genes and their read counts are generated from a negative binomial distribution under each condition a or b, nij~nb where µij and σij are the mean and variance, respectively. as in, we further let µij = e{nij } = qia dj under condition a and µij = e{nij } = qib dj under condition b, where qia and qib represent the true expression values of gene i under condition a and b, respectively, and where dj represents the sequencing library size of sample j. for the variance, we let σij2=μij+ϕi.μij where φi is the dispersion parameter of the negative binomial distribution. as a typical setting, 30% of the genes are simulated to be differentially expressed, among which 70% are set to be up-regulated. the library size factors are generated from the uniform distribution dj ~ u. we consider three different sample sizes |a| = |b| = and under each condition. simulation we generate read counts for genes using the same procedure as above in simulation except that the parameter values of qia, qib, and φi are randomly sampled with replacement from the experimental bottomly's dataset. thus we expect this setting is more realistic than the previous one in simulation maqc dataset. maqc dataset contains two rna sample types, stratagene's human universal reference rna and ambion's human brain reference rna. each sample type has seven replicates. in this dataset, genes have been assayed by the quantitative real-time polymerase chain reaction. as in, a gene is considered as differentially expressed if the log fold change ratio of its cycle threshold values exceeds or as non-differentially expressed if this log fold change ratio is smaller than. as a result, we obtain de genes and non-de genes from the qrt-pcr gold-standard to assess the performance of the proposed approach. griffith's dataset. gene expression is compared between two human colorectal cell lines, mip and mip/5-fu, of the fluorouracil -resistant and -nonresistant phenotype, respectively. qrt-pcr measurements were made for genes. a two-tailed t-test was applied to identify de and non-de genes with a cutoff point, which left de genes and non-de genes for performance evaluation. sultan's dataset. gene expression of two human cell lines, ramos b and hek 293t, were compared using rna-seq. in this dataset, there are two replicates for each cell line. see additional file for further details of these testing datasets. evaluation criteria we evaluate the performance of lfcseq from the following two aspects. first, we evaluate its ability to discriminate between de genes and non-de genes by ranking genes in order of significance for differential expression between conditions. with the gene ranking list, we plot a receiver operating characteristic curve and compute the area under the curve to measure the overall discriminating ability. then, lfcseq is compared with six other approaches in terms of auc without imposing any arbitrary cutoffs. for lfcseq, we rank genes in increasing order of the probability p. for three parametric approaches that assumed the negative binomial distribution, genes are ranked according to their estimated nominal p-values. for samseq, we use the false discovery rates estimated by a permutation plug-in method and, for noiseq and ebseq, we use their estimated probabilities of genes being differentially expressed for ranking. second, we evaluate the experimental results of lfcseq in terms of precision, sensitivity, and f-score. these evaluation metrics are defined as follows: pre = tp/, sen =tp/, and fs = × pre × sen/, where tp, fp, and fn are the number of true positives, the number of false positives and the number of false negatives, respectively. note that the metric f-score is the harmonic mean of sensitivity and precision and thus measure the overall differential expression inference performance of a method. in general, the higher the f-score, the better the inference performance. in order to compute precision and sensitivity scores, all the approaches used their respective default settings to call a list of de genes. specifically, lfcseq, noiseq, and ebseq used a probability cutoff of,, and, respectively. samseq used a fdr cutoff of, while edger, deseq, and sseq all used a p-value cutoff of after adjusted for multiple testing. performance on simulated data the experimental results of precision, sensitivity and f-scores are summarized in table and in table s in additional file for simulation and respectively. clearly, lfcseq has the best overall performance as it achieved the highest f-scores in all the tests. noiseq achieved relatively higher precision scores than lfcseq. however, its sensitivity scores are much lower than those of lfcseq so that its overall performance becomes inferior to lfcseq. for example, in simulation the f-scores of lfcseq are,, and for the tests with and replicates, respectively. they are,, and greater than the corresponding f-scores of noiseq. samseq did not call any de genes in the experiments with replicates. this is not surprising considering that the power of the wilcoxon test is generally low with a few replicates. however, it is not clear why sseq did not call any de genes either, while the other two similar parametric approaches edger and deseq performed relatively well in both precision and sensitivity. the numbers of replicates per condition are and respectively. the highest precision, sensitivity and f-scores achieved are highlighted in bold. the corresponding results obtained in simulation can be found in additional file in addition, we compared lfcseq with a simple hypergeometric test when the numbers of replicates per condition are and in simulation in the simple hypergeometric test, the null distribution for gene i is built on the randomly permuted samples of gene i between conditions a and b, instead of using the neighborhood n. from figure s and table s in additional file it can be seen that lfcseq performs better than shgt in the terms of both auc values and f-scores. performance on real data on maqc dataset, the curves of precision, sensitivity and f-scores obtained with varying number of replicates, as well as the precision-sensitivity curves obtained with replicates per condition, are shown in figure similar to the results in the previous simulation study, lfcseq achieved higher sensitivity and the comparable levels of precision with other methods. as a result, it has the highest f-scores and hence the best overall performance in all the tests. in comparison, noiseq provides higher precision than lfcseq, but its sensitivity scores are significantly lower than lfcseq's by up to 22%. on the other hand, samseq achieved comparably high sensitivity scores with lfcseq, but its precision scores are always the lowest among all the tested methods. it is interesting to note that samseq behaved differently in the simulation study, where its precision scores are instead higher than lfcseq's in most cases. as the number of replicates increases, noiseq maintains a relatively stable level of precision while all other approaches lose some precision. this result is in agreement with the observation in. in addition, the precision-sensitivity curves also clearly indicate that lfcseq is a high-performing approach for differential expression analysis of rna-seq data, as it yields the improved balance between precision and sensitivity. on griffith's dataset, the corresponding curves of precision, sensitivity and fscores are presented in figures s4-s in additional file overall, we observed a similar pattern of performance to that observed on maqc dataset. one noticeable difference is that although lfcseq still achieved the best overall performance in terms of f-score, there are only marginal improvements over the two parametric approaches edger and deseq. recall that only truly non-de genes were identified from griffith's limited rt-pcr data for the validation of prediction results. such a small true negative dataset is hardly sufficient to fully characterize the performance behavior of a method. on sultan's dataset, no gold-standard is available for performance validation. instead of computing precision and sensitivity scores, we plotted in figure s in additional file the fold changes of genes against their mean expression levels on the logarithmic scale. in those scatter plots, each red dot represents a gene being called de while each black dot represents a gene being called non-de. as we can see, lfcseq called de genes at both high and low expression ranges. however, noiseq called few de genes at low expression ranges, which might suggest that noiseq is biased against genes with low read counts and that its sensitivity could still be very low as we observed earlier. we also notice that sseq called a considerably less number of de genes than other approaches, which indicates that it is very conservative when calling de genes.  in this paper we proposed a new nonparametric approach for differential expression analysis of rna-seq data. it relies on the statistical tests of log fold changes of gene read counts between and within biological conditions. following the observation that the standard errors of log fold changes vary considerably with gene expression levels, we choose to create a gene-specific null probability distribution for each gene rather than a common null probability distribution for all genes. this is done by considering the gene neighborhood, which is defined as a set of genes at similar expression levels. as a result, the estimated probability of a gene being de depends only on the read counts of genes from its neighborhood. our experimental results demonstrate that the proposed approach lfcseq outperforms its competitors in better ranking the truly de genes ahead of non-de genes. it has the best overall performance as it achieved the highest f-scores in almost all our tests. the improvements over other methods are especially remarkable when the number of replicates is small. in such cases, those parametric methods based on negative binomial distribution, such as edger, deseq and sseq, could not estimate the distributional parameters accurately, while for the nonparametric samseq method, its wilcoxon statistic has a relatively low testing power. in this study, we applied a pre-specified probability cutoff of for our approach lfcseq. this cutoff generally works well, as shown in our experiments on both simulated data and real rna-seq data. however, it is certainly of interest to develop a data-driven cutoff selection method for a wide applicability of the approach. in addition, it is also interesting to formulate a framework to control the false discovery rate for our approach. we will explore these in future work. competing interests the authors declare that they have no competing interests. authors' contributions bql and xc conceived the idea. bql and xc contributed to the design of the study. bql processed the data and conducted simulation and real dataset experiments. bql, lfz and xc wrote the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 supplementary text and figures. this file contains related codes to use existing approaches, information and results for simulated and real datasets. click here for file acknowledgements the authors would like to thank ei-wen yang for his assistance in data analysis and prof tao jiang for his helpful comments. this work was partially supported by the singapore national medical research council grant and the ministry of education academic research fund. declarations the publication cost for this article was funded by a grant from the singapore ministry of education academic research fund. this article has been published as part of bmc genomics volume supplement 2014: proceedings of the 25th international conference on genome informatics: genomics. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/15/s 
through extensive simulation study and rna-seq real data analysis, we demonstrate that the proposed approach could well rank the differentially expressed genes ahead of non-differentially expressed genes, thereby achieving a much improved overall performance for differential expression analysis. however, violation of distributional assumptions or a poor estimation of parameters often leads to unreliable results. in comparative experiments, differential expression analysis is commonly performed on rna-seq data to identify genes/features that are differentially expressed between biological conditions. with the advances in high-throughput dna sequencing technologies, rna-seq has rapidly emerged as a powerful tool for the quantitative analysis of gene expression and transcript variant discovery. to test each gene for differential expression, lfcseq estimates a null probability distribution of count changes from a selected set of genes with similar expression strength. in contrast, the nonparametric noiseq approach relies on a null distribution estimated from all genes within an experimental condition regardless of their expression levels. most existing statistical methods for differential expression analysis are parametric and assume either poisson distribution or negative binomial distribution on gene read counts. in this paper, we introduce a new nonparametric approach called lfcseq that uses log fold changes as a differential expression test statistic. differential expressionnonparametricrna-seq15- december the 25th international conference on genome informatics tokyo, japan
 with the advances in high-throughput dna sequencing technologies, rna-seq has rapidly emerged as a powerful tool for the quantitative analysis of gene expression and transcript variant discovery. in comparative experiments, differential expression analysis is commonly performed on rna-seq data to identify genes/features that are differentially expressed between biological conditions. most existing statistical methods for differential expression analysis are parametric and assume either poisson distribution or negative binomial distribution on gene read counts. however, violation of distributional assumptions or a poor estimation of parameters often leads to unreliable results.  in this paper, we introduce a new nonparametric approach called lfcseq that uses log fold changes as a differential expression test statistic. to test each gene for differential expression, lfcseq estimates a null probability distribution of count changes from a selected set of genes with similar expression strength. in contrast, the nonparametric noiseq approach relies on a null distribution estimated from all genes within an experimental condition regardless of their expression levels.  through extensive simulation study and rna-seq real data analysis, we demonstrate that the proposed approach could well rank the differentially expressed genes ahead of non-differentially expressed genes, thereby achieving a much improved overall performance for differential expression analysis. differential expressionnonparametricrna-seq15- december the 25th international conference on genome informatics tokyo, japan 