# BADGE: A novel Bayesian model for accurate abundance quantification and differential analysis of RNA-Seq data

- Jinghua Gu
^{1}, - Xiao Wang
^{1}, - Leena Halakivi-Clarke
^{2}, - Robert Clarke
^{2}and - Jianhua Xuan
^{1}Email author

**15(Suppl 9)**:S6

https://doi.org/10.1186/1471-2105-15-S9-S6

© Gu et al.; licensee BioMed Central Ltd. 2014

**Published: **10 September 2014

## Abstract

### Background

Recent advances in RNA sequencing (RNA-Seq) technology have offered unprecedented scope and resolution for transcriptome analysis. However, precise quantification of mRNA abundance and identification of differentially expressed genes are complicated due to biological and technical variations in RNA-Seq data.

### Results

We systematically study the variation in count data and dissect the sources of variation into between-sample variation and within-sample variation. A novel Bayesian framework is developed for joint estimate of gene level mRNA abundance and differential state, which models the intrinsic variability in RNA-Seq to improve the estimation. Specifically, a Poisson-Lognormal model is incorporated into the Bayesian framework to model within-sample variation; a Gamma-Gamma model is then used to model between-sample variation, which accounts for over-dispersion of read counts among multiple samples. Simulation studies, where sequencing counts are synthesized based on parameters learned from real datasets, have demonstrated the advantage of the proposed method in both quantification of mRNA abundance and identification of differentially expressed genes. Moreover, performance comparison on data from the Sequencing Quality Control (SEQC) Project with ERCC spike-in controls has shown that the proposed method outperforms existing RNA-Seq methods in differential analysis. Application on breast cancer dataset has further illustrated that the proposed Bayesian model can 'blindly' estimate sources of variation caused by sequencing biases.

### Conclusions

We have developed a novel Bayesian hierarchical approach to investigate within-sample and between-sample variations in RNA-Seq data. Simulation and real data applications have validated desirable performance of the proposed method. The software package is available at http://www.cbil.ece.vt.edu/software.htm.

## Background

Next Generation Sequencing (NGS) technology has opened a new era for transcriptome analysis, which grants the ability to investigate novel biological problems, such as alternative splicing, differential isoforms, gene fusion, etc. By piling up millions of reads along the reference genome, RNA-Seq technology can obtain signals in a much larger dynamic range with much higher accuracy compared to traditional microarray based technologies. For RNA-Seq analysis, the most popular routine is to determine the expression of genes (abundance quantification) and to identify differentially expressed genes (DEGs).

Methods that quantify gene/isoform level expression mainly fall into two categories: Poisson count mode (e.g., Cufflinks [1], etc) and linear regression model (e.g., Isolasso [2], SLIDE [3], BASIC [4], etc.). The major challenge in accurate quantification of gene expression is that large systematic bias in sequencing counts has been observed due to multiple factors. In contrast to uniform assumption of read distribution, it has been reported that sequence counts show a variety of physical and chemical biases, including transcript length bias, GC-content bias, random hexamer priming bias, etc [5–7].

Differential analysis of RNA-Seq data has been focused on modeling variance among biological replicates or samples in the same phenotype group. EdgeR [8] is the first method that models the 'between-sample' variability by replacing the Poisson model with Negative Binomial model. Various statistical methods have been proposed to model the variance among samples in biological groups, aimed to improve overall fitting of count data or robustness against outliers [9–12].

