1. edgeR
The R Bioconductor package, edgeR,[12] provides statistical routines for detecting DE in RNA-seq data. The package is extremely flexible and can handle the count data irrespective of whether or not they are over-dispersed. If the data are over-dispersed, the NB model is used. Conversely, the Poisson model is used when there is no over-dispersion detected in the data. edgeR requires the data to be in either one of two formats: a single file containing a table of counts, with the first column containing the read (refer to 'tag' in the package) identifiers and the remaining columns containing the tag counts for each sample sequenced; or an individual file for each library, each with the first column for tag identifier and second column for counts.
Normalisation
The quantile-adjusted method is used to standardise total read counts (library sizes) across samples [11]. Samples are assumed to be independent and identically distributed from NB distribution (M*p, ϕ) (see details in the Model section), where ϕ is initially estimated from all the samples, M* is the geometric mean of original library sizes and p, the proportion of tag g in the sequenced sample, can then be estimated, providing values for M and ϕ. Linear interpolation of quantile function is used to equate the quantiles across samples. The process is updated with new ϕ and p until ϕ converges.
Model
edgeR is based on NB distribution. Let Y
gij
denote the observed data; where g is the gene (tag, exon, etc.), i is the experimental group and j is the index of samples. The counts can be modelled as Y
gij
~ NB(M
jpgi
, ϕ
g
), with mean μ
gi
= M
jpgi
and variance = , where M
j
represents the library size, (ie the sum of the counts of tags in a sample) and p
gi
represents the proportion of tag g of the sequenced sample for group i. ϕ
g
is the over-dispersion parameter (relative to Poisson) for accounting for biological or sample-to-sample variation. There are two options for ϕ
g
in the package; one is to use a common dispersion for all the tags and the other is to assume tagwise dispersion [3, 11]. For many applications, using a common dispersion will be adequate. For the tagwise dispersion, edgeR moderates the estimates towards a common dispersion. Moderation is determined using an empirical Bayes rule [3]. It is noted in general that tagwise dispersion penalises tags with great variability within groups. If the common dispersion estimate is much greater than 0, it indicates that there is more variability in the data than the Poisson model can account for, and NB distribution should be used. With ϕ
g
= 0, the NB distribution reduces to Poisson distribution.
Testing
edgeR employs an exact test for the NB distribution based on the normalised data. The test parallels with FET. The 'exactTest' function allows pairwise comparisons of groups. One of the objects produced by the function includes logFC, the log-fold change difference in the counts between the groups, and exact p-values. The results of the NB test can be accessed using the 'topTags' function, in which the adjusted p-values for multiple testing are reported using Benjamini and Hochberg's approach [13] as the default method of adjustment. Users can also supply their own desired adjustment method.
2. DEGseq
DEGseq is another R package specifically designed to identify DE from RNA-seq data [9]. The package includes two novel methods, the MA plot-based method (where M is the log ratio of the counts between two experimental conditions for gene g, and A is the two group average of the log concentrations of the gene) with a random sampling model (MARS) and the MA plot-based method with technical replicates (see details in the Models section), along with three existing methods: FET, the likelihood ratio test (LRT) and samWrapper. samWrapper was developed previously for microarray data analysis [14]. In this paper we focus our attention on MARS, FET and the LRT. Unlike edgeR, which allows over-dispersion, DEGseq assumes a binomial or Poisson distribution (which limits its application to data with no over-dispersion) and is extremely easy to use. The user needs only to specify the model, the normalisation methods and the input data, and the results will be saved to the user's designed folder. The input of the package is uniquely mapped reads (or tags), gene annotation of the corresponding genome and gene expression counts for each sample. The output includes a text file and a summary. The text file contains the original sample counts, p-value, and two q-values, indicating the expression difference between the two treatment groups for each gene, which are the adjusted p-values.
Normalisation
There are several choices for normalisation: 'none', 'loess' and 'median'. The recommended (or default) method is 'none.'
Models
(i) MARS
In the MARS model, RNA sequencing is modelled as a random sampling process [15]. Each read is sampled independently and uniformly from every possible nucleotide in the sample. The number of tags coming from a gene follows a binomial distribution, which can be approximated by a Poisson distribution. The model identifies and visualises DE genes based on the MA plot. It follows that M and A are both normally distributed, given that the samples from the two conditions are independent. The conditional distribution of M, given that A = a, is also normally distributed. Under the null hypothesis that the probabilities of the tag coming from a specific gene are the same between the two experimental conditions, a Z-score statistic for the gene can be calculated, and the p-value can be converted to indicate if the gene g is differentially expressed. The MA plot has been widely used to detect and visualise the intensity-dependent ratio of microarray data [16].
(ii) LRT
LRT was used by Marioni [1] to identify differentially expressed genes from sequencing data. In the dataset used, samples from each group are technical replicates. Let Y
gij
represent the number of reads mapped to gene g for the j-th lane (sample) from group i, as described earlier for the edgeR model. Y
gij
then can be modelled as a Poisson random variable with mean μ
gi
= M
jpgi
. Under the null hypothesis, the two groups A and B have the same value for gene g, p
gi
= p
j
and, under the alternative hypothesis, for samples from group A and for samples from group B. Poisson regression is then performed where the standard LRT is computed to test for differences in expression between the two groups.
(iii) FET
Under the random sample process, the number of reads from a gene follows a binomial distribution which can be approximated by a Poisson distribution. The group counts and the total counts are also Poisson distributed. FET is then used to calculate the probability of observing group counts as the observed or more extreme if the null hypothesis is true (no difference between the two groups in the expression of gene g) for each gene.
Testing
FET and large sample approximation, such as the LRT and Z score test, are the choices. Multiple testing was adjusted using the methods of either Benjamini and Hochberg [13] or Storey and Tibshirani [17].
3. baySeq
baySeq differs from the above two packages by employing an empirical Bayesian analysis approach to determine if there is DE between two different conditions [18]. It begins by assuming that the data follow a distribution, either Poisson or NB, which is defined by a set of underlying parameters. The prior for each parameter is estimated by first bootstrapping from the data, and then either applying the maximum likelihood method (assuming that the prior is from a Poisson or NB distribution), or applying quasi-likelihood methods. For each gene or tag, two scenarios (hypotheses) are envisaged: one where the expression pattern is the same across the two conditions (ie that there is no DE between the two conditions for the gene); the other where the expression patterns differ between the two conditions (ie that there is DE for the gene). Given the prior estimates and the likelihood of the distribution of the data, one can estimate the posterior likelihood under the two scenarios to determine if there is DE for that gene or tag. Two distributions are proposed for the data; one is Poisson and the other is NB. The Poisson model is faster, yet the NB model provides better fit for most RNA-seq data. baySeq recommends using the NB model in general. The required data format is the same as the edgeR package. Parallel processing is provided through the 'snow' package for faster processing.
Normalisation
No normalisation procedure is proposed in this package.
Models
Two models are used in the package; one is to assume a Poisson distribution on each tag that is Y
gij
~ (M
jpgi
,), where the prior for p
gi
is assumed to follow gamma distribution p
gi
~ Γ(α
gi
, βgi). This model is therefore named the Poisson-gamma approach. In general, a subset of data is taken initially and more than 5,000 iterations are recommended for the bootstrapping. Gamma parameters are calculated using maximum likelihood methods. The mean of the maximum likelihood estimates is taken to obtain a prior on p
gi
. An initial prior value needs to be provided for the program to start. The other model assumes that the data are NB distributed, Y
gij
~ NB(M
j
p
gi
, ϕ
g
). As there is no conjugate prior available for this distribution, a numerical solution for an empirical prior is required. The program first bootstraps from the data, with around 10,000 iterations suggested. The parameters for an empirical prior distribution are estimated using the quasi-likelihood approach. Given the prior and the likelihood of the data, the posterior likelihoods are then calculated. The program repeatedly bootstraps to improve the accuracy of the prior estimation, and the posterior likelihood is updated accordingly.
Testing
The estimated posterior likelihoods are reported on the natural logarithmic scale.