Source-linked AI summary

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

Peter Glaus, Antti Honkela, Magnus Rattray

arXiv:1109.0863v2q-bio.GNstat.AP

TL;DR

Transcript-level RNA-seq analysis must account for ambiguous read assignments, finite sampling, and biological variation when estimating expression and differential expression. BitSeq uses Bayesian posterior sampling for transcript expression and propagates that uncertainty into replicate-aware differential-expression analysis. The paper presents this approach as providing full posterior distributions and more reliable differential-expression estimates with biological replicates.

  • Problem

    RNA-seq transcript-expression and differential-expression analysis must account for ambiguous shared sequences, read-level uncertainty, and biological variation.

  • Method

    BitSeq uses a Bayesian generative model and MCMC sampling to estimate transcript-expression posteriors, then models biological variance across replicates for differential-expression analysis.

  • Results

    BitSeq provides full posterior distributions for transcript expression and produces more reliable differential-expression estimates with biological replicates.

  • Takeaways & Limitations

    Posterior expression distributions can support downstream analysis, while biological replication is important for accurate differential-expression assessment.

  • Takeaways & Limitations

    The model assumes reads are uniformly sequenced from molecule fragments, and its biological-variance hyperparameters are shared by transcripts with similar expression levels.

Abstract

from arXiv · show

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++.

1 Introduction

RNA-seq enables transcript-level expression analysis, but differential-expression analysis must handle ambiguous read origins and both technical and biological variation. BitSeq addresses these challenges with integrated probabilistic models and Bayesian inference.

  • Motivation: RNA-seq supports transcriptome analysis and expression estimation, but effective analysis remains challenging.Its applications include unknown transcript discovery, gene-expression estimation, and single-nucleotide polymorphism analysis.
  • Motivation: Differential-expression analysis is essential for comparing transcript or gene expression across conditions and studying biological differences.Applications include tissue comparisons, developmental changes, and microRNA target prediction.
  • Contribution: BitSeq models read generation and biological replication jointly to account for technical uncertainty and biological variation in differential-expression analysis.It retains posterior samples rather than point estimates, propagating uncertainty from read summarization into downstream analysis.
  • Background: Transcript-level analysis is supported because reads aligned to a gene approximately reflect transcript abundance, while shared exons create ambiguity that must be resolved probabilistically.Summing transcript expression levels can provide more precise gene-expression estimates than summing reads over all exons.
  • Motivation: Shared transcript sequences make read origins undecidable, requiring probabilistic transcript-expression estimation.Earlier EM approaches provided point estimates without measuring expression uncertainty, motivating Bayesian posterior inference.
  • Contribution: BitSeq uses an MCMC-based Bayesian model to infer transcript expression and expression changes between conditions, including reads that map to multiple genes.Its model estimates overall relative transcript expression across the transcriptome rather than restricting comparisons to transcripts within one gene.

2 Methods

BitSeq estimates transcript expression probabilistically from RNA-seq reads and assesses differential expression across replicates by propagating technical uncertainty while modelling biological variance.

  • Stage 1: Transcript expression estimation: The pipeline first estimates transcript expression from single-end or paired-end reads by sampling the posterior distribution over transcript expression levels.Expression may be summarized by posterior means when only point estimates are needed.
  • Stage 2: Differential expression: For differential expression, posterior samples from all replicates are treated as pseudo-data, condition means are inferred, and Bayesian model averaging combines evidence across samples.The approach models technical variance and regularizes variance estimates when replicate numbers are low.
  • Stage 1: Transcript expression estimation: Its generative model represents reads as conditionally independent observations linked to transcript abundance, transcript origin, alignment likelihood, and sequencing noise.A noise parameter accounts for unmapped and low-quality reads, while aligned reads may support paired-end fragment-length modelling.
  • Stage 1: Transcript expression estimation: MCMC inference uses collapsed Gibbs sampling after marginalizing transcript-abundance parameters, producing posterior samples that quantify transcript-expression uncertainty.Multiple chains are sampled and convergence is monitored with the bR statistic.
  • Stage 2: Differential expression: Biological variation is modelled hierarchically on log-expression levels, with per-transcript precision given an expression-level-dependent Gamma prior.Condition-specific replicate expression is represented with normal distributions and replicate-specific normalization constants.
  • Stage 2: Differential expression: Transcripts are ranked using the Probability of Positive Log-Ratio, the fraction of posterior samples in which one condition’s mean expression exceeds the other’s.This yields rankings for probable up-regulated and down-regulated transcripts.

3 Results and Discussion