Despite initial success to model uncertainties associated with sequencing counts from different aspects, there lacks a systematic effort to address variability in RNA-Seq data. We dissect the variance in sequencing counts along two dimensions: over-dispersion of read counts within the same sample (i.e., within-sample variation) and over-dispersion of read counts among individuals from the same biological group (i.e., between-sample variation). Within-sample variation typically leads to large variance of read counts among genomic loci (e.g., nucleotides or exons), which have similar expression level in the same sample. It is typically caused by technical artifacts such as uncorrected systematic bias and gene specific random effects. On the other hand, between-sample variation is mostly due to biological differences among samples under the same condition. To increase the accuracy of abundance estimation so as to improve DEG identification, immediate attention is needed to developing a unified model that takes care of both forms of variation in RNA-Seq data. We propose a computational method, namely Bayesian Analysis of Dispersed Gene Expression ('BADGE'), to model variability in RNA-Seq data. A full Bayesian model is employed to simultaneously account for within-sample variation and between-sample variation to improve inference. The proposed method has several novel contributions compared to existing methodologies: 1) a unified Bayesian causality model is developed for joint abundance estimation and DEG identification. The improved accuracy in profiling mRNA abundance can facilitate the identification of DEGs, which may in turn refine the parameter learning in abundance quantification. 2) A Poisson-Lognormal regression model is incorporated to model within-sample variation [13]. Instead of dealing with multiples sources of technical bias and variation separately, the proposed method can 'blindly' detect over-dispersion pattern within the individual sample. 3) Gamma-Gamma model [14] is used to model between-sample variation, which accounts for over-dispersion of read counts among multiple samples. BADGE is a unified computational method that extensively models variability in RNA-seq data to improve abundance quantification and DEG identification.

## Methods

### Bayesian Analysis of Dispersed Gene Expression (BADGE)

We have developed a computational method, namely Bayesian Analysis of Dispersed Gene Expression (BADGE), to model extensive variability in RNA-Seq data. BADGE explicitly models both between-sample variation and within-sample variation to improve abundance quantification and DEG identification. In this paper, we only focus on the gene level analysis, while the concept can be straight-forwardly generalized for genes with multiple transcripts (isoforms).

*y*

_{ g,i,j }represent observed counts that fall into the

*i*

^{th}(1 ≤

*i*≤

*I*

_{g}) exon region of gene

*g*(1 ≤

*g*≤

*G*) in sample

*j*(1 ≤

*j*≤

*J*), which follows Poisson distribution with mean

*γ*

_{ g,i,j }.

*I*

_{ g }is the number of exons in gene

*g. G*is the total number of genes.

*J*=

*J*

_{1}+

*J*

_{2}is the total number of samples, where

*J*

_{1}and

*J*

_{2}denote samples in condition 1 and 2, respectively. Within-sample over-dispersion indicates that

*γ*

_{ g,i,j }has unknown heterogeneity across the gene rather than taking constant value. A hierarchical Bayesian model is constructed to model within-sample variation of RNA-Seq data as follows:

where *β*_{
g,j
}is the true expression level of gene *g* for sample *j. x*_{
g,i
}is the length of the *i*^{th} exon weighted by the library size of sample *j. U*_{
g,i,j
}is the unknown within-sample variation parameter, which follows normal distribution with mean 0 and precision *τ*. 'Flat' prior is assigned for *τ* by setting its shape *a* = 1 and rate *b* = 0. Equations (1-4) are also known as the Poisson-Lognormal regression model with identity link function.

*y*

_{ g,i,j }exhibits over-dispersion, but also

*β*

_{ g,j }has variation across multiple samples in the same biological group. To model between-sample variation carried by

*β*

_{ g,j }, we adopt the Gamma-Gamma model that is widely used in microarray gene expression analysis [14] into the Poisson count model. Let

*j*

_{1}and

*j*

_{2}represent samples in condition 1 and 2.

*d*

_{ g }is the binary differential state of gene

*g*, where

*d*

_{ g }= 0 means gene

*g*is not differentially expressed;

*d*

_{ g }= 1, otherwise. The Gamma-Gamma model for RNA-Seq differential expression is given by:

*d*

_{ g }= 0, ${\lambda}_{g}^{\left(1\right)}={\lambda}_{g}^{\left(2\right)}={\lambda}_{g}$; if ${d}_{g}=1,{\lambda}_{g}^{\left(1\right)}\ne {\lambda}_{g}^{\left(2\right)}$.

*α*is the shape parameter for

*β*

_{ g,j }, which does not depend on differential state

*d*

_{ g }Moreover, we assume that the pooled rate parameter

*λ*

_{ g }from the gene population further follows Gamma distribution with shape parameter

*α*

_{ 0 }and rate parameter

