Fast and accurate approximate inference of transcript expression from RNA-seq data

Abstract

Motivation: Assigning RNA-seq reads to their transcript of origin is a fundamental task in transcript expression estimation. Where ambiguities in assignments exist due to transcripts sharing sequence, e.g. alternative isoforms or alleles, the problem can be solved through probabilistic inference. Bayesian methods have been shown to provide accurate transcript abundance estimates compared with competing methods. However, exact Bayesian inference is intractable and approximate methods such as Markov chain Monte Carlo and Variational Bayes (VB) are typically used. While providing a high degree of accuracy and modelling flexibility, standard implementations can be prohibitively slow for large datasets and complex transcriptome annotations.

Results: We propose a novel approximate inference scheme based on VB and apply it to an existing model of transcript expression inference from RNA-seq data. Recent advances in VB algorithmics are used to improve the convergence of the algorithm beyond the standard Variational Bayes Expectation Maximization algorithm. We apply our algorithm to simulated and biological datasets, demonstrating a significant increase in speed with only very small loss in accuracy of expression level estimation. We carry out a comparative study against seven popular alternative methods and demonstrate that our new algorithm provides excellent accuracy and inter-replicate consistency while remaining competitive in computation time.

Availability and implementation: The methods were implemented in R and C++, and are available as part of the BitSeq project at github.com/BitSeq. The method is also available through the BitSeq Bioconductor package. The source code to reproduce all simulation results can be accessed via github.com/BitSeq/BitSeqVB_benchmarking.

Fig. 1.

Graphical model of the RNA-seq mixture problem. Given a known Transcriptome $\mathbf{T}$ and some observed reads $\mathbf{R}$, the inference problem is for $\theta$ through the latent variables $\mathbf{Z}$

Fig. 2.

Ranking of methods for five replicates of simulated RNA-seq reads. WGE-Inter: inter-replicate consistency of within gene estimates, WGE-True: within gene estimates compared with the true values and Theta: estimated relative transcript expression compared with the true values. Scores have been normalized to unity per dataset. Alternative normalizations are available in supplementary material (Supplementary Fig. S9)

Fig. 3.

Run-time in hours (log-scale) for four synthetic data samples with five replicates per sample. The total number of simulated reads is shown in parenthesis

Fig. 4.

Five ENCODE pairs of replicates. (a) Alignment rates for transcriptome mapping (Bowtie1 and Bowtie2), genome mapping (Tophat 2.0.9 with Bowtie1 and Bowtie2), k-mers mapping (Sailfish) and pseudo-alignments (Kallisto). (b) Ranking of methods in terms of the Mean Absolute Error. (c) Run-time in hours (log-scale) with 24.6 M (mapped) reads per sample

Fig. 5.

Scatterplots of within gene estimates for one pair of replicates (SRR307907 and SRR307908) from the ENCODE data. The blue color corresponds to a smoothed color density representation of the scatterplot

Fig. 6.

A comparison of the first two moments of the approximate posterior expression in counts per transcript: (a) posterior mean (R² correlation is 0.999) (b) posterior standard deviation: the VB method significantly under-estimates the posterior variance (${\sigma }^2$). (c), (d) posterior mean-variance relationship in MCMC and VB respectively. Shading represents the number of transcripts in each region

Fig. 7.

Convergence comparison of Collapsed MCMC with standard VB algorithm and VB with Fletcher-Reeves conjugate gradient optimization. Expression estimates obtained by very long run of MCMC are used as a ground truth and average root mean square error over 10 runs was calculated, two standard deviations are used as error bars. The VB methods with several randomized initial conditions showed negligible differences in convergence