BitSeq was evaluated for transcript expression estimation and differential expression using MAQC, Xu et al., and simulated RNA-seq datasets. Results highlight posterior uncertainty, read-distribution correction, biological-replicate modelling, and strong overall DE performance.

  • Evaluation datasets: The evaluation used MAQC data with qRT-PCR measurements, Xu et al. data with technical and biological replicates, and simulated data with known differential expression.The MAQC dataset included 893 matching transcripts for qRT-PCR comparison, while the Xu et al. dataset supported analysis of biological replication.
  • Expression level inference: Posterior transcript-expression samples revealed negative correlation between transcripts sharing substantial sequence, demonstrating uncertainty from read-mapping ambiguity.The correlated transcripts were uc010oho.1 and uc010ohp.1; another transcript pair showed no visible correlation.
  • Expression estimation accuracy: BitSeq’s expression estimates were slightly more accurate than three alternatives on the MAQC qRT-PCR comparison, while bias correction improved performance over the uniform read model.BitSeq and Cufflinks used the same read-distribution bias correction; Cufflinks version 1.3.0 performed worse than version 0.9.3 in the reported comparison.
  • Expression estimation accuracy: Effective-length normalization produced BitSeq’s major accuracy improvement, and effective-length-adjusted measures such as RPKM were more suitable for DE analysis than normalized read counts.The comparison used bias-corrected effective length when converting relative fragment expression for qRT-PCR assessment.
  • Differential expression analysis: Accounting for biological variance reduced the probability of positive log ratio from 0.995 under naive replicate combination to 0.836 when condition means were inferred.The example used two biological replicates for each of two conditions and showed that naive combination produced overconfident conclusions.
  • Differential expression analysis: BitSeq was the most accurate overall method in simulated transcript-level DE analysis, followed by baySeq, edgeR, DESeq, and Cuffdiff.The trend was especially clear at lower expression levels, whereas DESeq and edgeR had slightly higher true-positive rates for the highest-expressed transcripts at larger false-positive rates.

4 Conclusion

The paper presents Bayesian models for transcript expression and differential expression that represent uncertainty in RNA-seq data. BitSeq provides posterior distributions, models biological variation across replicates, and produces more reliable DE estimates when biological replication is available.

  • Conclusion: BitSeq uses a Bayesian generative model of transcriptome sequencing and samples posterior transcript-expression distributions.The model incorporates read and alignment quality, non-uniform read distributions, and paired-end fragment-length distributions.
  • Conclusion: Posterior transcript-expression distributions support downstream differential expression analysis while propagating uncertainty from expression inference.The second stage assesses differences between conditions using inferred expression distributions.
  • Conclusion: With biological replicates, BitSeq models intrinsic noise and variation, producing more reliable condition estimates and fewer false differential expression calls.The paper emphasizes that accurate differential-expression assessment requires experimental designs containing biological replication.

A Methods

BitSeq models read-generation and biological-replication uncertainty probabilistically, estimating transcript expression with MCMC and differential expression from propagated sample-level uncertainty.

  • Read-generation model: Paired-end read probabilities combine strand, position, insert-length, and sequencing-reporting components for each candidate transcript.Without strand-specific protocols, strand probabilities are equal; fragment lengths are assumed log-normal.
  • Read-generation model: Position probabilities can incorporate transcript-specific bias, while uniform sampling reduces to a length-dependent distribution.The uniform form is P(p|m) = 1/(lm −lr +1).
  • Read-generation model: Sequence probabilities use base-quality scores and mismatches, converting Phred scores into per-base error probabilities.The error probability is perr,i = 10−Phredi/10.
  • Bias correction: Bias correction models sequence-specific and positional effects using variable-length Markov models and empirical read frequencies.Sequence-bias parameters are estimated from singly aligned reads, weighted by an initial expression estimate; positional probabilities use transcript-length and relative-position bins.
  • Effective length: Bias-corrected effective lengths account for fragment-length and read-distribution weights and can substantially improve expression-estimation accuracy.Effective-length normalization adjusts relative fragment expression for transcript-level reporting measures such as RPKM.

A.4.2 Hyper-parameter estimation across transcript groups

BitSeq estimates expression-dependent hyperparameters from transcript groups with similar expression, then smooths group estimates across empirical mean expression levels.

  • Group-based estimation: Hyperparameter estimation groups transcripts with similar expression and estimates Gamma-distribution parameters from their pseudo-data samples.The group mean sets µ0, while λ0 is fixed at 2.0.
  • Posterior inference: The hyperparameter posterior is sampled with Metropolis-Hastings random-walk MCMC because it lacks a standard form.The sampled parameters are α and β, with β serving as the Gamma-distribution rate.
  • Expression dependence: Lowess smoothing across transcript groups estimates each transcript’s mean hyperparameter from its empirical mean expression across conditions.This provides a transcript-specific hyperparameter estimate based on expression level.
  • Illustration: Figure 7 visualizes posterior expression densities for transcript pairs from the same gene using MCMC samples, alongside their transcript sequence profile.The sequencing data come from an miRNA-155 study.