*v*. We assign non-informative priors for hyper-parameters

*α, α*

_{ 0 }

*, d*

_{ g }(

*P*(

*d*

_{ g }= 1)=

*π*

_{ g }= 0.5) and

*v*(

*a*

_{ 0 }= 1,

*b*

_{ 0 }= 0). The sub-model defined by Equations (5-9) considers between-sample variation within the same group, which borrows knowledge from the entire population to improve parameter estimation of individual genes. Figure 1 gives the Bayesian hierarchical dependency graph for all the parameters involved in the BADGE method using plate notation. There are three plates in Figure 1: The inner plate denotes dependency among read counts within

*I*

_{ g }exons in gene

*g*; middle plate denotes dependency among

*J*samples; and the outmost plate represents

*G*genes. Observation

**y**, represented by shaded circle, is the raw exon level RNA-Seq count (for gene level count data,

*I*

_{ g }= 1), which depends on mRNA abundance

**β**, gene design matrix

**x**and within-sample over-dispersion parameter

**U**.

**β**further depends on group level between-sample over-dispersion Gamma parameter

**λ**and

*α*, and

**U**depends on global within-sample over-dispersion parameter

*τ*. Parameter

**λ**is determined by gene level differential state

**d**, and its priors

*ν*and

*α*

_{ 0 }

*. d*

_{ g }(differential state of gene

*g*) is assumed to follow P(dg = 1) =

*π*

_{ g }= 0.5. Hyper-parameters

*a*

_{0},

*b*

_{ 0 }

*, π, a*, and

*b*, shown in shaded square, are fixed to construct non-informative priors (see additional file 1).

### Estimate model parameters using Gibbs sampling

**y**(read count) and

**x**(exon length) is given by:

*β*

_{ g,j }

*, U*

_{ g,i,j }and τ can be sampled from their corresponding conditional distributions as:

We pool all the samples from *U*_{
g,i,j
} and get its estimate ${\hat{U}}_{g,i,j}=\frac{1}{T-{t}_{b}+1}{\displaystyle \sum _{{t}_{b}}^{T}}{U}_{g,i,j}^{t},$ where ${U}_{g,i,j}^{t}$ denotes sampled *U*_{
g,i,j
} at sample *t. t*_{
b
} is the last sample when burn-in stops. *T* is the total number of Gibbs samples. ${\hat{U}}_{g,i,j}$ will be passed to the Gamma-Gamma model to estimate parameters associated with DEG identification.

**β, λ**, α, α

_{0},

*v*and

**d**according to their conditionals. The posterior distribution of ${\beta}_{j}\left({\beta}_{{j}_{1}},{\beta}_{{j}_{2}}\right)$ can be sampled from:

*d*

_{ g }= 1,

*d*

_{ g }= 0,

*etal*. [14], the posterior distribution of

**d**given

**β**,

*α, α*

_{0}and ν can be derived as:

*α*

_{0}and

*α*are given by:

We use Gibbs sampling method [4, 13, 15, 16] to estimate the posterior distributions of individual parameters iteratively from their complex joint distribution. For parameters **β**, *ν, τ*, and **λ**, we use conjugate priors to sample from their conditional distributions with standard probability distributions (Gamma distribution). For parameters (**U**, *α*_{0} and *α*) that do not have conjugate priors, we use Metropolis-Hastings sampling to sample their posterior distributions. Please see more details about the implementation of Gibbs sampler in additional file 1.

## Results and discussion

### Over-dispersion of RNA-Seq counts in two dimensions: between-sample variation and within-sample variation

### Generate simulation data based on real RNA-Seq datasets

*et al*. [10] to first estimate model parameters from real datasets and then use them to generate sequencing counts based on human annotation file (version: GRCh37/hg19). Two RNA-Seq datasets were used in the study: 1) a mouse dataset with 10 C57BL/6J (B6) mouse samples and 11 DBA/2J (D2) mouse samples [20]; 2) 23 basal breast cancer samples from the TCGA project [18]. For the TCGA dataset, we divided the patients (which received chemotherapy treatment) into two groups: early re-current group (recurrence time < 2 yrs, 13 samples) and late recurrent group (recurrence time > 3 yrs, 10 samples). Figure 3 gives the trace plots of sampled model parameters. We used thin = 10 for the Gibbs sampling process, which means to record every 10

