an advantage of rna - seq over other gene expression profiling technologies is that it allows a comprehensive assay that does not require probes for targets to be specified in advance .
it has particularly been used for de novo detection of splice junctions and allows genome wide expression profiling of organisms with unknown genome sequence . by obtaining millions of short
reads from the population of interest and by mapping these reads to the reference genome , rna - seq produces read count data . with enough reads from a sample ,
it has the potential to detect and quantify biologically significant rnas with low and moderate abundances . before detecting biologically significant rnas , systematic technical variations due to experimental variability
various procedures for normalization of rna - seq have been proposed in literature , such as transcripts parts per million , trimmed mean of m values , and quantile normalization .
though these methods have been frequently used , no comparative analysis has been presented so far .
previous methods for identification of differential expressed genes include bloom et al .   who identified differential expression by taking log ratio of the transcript counts ; hoen et al .  
used a student 's   t - test and alternatively also applied a bayesian model of vncio et al . .
suggested to use poisson model ( and fisher 's exact test , or a likelihood ratio test as an approximation to it ) to test for differential expression .
recently published methods , edger   and deseq   use a negative binomial distribution to test for differential expression as it allows for over dispersion .
we also propose two statistical methods for inferring differential expression for rna - seq data .
methods for normalization , differential expression , along with the details of the dataset used to test the performance of our methods are detailed in the next section .
results along with a systematic comparison are presented on three real datasets and we conclude with a brief discussion .
conducted rna - seq experiment with liver and kidney of a single human male using illumina genome analyzer sequencing platform .
each tissue was sequenced in seven lanes , split across two runs of the machine and two different cdna concentrations ( 1.5  pm , 3  pm ) . for this work ,
we only use data sequenced at 3  pm concentration ( five lanes for each sample ) and 17708 ensembl transcripts that mapped with the array probes . 
 
marioni et al .   conducted rna - seq experiment with liver and kidney of a single human male using illumina genome analyzer sequencing platform .
each tissue was sequenced in seven lanes , split across two runs of the machine and two different cdna concentrations ( 1.5  pm , 3  pm ) . for this work ,
we only use data sequenced at 3  pm concentration ( five lanes for each sample ) and 17708 ensembl transcripts that mapped with the array probes . 

profiled mirna expression from the normal peripheral blood mononuclear cells from two different individuals and cancer cells of myeloid lineage , k562 ( chronic myelocytic leukemia ) and hl60 ( acute promyelocytic leukemia ) using solexa technology . 
 
profiled mirna expression from the normal peripheral blood mononuclear cells from two different individuals and cancer cells of myeloid lineage , k562 ( chronic myelocytic leukemia ) and hl60 ( acute promyelocytic leukemia ) using solexa technology . 

analyzed 6 globin reduced with 6 nonreduced human whole blood rna samples using a tag sequencing method on the illumina high - throughput sequencing platform . 
 
mastrokolias et al .   analyzed 6 globin reduced with 6 nonreduced human whole blood rna samples using a tag sequencing method on the illumina high - throughput sequencing platform .
  normalization is a procedure to remove nonbiological influence on biological data and to make data comparable across experiments , runs , and lanes .
various normalization procedures have been proposed in literature for rna - seq and here we evaluate three different normalization methods : ( 1 ) transcripts parts per million , ( 2 ) trimmed mean of m values , ( 3 ) quantile normalization . at present , transcripts parts per million ( tpm ) is a standard procedure to normalize rna - seq data . using this method , number of reads of a transcript / sequence
resulting normalized data is reported as reads ( or transcripts ) per million for each sample .
one of the major problems with rna - seq data is that while the total number of reads for a sample is known , the composition of the rna population is unknown .
thus , tpm normalization method has its limitations for datasets with marked different rna composition .
trimmed mean of m values ( tmm ) normalization has been suggested to remove rna compositional bias as tmm equates the overall expression levels of genes between samples by estimation of relative rna production levels or scale factors . another method in use is quantile normalization which has previously been applied for microarrays . in quantile normalization
, the distribution of read counts in each lane is matched to a reference distribution defined in terms of median counts across sorted lanes .
  we propose two methods for inferring differential expression across two biological conditions with technical replicates , each of which yields one test statistics per gene : ( i ) likelihood ratio method ( lrm ) ( casella and berger ) , ( ii ) bayesian method ( bm ) , an extension of technique due to audic and claverie   for more than 2 replicates within a condition .
let xj denotes the observed number of reads mapped to a gene in replicate j(j = 1,2,  m ) under condition-1 and let yj denotes the observed number of reads mapped to a gene in replicate j(j = 1,2,  n ) for condition-2 .
since the number of reads mapped to a gene represents a small ( less than 5% ) fraction of the total number of reads obtained after sequencing , we assume xj and yj to follow independent poisson distribution with different parameters .
methods are detailed for a gene and the same need to be applied for all genes .    for condition-1 ,
j=1,2,  ,m , 
 where j denotes the true expression level of gene in replicate j. as xj 's occur independently , the likelihood function of x1 , x2,  xm is given by 
 ( 2)l = l(1,2, 
,xm)=j=1mejjxjxj!. 
 to identify genes with similar read count across replicates , we test the null hypothesis h0 : 1 = 2 =  = m = ( say ,  ) against the alternative h1 : i  j for some i  j. under h0 , the maximum likelihood estimate ( mle ) of  is given by 
 ( 3)^=j=1mxjm = xm , 
 where x = j=1xj and under h1 , the mle of j is given by 
 ( 4)^j = xj ,  j=1,2,  ,m . 
 the likelihood ratio for testing 1 = 2 =  = m =
( say ,  ) for condition-1 is given by 
 ( 5)1=sup  h0  lsup  h1  l=(x / m)xj=1m(xj)xj=(x)xmxj=1m(xj)xj . 

similarly , for condition-2 , yj follows poisson distribution with parameters j , j = 1,2 ,  , n. as derived above , the likelihood ratio for testing
1 = 2 =  , n = ( say ,  ) for condition-2 is given by 
 ( 6)2=(y)ynyj=1n(yj)yj , 
 where y = j=1yj .
for identifying differentially expressed genes across the two conditions , for a gene , define x = j=1xj and y = j=1yj to be independent poisson random variables with parameters m and n , respectively , and test if . the joint likelihood of the two conditions is given as 
 ( 7)l = l(,    x1,x2,  ,xm;y1,y2,  ,yn)=j=1mexjxj!j=1neyjyj ! , 
 and the unconditional mle 's of  and  are given by x / m and y /
n , respectively , mle of  under the hypothesis  =  is ( x + y)/(m + n ) .
the likelihood ratio for testing  =  is given by 
 ( 8)3=(mm+n)x(nm+n)y(x+y)x+yxxyy . 

we reject the null hypothesis for the small values of the statistic , 3 .    back in 1997 ,
the method of audic and claverie was used to establish the probability distribution governing the occurrence of the same rare event in repeated experiments and was applied for the analysis of digital gene expression profiles .
it was then described for only 2 replicates which we have attempted to extend to 3 or more replicates and apply to rna - seq data . as defined before
, x1 represents the number of reads mapped to a gene in replicate 1 of the condition-1 and follows poisson distribution 
 ( 9)p(x1)=ex1x1 ! , 
 where  denotes the actual number of reads mapped to the gene .
let x2 represents the number of reads mapped to a gene in replicate 2 of the condition-1 .
then , 
 ( 10)p(x2    x1)=0p(d=    x1)p(x2    d=)d , 
 where p(d = x1 ) in above equation is the posterior probability of  given x1 occurrences of a gene in an experiment and p(x2d =  ) = e/x2 ! is the probability of drawing x2 observations from poisson distribution with parameter . using bayes theorem , vncio et al
.   showed that , 
 ( 11)p(x2    x1)=(x1+x2)!x1!x2!2(x1+x2 + 1 ) , 
 where the prior distribution p(d =  ) is taken as uniform distribution over the interval [ 0 ,  ] .
we extended the above results when the condition is replicated thrice and 
 ( 12)p(x3    x1,x2)=0p(d=    x1,x2)p(x3    d=)d. 
 from bayes theorem , 
 ( 13)p(d=    x1,x2)=p(x1,x2    d=)p(d=)0p(x1,x2    d=)p(d=)d. 
 again , using uniform prior for  , we get 
 ( 14)p(d=    x1,x2)=2x1+x2 + 1e2x1+x2(x1+x2 ) ! , 
 which is a gamma random variable with scale parameter 2. this gives 
 ( 15)p(x3    x1,x2)=2x1+x2 + 1(x1+x2+x3)!(x1+x2)!x3!3x1+x2+x3 + 1 . 

therefore , 
 ( 16)p(x3    x1,x2)p(x2    x1)=(x1+x2+x3)!x1!x2!x3!3x1+x2+x3 + 1 . 
 similarly , if the condition is replicated m times , we consider the following probability .
( 17)p~(x1,x2,  ,xp,  ,xm )  = p(xm    x1,x2,  ,xp, 
x1 )  = ( x1+x2++xp++xm)!x1!x2!xm!mx1+x2++xm+1 . 
 in order to find genes with similar read counts within a condition
, we find two numbers a , b such that 
 ( 18)xp=0ap~(x1,x2,  ,xp,  ,xm)=,xp = bp~(x1,x2,  ,xp,  ,xm)=. 
 equation ( 18 ) implies that if the observation xm of the mth replicate lies in the interval [ a , b ] then we conclude with probability ( 1  2 ) that there are no systematic differences between the replicates .
similarly , the results can be derived for n replicates of a gene in condition-2 ( i.e. , yj , j = 1,2 ,  , n ) . for identifying differential expression across two conditions ,
define x1 = j=1xj , y1 = j=1yj to be independent poisson random variables with parameters m and n , respectively , and use ( 11 ) . under the bayesian method
, we can only identify genes that are different across two conditions if the number of replicates for the two conditions are the same ( i.e. , m = n ) .
  we assessed the variability within technical replicates using dataset 1 which comprises of liver and kidney tissue , each with five technical replicates and 17708 ensembl transcripts .
boxplots of unnormalized data from both liver and kidney samples are shown in figure 1(a ) .
we evaluate this variability statistically using a likelihood ratio method detailed in the previous section .
the analysis was performed at 1% , 2.5% , 5% , and 10% levels while considering two , three , four , and five replicates on the un - normalized data from kidney . as shown in table 1
, there is a decrease in the percentage of genes with similar counts as the number of replicates increases , which is expected ; however , the decreases is only marginal .
the percentage of genes with similar counts also decrease with the increase in the levels .
thus , dataset 1 is highly reproducible with few systematic differences among the replicates .  
we assess the impact of all three normalization methods using the likelihood ratio method at 1% , 2.5% , 5% , and 10% levels .
we used data from liver tissue with five replicates without normalization , with tmm , quantile , and tpm normalization .
it can be seen from table 2 that the percentage of genes with similar counts increased after tmm and quantile normalization and , thus , reduction in variability after normalization .
a gain of 2% is achieved after tmm or quantile normalization while the performance of tpm normalization was found to be poor .
figures 1 , 2 , and 3 represents boxplots of un - normalized , normalized after tpm , tmm and quantile for datasets 1 , 2 , and 3 , respectively .  
we compared the two proposed methods for inferring differentially expressed ( de ) genes : likelihood ratio method and bayesian method on datasets 2 and 3 .
for comparison between any two biological conditions , the read count values from the conditions can be categorized under three categories .
( 2 ) when one sample has zero or low counts and a reasonable count in the other .
this is an interesting biological phenomena where a gene is not expressed in one of the conditions .
for the quantile normalized normal versus hl60 data ( dataset 2 ) , 19 mirnas are absent in either of the two samples and present with a reasonable count for the other and 155 mirnas were present with read count of at least 5 in both the samples . using the likelihood ratio method at 1% level of significance , all 19 mirnas absent in either of the two conditions were identified as de and out of the 155 mirnas , 57 were identified as de . using the bayesian method at 1% level of significance , mirnas absent in either of the two conditions were also identified as de and out of the 155 mirnas , 58 were identified as de .
nearly same mirnas , except one , were identified as de using both the methods .
we also analyzed this dataset using deseq and edger and they did not identify mirnas absent in one of the two conditions .
of the 155 mirnas , deseq identified 3 mirnas as de with p value 0.01 and edger identified 4 mirnas as de with p value 0.01 .
see additional file 1 , 2 , and 3 in supplementary material available online at doi:10.5402/2012/817508 for detailed analysis and table 3 for a systematic comparison between methods for all three datasets . from additional
file 1 in supplementary material , it is clear that likelihood ratio method and bayesian method give very similar results for normal versus hl60 and normal versus k562 datasets ( dataset 2 ) .
both methods identified all mirnas previously identified as differentially expressed in vaz et al . .
however , deseq and edger could not identify most of the de mirnas reported in vaz et al . .
few mirnas experimentally verified using rnase protection assay ( rpa ) and real - time rt - pcr in vaz et al .  
in addition , we also identified differential expression of mir-181a family of hl60 , previously reported in .
for globin reduced versus non - reduced data ( dataset 3 ) , likelihood method reports 2513 significant genes at 1% level of significance , bayesian method reports 2344 at 1% level of significance , deseq reports 1505 with p value 0.01 and edger reports 2987 genes with p value 0.01 . from these numbers alone , it is difficult to comment on the performance of any method .  
figure 4 demonstrates the distribution of expression strength of the significant gene list obtained from likelihood ratio method , deseq , edger , and all genes .
one would expect the distribution of the significant gene lists to roughly follow the expression strength distribution for all genes . for likelihood ratio method and deseq ,
edger seems to be identifying genes from all expression strengths and thus not reflection biolog but the rigidity of its error models .
few genes experimentally verified in mastrokolias et al .   using qpcr ( i.e. , cxorf25 , hba1 , hba2 , hbd , hbb ) were obtained with high fold values in our analysis .
table 4 shows how a confidence interval was evaluated in bayesian method for quantile normalized normal - hl60 data .
hsa - let-7 g has a read count of 15117 in normal ( condition-1 ) and 6236 in hl60 ( condition-2 ) . using ( 18 ) , for one replicate , we estimated the lower and upper bound of the confidence interval around normal as 1386 and 1644 .
read count of 6236 for hsa - let-7 g in hl60 lies well outside the estimated confidence interval ( 1386 , 1644 ) .
thus , the read count in normal and hl60 are significantly different and reported in table 3 as t(i.e .
  we assessed three different types of normalizations and showed that though illumina data is highly replicable before normalization , normalization further reduces the technical variability , likelihood ratio method was used to statistically evaluate variation across replicates .
we also presented two methods for finding differentially expressed genes for rna - seq data with or without replicates , likelihood ratio method is a general method that does not impose any restriction on the equality of the number of replicates across the two conditions .
bayesian method on the other hand can only be applied if there is equality on the number of replicates for the two conditions being compared .
the performance of both the methods was compared to deseq , edger . for small rna dataset ,
likelihood ratio method and bayesian method perform similarly but better than edger and deseq . for dataset 3
, the distribution of the significant gene lists from likelihood ratio method and deseq roughly follows the expression strength distribution for all genes .
however , this was not true for edger . for both likelihood ratio method and bayesian method
, we assume that the underlying distribution for observed number of reads to be poisson .
poisson distribution is intuitively appealing and mathematically easy to handle but with a limitation that the mean and variance of poisson random variable are the same . to avoid this
, the efficiency of the proposed methods in identifying differentially expressed genes , their mathematical convenience , and simplicity should make these methods extremely useful .
 rna - seq is increasingly being used for gene expression profiling . in this approach , next - generation sequencing ( ngs ) platforms 
 are used for sequencing . due to highly parallel nature 
 , millions of reads are generated in a short time and at low cost . 
 therefore analysis of the data is a major challenge and development of statistical and computational methods is essential for drawing meaningful conclusions from this huge data . in here 
 , we assessed three different types of normalization ( transcript parts per million , trimmed mean of m values , quantile normalization ) and evaluated if normalized data reduces technical variability across replicates . 
 in addition , we also proposed two novel methods for detecting differentially expressed genes between two biological conditions : ( i ) likelihood ratio method , and ( ii ) bayesian method . 
 our proposed methods for finding differentially expressed genes were tested on three real datasets . 
 our methods performed at least as well as , and often better than , the existing methods for analysis of differential expression . 
1. Introduction
2. Material and Methods 
3. Normalization
4. Differential Expression
5. Results
6. Discussions