B Results

For transcripts sharing substantial sequence, BitSeq’s posterior expression distributions can be correlated, including negative correlation when reads cannot be assigned confidently between isoforms.

  • Posterior dependence: Transcripts sharing a large proportion of sequence exhibit correlation in their expression posterior distributions.The example examines three transcripts from the only UCSC Known Gene annotation for gene Q6ZMZ0.
  • Posterior dependence: The expression samples of uc001bwm.3 and uc010ohp.1 are negatively correlated because some reads cannot be assigned to a specific transcript.The posterior distribution represents this assignment ambiguity rather than selecting one transcript deterministically.

B.2 Read distribution bias correction

Bias correction improves transcript-expression estimation primarily through effective-length normalization that accounts for fragment-length distributions and read-distribution bias.

  • Evaluation: Four expression-estimation methods were compared using Pearson R2 correlations between 893 transcript estimates and TaqMan qRT-PCR measurements.Methods were evaluated on three technical replicates separately and on their averaged expression levels.
  • Effective-length normalization: BitSeq’s bias-correction benefit comes mainly from normalizing relative fragment expression with effective transcript length.This normalization converts fragment expression into transcript expression or RPKM.
  • Normalization variants: BitSeq was evaluated with transcript length, fragment-length-adjusted effective length, and bias-weighted effective length normalization.The bias-weighted version accounts for fragment-length distribution and read-distribution bias.
  • Results and scope: The bias-weighted effective-length approach provided the best results on this dataset.It combines fragment-length distribution with read-distribution bias weights, although further investigation is limited by few datasets with known expression for paired-end reads.

B.3 Assessing transcript expression inference using simulated data

BitSeq expression estimates were evaluated against known simulated transcript abundances and compared with Cufflinks, RSEM, and MMSEQ. Its overall accuracy was close to the best-performing methods, while within-gene performance favored BitSeq under stringent read-count criteria.

  • 10M simulated paired-end reads were generated from known expression estimates to assess BitSeq against Cufflinks, RSEM, and MMSEQ.The simulation used a fragment-size distribution inferred from SRA012427 and UCSC NCBI37/hg19 knownGene transcripts.
  • RSEM achieved the highest overall correlation with ground truth at R2 = 0.998, followed by MMSEQ at R2 = 0.997 and BitSeq at R2 = 0.995.The coefficient covered 46,841 transcripts with at least one generated read.
  • Cufflinks performed substantially worse than the other methods, with R2 = 0.307 for its latest version and R2 = 0.784 for version 0.9.3.
  • BitSeq performed best for the narrow within-gene transcript set, while RSEM and MMSEQ had comparable results.The narrow criterion included transcripts whose genes had at least 100 ground-truth reads.
  • Under the less stringent within-gene criterion, BitSeq retained very good ground-truth correlation while the other methods deteriorated.The two criteria used genes with at least 10 reads and at least 100 reads in the ground truth.

B.4 Biological variance of RNA-seq data

RNA-seq expression estimates show both technical and biological variation across replicates. Technical uncertainty decreases with expression, whereas biological variance becomes relatively more important at high expression and mapping ambiguity and limited depth remain informative at low expression.

  • The RNA-seq dataset contained technical and biological replicates for each condition, and significant biological variation was observed within conditions.
  • Higher expression reduced expression-estimation variation, as expected.The comparison examined standard deviation across expression levels.
  • At high expression levels, technical-replicate variance approached within-experiment variance, while biological variance became relatively more significant.
  • MCMC-estimated intrinsic variance from mapping ambiguity and limited read depth provides useful information for assessing replicate variance in the low-expression regime.

B.5 Assessing DE performance with simulated data

BitSeq’s differential-expression analysis was evaluated against four alternatives using simulated data and replicate-aware comparisons. Performance generally improved with larger fold changes, with BitSeq and baySeq usually strongest except at the smallest fold change.

  • ROC curves were averaged over 5 runs with standard-deviation error bars for transcripts averaging at least one ground-truth read and fold changes from 1.5 to 3.5.
  • The evaluation compared BitSeq with Cuffdiff, DESeq, edgeR, and baySeq using expression estimates converted into read counts for the gene-level methods.The comparison used default parameter settings for each method.
  • Larger fold changes clearly improved method performance, and BitSeq and baySeq generally outperformed the other methods.At fold change 1.5, baySeq had the lowest true positive rate, while edgeR and DESeq outperformed BitSeq in half of the spectrum.
  • Cuffdiff performed worst among the five algorithms in all differential-expression experiments.Returning to Cufflinks version 0.9.3 did not improve its differential-expression results.
Loading 1109.0863v2…