^{th}sample. It took BADGE several hours to estimate posterior distributions using 10,000 iterations on real dataset. We have also plotted auto-correlation curves of each parameter in supplementary materials (additional file 1). We further explored the variability of RNA-Seq data by close examination of estimated model parameters. For within-sample variability, two typical values of over-dispersion prior parameter

*τ*have been estimated, where

*τ*= 0.44 in mouse dataset and

*τ*= 1.78 in TCGA dataset. In our Poisson-Lognormal regression model,

*τ*is the precision parameter (inverse of Gaussian variance

*σ*

^{2}) that controls overall degree of within-sample over-dispersion. Smaller

*τ*indicates larger variation of read counts across the transcriptome. Small value of

*τ*observed in real datasets strongly supported our motivation that read counts along one transcript must be corrected for improved abundance estimation. Between-sample variability is jointly determined by

*α, α*

_{0}, and

*ν*. Small value in

*ν*(0.001 in mouse and 0.2 in TCGA sample) yield large variation of VAR (

*λ*), where VAR (

*β*) increases when

*λ*takes small value. We used the model parameters estimated from real datasets to generate read counts for simulation. Please refer to supplementary materials in additional file 1 for more detailed description of simulation data.

### Performance comparison for abundance estimation on simulation data

*τ*from real datasets to generate synthetic data and studied the performance of the BADGE method for abundance estimation. We compared our method with four commonly used methods for RNA-Seq normalization, which were: reads-per-kilobase-per-million (RPKM), DESeq, trimmed mean of M-values (TMM) and upper quantile normalization (UQUA). RPKM [22] is calculated through normalizing reads by length of genomic features (genes and exons) and total library size. DESeq normalization is implemented by DESeq (1.14.1) [9], which is a differential gene identification method based on negative binomial model. TMM was originally developed in edgeR [8] and later included into BioConductor package NOISeq [23]. NOISeq (2.0.0) also has separate implementations of RPKM and UQUA methods, which were used here for performance comparison. Based on estimated parameters from real RNA-Seq datasets (mouse and TCGA breast cancer data), we selected typical model parameters to simulate count data: σ = 0.75, 1 and 1.5; ν = 0.1 and 0.001; α = 2 and α

_{0}= 0.5. We used the average correlation of normalized counts with ground truth gene expression to measure the accuracy of abundance quantification across multiple samples. Figure 4 gives the average correlation for all competing methods under different model settings. From Figure 4, we see that the BADGE method had robust performance under different over-dispersion settings by maintaining a performance measure (correlation to ground truth) very close to 1. Among the rest of the normalization methods, DESeq, TMM and UQUA achieved comparable performance across multiple parameter settings, while RPKM had the least favourable performance in all scenarios. Our computational result is quite consistent with the observation by Dillies

*et al*. [24] that DESeq and TMM (edgeR) are much better normalization methods than RPKM.

### Performance comparison for differentially expressed gene (DEG) identification on simulation data

*τ*(or standard deviation σ), and we set

*τ*= 1.78 (i.e.,

*σ*= 0.75, which was learned from TCGA basal dataset) and

*τ*= 0.44 (i.e.,

*σ*= 1.5, which was learned from mouse dataset.). According to estimated parameters from real datasets, we set

*α*= 2,

*α*

_{0}= 0.5, while varied ν between 0.001 and 0.1, which was consistent with the parameter settings in abundance estimation. For each parameter set, we randomly selected 10 genesets from hg19 annotation file to evaluate the variance of the performance. Area-under-the-curve (AUC) of the receiver operating characteristic (ROC) curve was used as performance measurement. Tables 1, 2, 3 give the AUC values for each method along with standard deviations (listed in parentheses) of AUCs across 10 experiments.

Performance comparison using AUC for DEG identification (highly differentially expressed genes)

