Identifying differentially expressed transcripts from RNA-seq data with biological variation

Motivation: High-throughput sequencing enables expression analysis at the level of individual transcripts. The analysis of transcriptome expression levels and differential expression estimation requires a probabilistic approach to properly account for ambiguity caused by shared exons and finite read sampling as well as the intrinsic biological variance of transcript expression. Results: We present BitSeq (Bayesian Inference of Transcripts from Sequencing data), a Bayesian approach for estimation of transcript expression level from RNA-seq experiments. Inferred relative expression is represented by Markov chain Monte Carlo (MCMC) samples from the posterior probability distribution of a generative model of the read data. We propose a novel method for differential expression analysis across replicates which propagates uncertainty from the sample-level model while modelling biological variance using an expression-level-dependent prior. We demonstrate the advantages of our method using simulated data as well as an RNA-seq dataset with technical and biological replication for both studied conditions. Availability: The implementation of the transcriptome expression estimation and differential expression analysis, BitSeq, has been written in C++.

💡 Research Summary

The paper introduces BitSeq, a Bayesian framework for transcript‑level expression estimation and differential expression (DE) analysis from RNA‑seq data that explicitly accounts for both technical uncertainty and biological variation across replicates. The authors first model the read generation process as a hierarchical probabilistic graphical model: each read rₙ is associated with a latent transcript indicator Iₙ (categorical with Dirichlet‑distributed relative abundances θ) and a noise indicator Zₙᵃᶜᵗ (Bernoulli with Beta‑distributed noise rate θᵃᶜᵗ). Alignment probabilities P(rₙ|Iₙ=m) are pre‑computed using Bowtie, preserving all multi‑mapping possibilities and optionally correcting for positional and sequence bias. Because the posterior distribution of θ is analytically intractable, the authors marginalize θ and employ a collapsed Gibbs sampler to draw samples of the transcript assignments Iₙ, from which posterior samples of θ are obtained. Multiple MCMC chains are run and convergence is monitored with the Gelman‑Rubin b̂_R statistic.

In the second stage, the posterior samples from each sample are treated as “pseudo‑data” for each replicate. For each condition c and replicate r, the log‑expression y_{crm} of transcript m is modeled as Normal(μ_{cm}+n_{cr}, 1/λ_{cm}), where μ_{cm} is the condition‑specific mean, n_{cr} a replicate‑specific offset, and λ_{cm} a precision term. The μ’s follow a hierarchical log‑normal model, and the prior on the between‑replicate variance is made expression‑level dependent via a non‑parametric regression (e.g., spline) learned from all transcripts. Because the pseudo‑data are already corrupted by technical noise, the hierarchical model only needs to capture biological variance, and inference remains conjugate and analytically tractable. Posterior distributions of condition means are combined across pseudo‑data using Bayesian model averaging, yielding the Positive Log‑ratio Probability (PPLR) for each transcript—the probability that the log‑fold change between conditions is positive. Transcripts are ranked by PPLR to produce a DE list.

The authors evaluate BitSeq on simulated data with varying expression levels and replicate numbers, and on a real RNA‑seq dataset containing both technical and biological replicates for two conditions. Compared with EM‑based methods (e.g., Cufflinks, eXpress) and simpler Bayesian approaches, BitSeq achieves lower false discovery rates, especially for low‑expressed transcripts, and provides well‑calibrated uncertainty estimates. The method’s strengths include full preservation of multi‑mapping information, explicit noise modeling, propagation of posterior uncertainty from the first to the second stage, and an expression‑dependent variance prior that stabilizes variance estimates when replicates are few. Limitations are the computational cost of MCMC sampling and the need for careful convergence diagnostics; the authors suggest future extensions such as variational inference or GPU acceleration to improve scalability.

In summary, BitSeq offers a comprehensive Bayesian solution for transcript‑level differential expression analysis that jointly models read assignment ambiguity, technical noise, and biological replication, delivering more accurate and reliable DE calls than existing pipelines. The software is implemented in C++ and made publicly available, facilitating adoption and further development by the community.