Improving RNA-Seq expression estimation by modeling isoform- and exon-specific read sequencing rate
- Xuejun Liu^{1}Email author,
- Xinxin Shi^{1},
- Chunlin Chen^{1} and
- Li Zhang^{1}
Received: 24 November 2014
Accepted: 24 September 2015
Published: 16 October 2015
Abstract
Background
The high-throughput sequencing technology, RNA-Seq, has been widely used to quantify gene and isoform expression in the study of transcriptome in recent years. Accurate expression measurement from the millions or billions of short generated reads is obstructed by difficulties. One is ambiguous mapping of reads to reference transcriptome caused by alternative splicing. This increases the uncertainty in estimating isoform expression. The other is non-uniformity of read distribution along the reference transcriptome due to positional, sequencing, mappability and other undiscovered sources of biases. This violates the uniform assumption of read distribution for many expression calculation approaches, such as the direct RPKM calculation and Poisson-based models. Many methods have been proposed to address these difficulties. Some approaches employ latent variable models to discover the underlying pattern of read sequencing. However, most of these methods make bias correction based on surrounding sequence contents and share the bias models by all genes. They therefore cannot estimate gene- and isoform-specific biases as revealed by recent studies.
Results
We propose a latent variable model, NLDMseq, to estimate gene and isoform expression. Our method adopts latent variables to model the unknown isoforms, from which reads originate, and the underlying percentage of multiple spliced variants. The isoform- and exon-specific read sequencing biases are modeled to account for the non-uniformity of read distribution, and are identified by utilizing the replicate information of multiple lanes of a single library run. We employ simulation and real data to verify the performance of our method in terms of accuracy in the calculation of gene and isoform expression. Results show that NLDMseq obtains competitive gene and isoform expression compared to popular alternatives. Finally, the proposed method is applied to the detection of differential expression (DE) to show its usefulness in the downstream analysis.
Conclusions
The proposed NLDMseq method provides an approach to accurately estimate gene and isoform expression from RNA-Seq data by modeling the isoform- and exon-specific read sequencing biases. It makes use of a latent variable model to discover the hidden pattern of read sequencing. We have shown that it works well in both simulations and real datasets, and has competitive performance compared to popular methods. The method has been implemented as a freely available software which can be found at https://github.com/PUGEA/NLDMseq.
Keywords
Transcript expression Gene expression RNA-Seq data analysis Latent dirichlet allocation Probabilistic modelBackground
RNA-Seq, based on high-throughput sequencing technology, is nowadays a commonly used method to study transcriptome. An RNA-Seq experiment typically produces millions or billions of short sequenced reads. By counting the reads aligned to a reference transcriptome, the expression of the related transcripts can be calculated [1]. As alternative splicing, where various protein isoforms are yielded due to the different exon constitutions for a gene, receives more and more interest in the research of biomedicine, RNA-Seq has been widely used to quantify gene and isoform expression.
As read counts are proportional to the mRNA fragment abundance of a gene, the reads per kilobase of transcript per million mapped reads (RPKM) has been commonly used to represent gene expression [2]. Due to the existence of alternative splicing, multiple isoforms can share exons of a gene (we use isoforms to represent splice variants hereafter for simplicity). Therefore, there are a large number of reads which cannot be uniquely aligned to the isoform they originate from. For this reason, reads mapped to the shared exon need to be deconvoluted before computing isoform expression. Inferring the mechanism of dispatching mapped reads to their origins is a difficulty in isoform expression calculation from RNA-Seq data. A number of approaches make use of the additive property of Poisson distribution to address this problem [3–5]. Given a transcriptome assembly, the sequencing rate for a read count can be factorized as a sum of contributions from all isoforms to which this read can be mapped. Another group of methods employ latent variable models to discover the unknown isoform expression [6–8]. By introducing discrete latent variables to representing isoforms each read originates from, the sequencing process of each individual read is modeled.
In addition to read mapping ambiguity, another difficulty of expression calculation from RNA-Seq data is the non-uniformity of read distribution along the reference genome and transcriptome due to positional, sequencing, mappability and other undiscovered sources of biases [9]. The non-uniform read distribution violates the assumption of the direct RPKM representation and many expression calculation approaches, such as Poisson-based models. In order to obtain accurate expression estimates, these biases have to be accounted in expression calculation. Consequently, many bias correction strategies have been proposed to relieve the influence of non-uniformity of read distribution [5, 9–13]. Most of these approaches estimates the overall biases based on surrounding sequence contents or empirical data, and bias models are mostly shared by all genes. For example, the commonly used method, Cufflinks [12, 14], used a variable length Markov model to learn sequence-specific and positional biases based on the surrounding sequences of all reads in a set of empirical data. This model was then shared for all genes to make bias correction. Therefore, if two reads have similar nucleotide constitution and similar relative positions along transcriptome reference, they would have the similar biases under this model assumption. A recent study found that the read distribution is isoform-specific and a Poisson-based model, Sequgio, is proposed to jointly estimate isoform expression and isoform-specific read sequencing bias [15].
In this paper, we propose a latent variable model, NLDMseq (normalized latent Dirichlet-Multinomial model for RNA-Seq data), to estimate gene and isoform expression. Given known annotations, this model introduces latent variables to represent the unknown isoforms, from which reads originate, and the proportion of multiple spliced isoforms. The isoform- and exon-specific read sequencing rates are modeled as parameters which can be estimated from data. We utilize the replicate information of multiple lanes for a single library to identify read sequencing rates. NLDMseq benefits from the latent variable model that would help to discover the hidden pattern of read sequencing. Another advantage of NLDMseq is that it models normalized read counts for each exon rather than each individual read as many latent variable models do [6–8]. Therefore, NLDMseq does not need to deal with the nucleotide constitution of each read, and thus achieves fast computation. We employ simulation and real data to verify the performance of our method in terms of accuracy in the calculation of gene and isoform expression.
Methods
Latent dirichlet-multinomial model
It has been found by previous studies that there exist similarities among the patterns of count variation for many genes or isoforms among samples [10, 15]. There is therefore a need to model read distribution bias for each individual gene or isoform. To reach this goal, it would be ideal to model the sequencing rate for each individual read of a particular gene or isoform. However, this would bring tremendous computation in light of the huge amount of read data. Also, for many lowly expressed genes overfitting would be a problem with little data available. Considering computational efficiency of the model, we count read number for each exon in the transcriptome reference and take exon as the basic unit of data to model exon-specific read sequencing rate for each isoform. As long reads and paired-end reads become more and more popular in RNA-seq protocal, we use fragment to represent any transcript sequence found in sample and are then interested in fragment number for each exon hereafter. We aim to model the stochastic mechanism in generating the “exon corpora” of RNA-seq data for each gene.
We denote e as an indicator vector with E entries, among which there is only a single entry with value 1 and all the others are zero. Thus, each vector represents an exon. For each observed exon e, we assume there is an associated latent variable t indicating the unknown isoform which generates this exon. The K-vector t has value one for the kth entry if the exon is generated by isoform k, and 0 for others. For the purpose of accommodating the random noise in fragment sequencing and fragments generated from undiscovered isoforms, we introduce an additional isoform which contains all exons of this gene. If the possibility of the generation for a particular exon from all known isoforms is considerably low, it may be assigned to this special isoform automatically. Figure 1 shows the construction of gene structure and fragment counting strategy of an example gene. A K-vector variable θ is used to denote the percentage of the K isoforms, where θ _{ i }>0 and \(\sum _{i}\theta _{i}=1\). We assume a Dirichlet prior on θ.
- 1.
Generate θ _{ l } once for lane l, θ _{ l }∼Dirichlet(α).
- 2.
For lane l, generate t _{ ln }, t _{ ln }∼Multinomial(1,θ _{ l }).
- 3.
Generate exon e _{ ln } given t _{ ln }, \(e_{\textit {ln}}\sim \text {Multinomial}(1,\beta ^{t_{\textit {ln}}})\phantom {\dot {i}\!}\).
We notice that Sequgio proposed in [15] also modeled the isoform-specific read sequencing bias. Our model differs from Sequgio in the fact that Sequgio uses the simple Poisson distribution to model the generated reads and does not model the variability in the selection of isoform which generates each read. Also, we do not share the isoform-specific read sequencing rate β across samples as Sequgio in case that biological variation may violate the conserved read count distribution for some genes. Biological samples are therefore processed individually with our approach. We aim to provide accurate expression measurements at this stage and intend to handle the biological replicate noise in the downstream analysis (see the section of “Application to DE analysis”). NLDMseq is also distinct from other latent variable models, such as RSEM [6], MISO [7] and Bitseq [8], which explicitly modeled the generation process of each read by examining the the start sequencing position and nucleotide components of each read. These approaches either assume a uniform read distribution along isoform sequence [7] or need to deal with nucleotide constitution in order to correct the bias in read distribution [6, 8]. Consequently, this would increase the computation cost of models. In contrast, our method models the frequency of the existence of each exon and this is able to reduce the computation load of our model. Meanwhile, the isoform- and exon-specific read sequencing rate is also obtained and this accounts for the non-uniformity of the primary read distribution.
Variational EM solution
where < · > represents expectation and q denotes the variational distribution of θ _{ l }. We omit the variational parameters of q distribution for simplification. The Eqs. (1) and (2) are updated iteratively until convergence is obtained.
If isoform k includes exon j, β _{ kj } is updated using Eq. (5), otherwise it is confined to be zero. The hyperparameter α cannot be solved analytically and therefore can be worked using the efficient Newton-Raphson method which was used in [16]. The E-step and M-step are repeated until a stable solution is achieved. Readers can refer to [16] for details of the variational EM algorithm.
Gene and isoform expression representation
where N _{ s } is the total number of fragments obtained for this sample. Gene expression can be represented by \(\sum _{k}f_{k}\). In order to obtain a level of measurement uncertainty for both isoform and gene expression, we draw M samples from the distribution p(θ|α) and calculate isoform and gene expression for each sample θ ^{ m }, where 1≤m≤M. Thus, we have M estimates for each isoform and gene expression. Using these samples, we can obtain the variance of the estimated expression which can be propagated into downstream analysis to obtain improved results. For single-isoform genes, the expression is calculated directly using FPKM assuming the uniform distribution of the reads.
Software
We have implemented NLDMseq in a free Python/C tool, which is available from https://github.com/PUGEA/NLDMseq, for public use. Any aligner that can align the raw experimental reads to reference transcriptome can be used for the preprocess of RNA-seq data for NLDMseq. We use Bowtie 2 [17] throughout this paper since it is sensitive to gaps and good at aligning long reads, and is used by many approaches, such as Bitseq [8] and MMSEQ [4]. After the alignment, our software can be used to calculate gene and isoform expression. Our software includes two parts. One is Python scripts which are used to expand gene models and obtain the fragment count for each exon from the alignment output file. The other is C codes, which solve the NLDMseq model and compute expression measurements. For simplicity of usage, the software can be called in a single run. For details of using our software, readers can refer to the documentation and examples on the software’s website.
Datasets
We use the well-studied Microarray Quality Control (MAQC) dataset [18] to validate gene expression estimation from NLDMseq. Gene expression from high-quality RNA samples is measured in MAQC project to assess the comparability across multiple platforms mainly including various microarray and next-generation sequencing technologies. Two RNA samples, the universal human reference (UHR) RNA and the human brain reference (HBR) RNA, from Illumina platform are selected in our work to verify our new method, including single-end data (SRA010153) and paired-end data (SRA012427). The SRA010153 data contains two samples, HBR and UHR, while the SRA012427 data contains a single sample UHR. Besides genes in the whole human genome measured by microarray and RNA-Seq, around one thousand genes have been measured by qRT-PCR experiments. These genes can be used as gold standard to validate gene expression estimation obtained from other platforms. This dataset has been widely used as the benchmark to verify various analysis methods for microarray and RNA-Seq technologies [19–21]. Among the qRT-PCR validated genes, we used the Ensembl annotation data (GRCh37/hg19) and obtained 740 matching multi-isoform genes (see Additional file 5: Table S1). These genes are used to evaluate the performance of our approach at gene level. We further use the method in [20] to filter 198 differential expression (DE) genes and 81 non-DE multi-isoform genes (see Additional file 2: Table S2) with high confidence according to qRT-PCR measurements. Data of these 279 qRT-PCR validated genes is used as a gold standard to evaluate the sensitivity and the specificity of our method applied to DE analysis.
A real study on molecular injury in response to tobacco smoke exposure and lung cancer pathogenesis (SELC) [22] is used to further evaluate the performance of our approach on gene expression calculation. This experiment involves four phenotypes, healthy never smokers (NS), current smokers (S), smokers with (C) and without (NC) lung cancer undergoing lung nodule resection surgery. RNA-seq libraries were prepared and sequenced to obtain 22 million 75 nt paired-end reads for each sample by using the standard Illumina mRNA-seq protocol. The study considered differential gene expression in the comparisons, S vs. NS and C vs. NC. The RNA-seq measurements of eight genes were validated using qRT-PCR data. We aligned the reads to Ensembl annotation (GRCh37/hg19) and found one of the eight validated genes, NFKB1A1, not be annotated. We therefore exclude this gene and use the obtained fold changes of the other seven genes in expression comparison.
We also use a publicly available human breast cancer (HBC) dataset [23] to verify the performance of NLDMseq on isoform expression estimation. This dataset includes two samples, human breast cancer cell line (MCF-7) and normal cell line (HME). Four multi-isoform genes (TRAP1, ZNF581, WISP2 and HIST1H2BD) were validated by qRT-PCR experiments [24]. Each gene has two isoforms which have been interrogated for both cell lines. According the experiment results in [24], the regulations of these eight isoforms in the two samples and the regulations of the two isoforms for the same gene in the same sample can be calculated. We therefore obtain 16 qRT-PCR validated regulations which can be used to verify the performance of various approaches on isoform expression calculation (see Additional file 7: Table S3).
We generated simulated data using our model based on the calculated parameters from the MAQC dataset. We select all 1,604 multi-isoform genes and 12,064 isoforms of chromosome 1 from Ensembl human genome assembly GRCh37/hg19 and simulate two samples with single-end and paired-end data respectively. For the single-end data, we generate 3,673,804 single-end reads with seven lanes using parameters \(\hat {\beta }\) calculated from the single-end data for HBR sample. For the paired-end data we generate 125,126 reads with three lanes using parameters \(\hat {\beta }\) calculated from the paired-end data for UHR sample. For both cases, the values of parameters \(\hat {\alpha }\) are randomly generated. For each simulated read, we first sample θ _{ l } from \(\text {Dirichlet}(\hat {\alpha })\). Second, for each lane generate t _{ ln } from Multinomial(1,θ _{ l }). Next, given the generated t _{ ln } generate exon e _{ ln } from \(\text {Multinomial}(1,\beta ^{t_{\textit {ln}}})\phantom {\dot {i}\!}\). Finally, uniformly sample the starting position along exon e _{ ln } for a single-end read. To generate a pair of paired-end reads, we sample the length of the sequenced fragment from N(206,19.6) and then uniformly sample the starting position of this fragment along exon e _{ ln }. With the sampled fragment position and length, the paired reads for both ends can be simulated.
Results and discussion
Chart displaying the various study goals (row names) and the datasets (column names) used to accomplish each goal. Details of each dataset can be found in the section of Datasets
MAQC | SELC | HBC | Simulation | |
---|---|---|---|---|
Compare gene expression calculation | √ | √ | √ | |
Compare isoform expression calculation | √ | √ | ||
Compare DE detection performance | √ | |||
Compare computation efficiency | √ |
Application to real data
We evaluate the proposed approach, NLDMseq, on the estimation of gene and isoform expression in terms of accuracy using three real datasets and considering both single-end and paired-end data.
Estimated isoform- and exon-specific read sequencing rate
Accuracy on gene expression estimation
The Pearson correlation coefficients of estimated expression with qRT-PCR measurements using MAQC dataset
Dataset | Cufflinks | RSEM | MMSEQ | NLDMseq |
---|---|---|---|---|
SRA010153 (HBR) | 0.8033 | 0.8117 | 0.8000 | 0.8442 |
SRA010153 (UHR) | 0.8238 | 0.8265 | 0.8356 | 0.8585 |
SRA012427 (UHR) | 0.8107 | 0.8345 | 0.8430 | 0.8481 |
Comparison results of various methods for SELC dataset on calculated log _{2} fold change of gene expression
Comparison | Gene | qRT-PCR | Cufflinks | RSEM | MMSEQ | NLDMseq |
---|---|---|---|---|---|---|
S vs. NS | S100A8 | +1.15 | +1.13(0.02) | +2.72(1.38) | +3.20(1.78) | +1.65(0.44) |
S100A9 | +0.83 | +1.69(1.03) | +1.79(1.16) | +1.86(1.24) | +1.75(1.10) | |
CYP4F2 | +1.27 | +4.12(2.24) | +3.59(1.82) | +9.51(6.49) | +2.09(0.64) | |
C vs. NC | CCL20 | +4.01 | +6.00(0.50) | +4.66(0.16) | +5.32(0.32) | +7.08(0.76) |
IL8 | +1.50 | +1.51(0.01) | +1.21(0.19) | +1.67(0.11) | +1.05(0.30) | |
SCGB3A1 | +0.51 | +1.65(2.23) | +1.76(2.45) | +1.85(2.63) | +1.82(2.57) | |
SCGB1A1 | −0.48 | +1.37(3.89) | +1.56(4.26) | +1.02(3.13) | +1.62(4.39) | |
Average AER | NA | 1.42 | 1.63 | 2.24 | 1.46 |
Accuracy on isoform expression estimation
Comparison results of various methods for HBC dataset on calculated log _{2} fold change of isoform expression
Comparison | qRT-PCR | Cufflinks | RSEM | MMSEQ | NLDMseq |
---|---|---|---|---|---|
Case 1 | −0.4 | −0.81(1.03) | −0.41(0.02) | −0.74(0.85) | −0.91(1.28) |
Case 2 | −0.5 | −0.73(0.46) | −0.83(0.66) | −0.57(0.14) | −0.59(0.18) |
Case 3 | −0.9 | +4.52(6.02) | +4.91(6.46) | +4.42(6.02) | −0.34(0.62) |
Case 4 | −1.0 | +3.36(4.36) | +4.44(5.44) | +4.59(5.59) | −0.04(0.96) |
Case 5 | −0.3 | −1.74(4.80) | −0.93(2.10) | −1.29(3.30) | −2.41(7.03) |
Case 6 | −1.0 | −1.11(0.11) | −0.71(0.29) | −1.08(0.08) | −1.29(0.29) |
Case 7 | +1.2 | +1.12(0.07) | +1.34(0.12) | +1.47(0.23) | +0.37(0.69) |
Case 8 | +1.0 | +1.26(0.26) | +1.56(0.57) | +1.67(0.67) | +1.37(0.37) |
Case 9 | −5.6 | −5.70(0.02) | −5.40(0.04) | −5.50(0.02) | −7.76(0.39) |
Case 10 | −4.5 | −4.67(0.04) | −4.71(0.05) | −4.83(0.07) | −5.61(0.25) |
Case 11 | +0.4 | +0.52(0.30) | 0.01(0.98) | +0.34(0.15) | +0.26(0.35) |
Case 12 | +1.5 | +1.71(0.14) | +0.81(0.46) | +1.00(0.33) | +2.38(0.59) |
Case 13 | −4.7 | −2.74(0.42) | −2.92(0.38) | −4.28(0.09) | −1.07(0.77) |
Case 14 | −5.2 | −5.89(0.13) | −4.51(0.13) | −4.70(0.10) | −4.26(0.18) |
Case 15 | −5.4 | +3.57(1.66) | +1.83(1.34) | +1.79(1.33) | −0.21(0.96) |
Case 16 | −5.9 | +1.11(1.19) | +0.21(1.04) | +1.37(1.23) | −3.33(0.44) |
No. of inconsistency | NA | 4(1.31) | 4(1.25) | 4(1.26) | 0(0.96) |
Simulation study
Application to DE analysis
The calculated gene and isoform expression from NLDMseq can be used in the downstream analysis of RNA-seq, such as DE gene and isoform detection. We apply NLDMseq to MAQC dataset to show the usefulness of our method in DE analysis and compare it with the other alternatives. The 279 qRT-PCR validated DE genes with high certainty are taken as the golden standard to verify DE analysis using our method. To reach this goal, we combine NLDMseq with PPLR [25] approach, which is a hierarchical Bayesian model previously applied to microarray analysis and has been implemented in the R package puma [26]. PPLR combines the biologically replicated expression measurements and considers measurement uncertainty associated with these estimates. This approach can be applied to the expression estimated from both microarray and RNA-seq. As we introduced in the section of ‘Methods’, NLDMseq is able to produce a level of measurement uncertainty associated with expression estimate. We propagate this measurement uncertainty into PPLR to perform DE analysis. Details of the usage of PPLR can be found in the documentation of puma. Cufflinks, RSEM and MMSEQ are combined with the embedded DE analysis approaches, Cuffdiff [27], EBSeq [28], and MMDiff [29], respectively.
Comparison on computation time
The expression computation time (in minutes) of various methods for the NC phenotype of the SELC dataset. All methods use parallel computing involving four threads. We show the running time for expression computation and exclude the time for read alignment. Computation time is obtained on a 3.2 GHz quad-core Intel machine with 16G RAM
Method | Cufflinks | RSEM | MMSEQ | NLDMseq |
---|---|---|---|---|
Computation time | 70 | 44 | 50 | 26 |
Conclusion
In this contribution, we have presented a latent variable model, NLDMseq, to accurately estimate gene and isoform expression from RNA-Seq data given a known annotation. NLDMseq handles the two major difficulties in expression calculation in RNA-Seq data analysis, read mapping ambiguity and non-uniformity of read distribution, by adopting a latent model to discover the hidden pattern of read sequencing and modeling the isoform- and exon-specific read sequencing rates. Unlike many latent variable models which simulate the sequencing of each individual read, NLDMseq models normalized read counts for each exon without dealing with nucleotide constitution of each read. This significantly reduces the computation load of the model. We have used the replicate information of multiple lanes for a single library run to achieve read sequencing rate for each isoform-specific exon which shows the particular sequencing bias for this exon. We have constructed junction-exons to deal with reads mapped to the region of exon junction. For accommodation of noise reads and reads mapped to undiscovered isoforms, we have considered a special isoform which contains all the exons of a gene. We have found that the fraction of this special isoform is usually lower than that of known isoforms indicating that most reads are useful in calculation of expression for known isoforms and there are few reads which are due to noise and undiscovered isoforms.
We have employed simulation and qRT-PCR validated real data to verify the performance of our method in terms of accuracy in the calculation of gene and isoform expression. Results have shown that NLDMseq is comparable to and in some cases outperforms the other competing methods. We have also applied NLDMseq to DE analysis by combining a previously designed DE analysis method, PPLR, which has been successfully applied to microarray analysis. This DE analysis approach was tested on a well-studied qRT-PCR validated benchmark. Results have shown that our approach have presented competitive performance compared with other popular methods. We also used a dataset to examine the running time of our approach. The comparison results indicate that our approach is computationally efficient compared with the other alternatives.
One advantage of our method is that isoform- and exon-specific sequencing rate can be estimated to account for different sequencing bias for each isoform-specific exon. However, we find the distribution of the observed read counts is over-dispersed compared with the obtained sequencing rates as shown in subplots g and h in Fig. 3. To address this problem, a Dirichlet prior can be put over the sequencing rates for isoform-specific exons in the future work. We expect this Dirichlet-Multinomial model of observed counts would further improve the performance of our model.
We have used an example of detecting DE genes to show the application of NLDMseq to DE analysis. As for the detection of differential alternative splicing, NLDMseq can also be combined with available DE analysis method, such as PPLR, to find DE isoforms by using the isoform expression calculated in Eq. (7). Moreover, by using the posterior distribution of the percentage of isoforms, {θ _{ l }}, obtained from the variational EM estimation, it is also possible to detect differential isoform usage, where multiple isoforms of a single gene are expressed, but at different proportions between two groups of samples. This also provides a possibility of identifying alternative splicing.
Declarations
Acknowledgements
We acknowledge the support from NSFC (61170152), Qing Lan Project and the Fundamental Research Funds for the Central University (CXZZ11_0217).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.
Authors’ Affiliations
References
- Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009; 10:57–63.View ArticlePubMedPubMed CentralGoogle Scholar
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008; 5:621–8.View ArticlePubMedGoogle Scholar
- Jiang H, Wong WH. Statistical inferences for isoform expression in RNA-Seq. Bioinformatics. 2009; 25:1026–32.View ArticlePubMedPubMed CentralGoogle Scholar
- Turro E, Su SY, Conçalves Â, Coin LJ, Richardson S, Lewin A. Haplotype and isoform specific expression estimation using multi-mapping RNA-Seq reads. Genome Biol. 2011; 12:3.View ArticleGoogle Scholar
- Wu Z, Wang X, Zhang X. Using non-uniform read distribution models to improve isoform expression inference in RNA-Seq. Bioinformatics. 2011; 27:502–8.View ArticlePubMedGoogle Scholar
- Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN. RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics. 2010; 26:493–500.View ArticlePubMedGoogle Scholar
- Katz H, Wang ET, Airoldi EM, Burge CB. Analysis and design of RNA sequencing experiments for identifying isoform regulation. Nat Methods. 2010; 12:1009–15.View ArticleGoogle Scholar
- Glaus P, Honkela A, Rattray M. Identifying differentially expressed transcripts from RNA-Seq data with biological variation. Bioinformatics. 2012; 28:1721–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Li W, Jiang T. Transcriptome assembly and isoform expression level estimation from biased RNA-Seq reads. Bioinformatics. 2012; 28:2914–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Li L, Jiang H, Wong WH. Modeling non-uniformity in short-read rates in RNA-Seq data. Genome Biol. 2010; 11:50.View ArticleGoogle Scholar
- Srivastava S, Chen L. A two-parameter generalized Poisson model to improve the analysis of RNA-Seq data. Nucleic Acids Res. 2010; 38:170.View ArticleGoogle Scholar
- Roberts A, Trapnell C, Donaghey J, Rinn JL, Pachter L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. 2011; 12:22.View ArticleGoogle Scholar
- Jones DC, Ruzzo WL, Peng X, Katze MG. A new approach to bias correction in RNA-Seq. Bioinformatics. 2012; 28:921–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and abundance estimation from rna-seq reveals thousands of new transcripts and switching among isoforms. Nat Biotechnol. 2010; 28:511–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Suo C, Calza S, Salim A, Pawitan Y. Joint estimation of isoform expression and isoform-specific read distribution using multi-sample RNA-Seq data. Bioinformatics. 2014; 30:506–13.View ArticlePubMedGoogle Scholar
- Blei DM, Ng AY, Jordan MI. Latent Dirichlet allocation. J Mach Learn Res. 2003; 3:993–1022.Google Scholar
- Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012; 9:357–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Consortium M. The MicroArray Quality Control (MAQC) project shows inter- and intraplatform reproducibility of gene expression measurements. Nat Biotechnol. 2006; 24:1151–61.View ArticleGoogle Scholar
- Bemmo A, Benovoy D, Kwan T, Gaffney DJ, Jensen RV, Majewski J. Gene expression and isoform variation analysis using Affymetrix Exon Arrays. BMC Genomics. 2008; 9:529.View ArticlePubMedPubMed CentralGoogle Scholar
- Bullard J, Purdom E, Hansen K, Dudoit S. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinf. 2010; 11:94.View ArticleGoogle Scholar
- Rapaport F, Khanin R, Liang Y, Pirun M, Krek A, Zumbo P, et al. Comprehensive evaluation of differential gene expression analysis methods for RNA-Seq data. Genome Biol. 2013; 14:95.View ArticleGoogle Scholar
- Beane J, Vick J, Schembri F, Anderlind C, Gower A, Campbell J, et al. Characterizing the impact of smoking and lung cancer on the airway transcriptome using RNA-seq. Cancer Prev. 2011; 4:803–17.View ArticleGoogle Scholar
- Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, et al. Alternative isoform regulation in human tissue transcriptomes. Nature. 2008; 456:470–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Kim H, Bi Y, Pal S, Gupta R, Davuluri RV. IsoformEx: isoform level gene expression estimation using weighted non-negative least squares from mRNA-Seq data. BMC Bioinf. 2011; 12:305.View ArticleGoogle Scholar
- Liu X, Milo M, Lawrence ND, Rattray M. Probe-level measurement error improves accuracy in detecting differential gene expression. Bioinformatics. 2006; 22:2107–13.View ArticlePubMedGoogle Scholar
- Liu X, Gao Z, Zhang L, Rattray M. puma 3.0: improved uncertainty propagation methods for gene and transcript expression analysis. BMC Bioinf. 2013; 14:39.View ArticleGoogle 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. 2012; 31:46–53.View ArticlePubMedGoogle Scholar
- Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BMG, et al. EBSeq: An empirical bayes hierarchical model for inference in RNA-Seq experiments. Bioinformatics. 2013; 29:1035–43.View ArticlePubMedPubMed CentralGoogle Scholar
- Turro E, Astle WJ, Tavaré S. Flexible analysis of RNA-seq data using mixed effects models. Bioinformatics. 2014; 30:180–8.View ArticlePubMedGoogle Scholar