σ (τ) |
| BADGE | DESeq | edgeR | DSS | EBSeq |
---|---|---|---|---|---|---|

0.75 | 0.001 |
| 0.921 | 0.919 | 0.881 | 0.928 |

(1.78) | (0.018) | (0.017) | (0.018) | (0.055) | (0.017) | |

0.1 |
| 0.894 | 0.895 | 0.889 | 0.906 | |

(0.018) | (0.018) | (0.020) | (0.024) | (0.021) | ||

1.5 | 0.001 |
| 0.881 | 0.875 | 0.809 | 0.890 |

(0.44) | (0.014) | (0.029) | (0.032) | (0.082) | (0.030) | |

0.1 |
| 0.868 | 0.865 | 0.858 | 0.875 | |

(0.023) | (0.027) | (0.025) | (0.027) | (0.022) |

Performance comparison using AUC for DEG identification (moderately differentially expressed genes)

σ (τ) |
| BADGE | DESeq | edgeR | DSS | EBSeq |
---|---|---|---|---|---|---|

0.75 | 0.001 |
| 0.724 | 0.753 | 0.739 | 0.750 |

(1.78) | (0.071) | (0.070) | (0.081) | (0.197) | (0.069) | |

0.1 |
| 0.768 | 0.760 | 0.852 | 0.781 | |

(0.076) | (0.051) | (0.066) | (0.071) | (0.044) | ||

1.5 | 0.001 |
| 0.691 | 0.699 | 0.790 | 0.725 |

(0.44) | (0.056) | (0.041) | (0.042) | (0.088) | (0.035) | |

0.1 |
| 0.743 | 0.729 | 0.798 | 0.748 | |

(0.080) | (0.070) | (0.081) | (0.074) | (0.058) |

Performance comparison using AUC for DEG identification (weakly differentially expressed genes)

σ (τ) |
| BADGE | DESeq | edgeR | DSS | EBSeq |
---|---|---|---|---|---|---|

0.75 | 0.001 |
| 0.656 | 0.680 | 0.603 | 0.687 |

(1.78) | (0.099) | (0.058) | (0.084) | (0.133) | (0.031) | |

0.1 |
| 0.659 | 0.663 | 0.695 | 0.684 | |

(0.092) | (0.043) | (0.062) | (0.075) | (0.049) | ||

1.5 | 0.001 |
| 0.641 | 0.647 | 0.635 | 0.669 |

(0.44) | (0.080) | (0.028) | (0.071) | (0.134) | (0.029) | |

0.1 |
| 0.687 | 0.680 | 0.724 | 0.682 | |

(0.109) | (0.072) | (0.063) | (0.104) | (0.073) |

We simulated RNA-Seq gene expression with three different scenarios: genes that were highly differentially expressed, moderately differentially expressed and weakly differentially expressed (see supplementary materials in additional file 1 for more information). From Tables 1, 2, 3, we see that the BADGE method consistently outperformed existing methods in different parameter settings (highlighted in bold). For highly differentially expressed genes (Table 1), the performance of the other methods degraded as we decreased *τ*, while BADGE was able to maintain good performance by employing a Poisson-Lognormal model to account for within-sample variability. For genes that were weakly differentially expressed (Table 3), BADGE achieved a maximum improvement of AUC up to 1.3, compared to the second best method EBSeq (Table 3, *τ* = 0.44, *ν* = 0.001).

### Performance comparison for differentially expressed gene (DEG) identification on Sequencing Quality Control (SEQC) data

On the SEQC dataset, BADGE had the best performance among all five methods by achieving an AUC very close to 0.9. The second best method was DSS with an AUC about 0.85. DESeq (AUC = 0.7624) and edgeR (AUC = 0.7675) had very close performance, which was consistent with the previous results reported by Rapaport *et al*. [25]. EBSeq, on the other hand, had the least favourable performance (AUC = 0.71) on this specific dataset and it failed to detect the most strongly differentially expressed genes: its sensitivity was less than 0.1 when its specificity was about 0.9. By close examination of the 'left' ROC curves of the five methods, we can further infer that BADGE should have significantly better precision than the competing methods, whose sensitivity went all the way up to 0.7 before any sacrifice in specificity.

### 'Blind' estimation of hidden heterogeneity in RNA-Seq data

*U*

_{ g,i,j }(gene

*g*, exon

*i*, sample

*j*), while the overall degree of variation is controlled by

*τ*. Based on sampled parameters from real datasets, we further investigated

*U*

_{ g,i,j }to see how systematic artifacts in RNA-Seq (such as transcript length bias and GC content) can be de-convoluted by BAGDE method. Figure 6 shows the histogram of Pearson's correlation between estimated

*U*

_{ g,i,j }and: 1. transcript location; 2. GC content. From Figure 6 (a), we see strong positive correlation between estimated over-dispersion parameter

*U*

_{ g,i,j }and transcript location, which indicates that most of the genes had biased expression towards 3'-end of the transcript in our dataset, while about 15 percent of the genes had low correlation (<0.5). In addition, a significant correlation between

*U*

_{ g,i,j }and GC-content bias (extremely low or high GC content are associated with low abundance [26]) were also observed in Figure 6 (b). These observations support that BADGE can correctly estimate hidden sources of variation in RNA-Seq data 'blindly' without using transcript location information or sequence information.

## Conclusions

Large variation in RNA-Seq data has become the major obstacle against accurate estimation of gene expression and DEG identification. Much effort has been made to model variation across biological replicates, while limited attention is paid to tackle extensive over-dispersion observed in sequencing counts. For short-read sequencing technologies (e.g., Illumina), multiple sources of systematic bias have been identified, including transcript length bias, GC-content bias, etc. However, in-depth investigation of real RNA-Seq datasets has revealed the following complications: 1) Sequencing bias not only changes from one gene to another, but also varies among samples (Figure 2(c) and (d)); 2) Gene expression is jointly influenced by multiple bias factors, which leads to large variation across the entire transcriptome (Figure 6). However, current research activities have been focused on addressing individual bias corrections, which lacks a unified effort to account for total variability in RNA-Seq data. Therefore, we propose the BADGE method to extensively model both within-sample variability (bias and random variance) and between-sample variability (biological variations among replicates or within the same phenotype group) to improve quality of inference.

## Declarations

The publication costs for this article were funded by National Institutes of Health (NIH) [CA149653].

## Declarations

### Acknowledgements

This work is supported by National Institutes of Health (NIH) [CA149653, CA149147, CA164384, and NS29525-18A, in part].

This article has been published as part of BMC Bioinformatics Volume 15 Supplement 9, 2014: Proceedings of the Fourth Annual RECOMB Satellite Workshop on Massively Parallel Sequencing (RECOMB-Seq 2014). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/15/S9.

## Authors’ Affiliations

## References

- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28 (5): 511-515.PubMed CentralView ArticlePubMedGoogle Scholar
- Li W, Feng J, Jiang T: IsoLasso: a LASSO regression approach to RNA-Seq based transcriptome assembly. J Comput Biol. 2011, 18 (11): 1693-1707.PubMed CentralView ArticlePubMedGoogle Scholar
- Li JJ, Jiang CR, Brown JB, Huang H, Bickel PJ: Sparse linear modeling of next-generation mRNA sequencing (RNA-Seq) data for isoform discovery and abundance estimation. Proc Natl Acad Sci USA. 2011, 108 (50): 19867-19872.PubMed CentralView ArticlePubMedGoogle Scholar
- Zheng S, Chen L: A hierarchical Bayesian model for comparing transcriptomes at the individual transcript isoform level. Nucleic Acids Res. 2009, 37 (10): e75-PubMed CentralView ArticlePubMedGoogle Scholar
- Wu Z, Wang X, Zhang X: Using non-uniform read distribution models to improve isoform expression inference in RNA-Seq. Bioinformatics. 2011, 27 (4): 502-508.View ArticlePubMedGoogle Scholar
- Hansen KD, Irizarry RA, Wu Z: Removing technical variability in RNA-seq data using conditional quantile normalization. Biostatistics. 2012, 13 (2): 204-216.PubMed CentralView ArticlePubMedGoogle Scholar
- Hansen KD, Brenner SE, Dudoit S: Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res. 2010, 38 (12): e131-PubMed CentralView ArticlePubMedGoogle Scholar
- Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26 (1): 139-140.PubMed CentralView ArticlePubMedGoogle Scholar
- Anders S, Huber W: Differential expression analysis for sequence count data. Genome Biol. 2010, 11 (10): R106-PubMed CentralView ArticlePubMedGoogle Scholar
- Wu H, Wang C, Wu Z: A new shrinkage estimator for dispersion improves differential expression detection in RNA-seq data. Biostatistics. 2013, 14 (2): 232-243.PubMed CentralView ArticlePubMedGoogle Scholar
- Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, Haag JD, Gould MN, Stewart RM, Kendziorski C: EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013, 29 (8): 1035-1043.PubMed CentralView ArticlePubMedGoogle Scholar
- Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L: Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013, 31 (1): 46-53.View ArticlePubMedGoogle Scholar
- Hu M, Zhu Y, Taylor JM, Liu JS, Qin ZS: Using Poisson mixed-effects model to quantify transcript-level gene expression in RNA-Seq. Bioinformatics. 2012, 28 (1): 63-68.PubMed CentralView ArticlePubMedGoogle Scholar
- Wei Z, Li H: A Markov random field model for network-based analysis of genomic data. Bioinformatics. 2007, 23 (12): 1537-1544.View ArticlePubMedGoogle Scholar
- Lawrence CE, Altschul SF, Boguski MS, Liu JS, Neuwald AF, Wootton JC: Detecting subtle sequence signals: a Gibbs sampling strategy for multiple alignment. Science. 1993, 262 (5131): 208-214.View ArticlePubMedGoogle Scholar
- Sabatti C, James GM: Bayesian sparse hidden components analysis for transcription regulation networks. Bioinformatics. 2006, 22 (6): 739-746.View ArticlePubMedGoogle Scholar
- Wang Z, Gerstein M, Snyder M: RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009, 10 (1): 57-63.PubMed CentralView ArticlePubMedGoogle Scholar
- Comprehensive molecular portraits of human breast tumours. Nature. 2012, 490 (7418): 61-70.Google Scholar
- Cheung VG, Nayak RR, Wang IX, Elwyn S, Cousins SM, Morley M, Spielman RS: Polymorphic cis- and trans-regulation of human gene expression. PLoS Biol. 2010, 8 (9):Google Scholar
- Bottomly D, Walter NA, Hunter JE, Darakjian P, Kawane S, Buck KJ, Searles RP, Mooney M, McWeeney SK, Hitzemann R: Evaluating gene expression in C57BL/6J and DBA/2J mouse striatum using RNA-Seq and microarrays. PLoS One. 2011, 6 (3): e17820-PubMed CentralView ArticlePubMedGoogle Scholar
- Edgren H, Murumagi A, Kangaspeska S, Nicorici D, Hongisto V, Kleivi K, Rye IH, Nyberg S, Wolf M, Borresen-Dale AL: Identification of fusion genes in breast cancer by paired-end RNA-sequencing. Genome Biol. 2011, 12 (1): R6-PubMed CentralView ArticlePubMedGoogle Scholar
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5 (7): 621-628.View ArticlePubMedGoogle Scholar
- Tarazona S, Furio-Tari P, Ferrer A, Conesa A: NOISeq: Exploratory analysis and differential expression for RNA-seq data. R package version 200. 2012Google Scholar
- Dillies MA, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, Keime C, Marot G, Castel D, Estelle J: A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2013, 14 (6): 671-683.View ArticlePubMedGoogle Scholar
- Rapaport F, Khanin R, Liang Y, Pirun M, Krek A, Zumbo P, Mason CE, Socci ND, Betel D: Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. Genome Biol. 2013, 14 (9): R95-PubMed CentralView ArticlePubMedGoogle Scholar
- Benjamini Y, Speed TP: Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic Acids Res. 2012, 40 (10): e72-PubMed CentralView ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.