Improving RNA-Seq expression estimates by correcting for fragment bias - PubMed (original) (raw)
Improving RNA-Seq expression estimates by correcting for fragment bias
Adam Roberts et al. Genome Biol. 2011.
Abstract
The biochemistry of RNA-Seq library preparation results in cDNA fragments that are not uniformly distributed within the transcripts they represent. This non-uniformity must be accounted for when estimating expression levels, and we show how to perform the needed corrections using a likelihood based approach. We find improvements in expression estimates as measured by correlation with independently performed qRT-PCR and show that correction of bias leads to improved replicability of results across libraries and sequencing technologies.
Figures
Figure 1
Overview of a typical RNA-Seq experiment. RNA is initially fragmented (1) followed by first-strand synthesis priming (2), which selects the 3' fragment end (in transcript orientation), to make single stranded cDNA. Double stranded cDNA created during second-strand synthesis (3), which selects the 5' fragment end, is then size selected (4) resulting in fragments suitable for sequencing (5). Sequenced reads are mapped to opposite strands of the genome (6), and in the case of known transcript or fragment strandedness, the read alignments reveal the 5' and 3' ends of the sequenced fragment (see Supplementary methods in Additional file 3). All arrows are directed 5' to 3' in transcript orientation.
Figure 2
Nucleotide distribution surrounding fragment ends and calculation of bias weights. (a) Sequence logos showing the distribution of nucleotides in a 23 bp window surrounding the ends of fragments from an experiment primed with 'not not so random' (NNSR) hexamers [11]. The 3' end sequences are complemented (but not reversed) to show the sequence of the primer during first-strand synthesis (see Figure 1). The offset is calculated so that zero is the 'first' base of the end sequence and only non-negative values are internal to the fragment. Counts were taken only from transcripts mapping to single-isoform genes. (b) Sequence logo showing normalized nucleotide frequencies after reweighting by initial (not bias corrected) FPKM in order to account for differences in abundance. (c) The background distribution for the yeast transcriptome, assuming uniform expression of all single-isoform genes. The difference in 5' and 3' distributions are due to the ends being primed from opposite strands. Comparing (c) to (a) and (b) shows that while the bias is confounded with expression in (a), the abundance normalization reveals the true bias to extend from 5 bp upstream to 5 bp downstream of the fragment end. Taking the ratio of the normalized nucleotide frequencies (b) to the background (c) for the NNSR dataset gives bias weights (d), which further reveal that the bias is partially due to selection for upstream sequences similar to the strand tags, namely TCCGATCTCT in first-strand synthesis (which selects the 5' end) and TCCGATCTGA in second-strand synthesis (which selects the 3' end). Although the weights here are based on independent frequencies, we found correlations among sites in the window and take these into account in our full model to produce more informative weights (see Supplementary methods in Additional file 3). A similar figure to this for the standard Illumina Random Hexamer protocol and plots similar to (d) for all datasets in the paper can be found in Figures S1 and S2 of Additional file 1 respectively.
Figure 3
Bias correction within transcripts. An example showing the effect of bias correction on the read counts for human transcript
NM_004684
. The top panel shows raw read counts (number of 3' ends of fragments at each location), and the bottom panel shows the product of the bias parameters (total bias weight defined in the Supplementary methods in Additional file 3) at the same locations. We correctly identify bias at different positions and can therefore correct for the non-uniformity. Note that the bias parameters were learned from the entire dataset excluding reads mapped to this transcript in order to cross-validate our results. The RNA-Seq for the experiment was performed with the NSR protocol [21], which is why 3' counts were used instead of 5'.
Figure 4
Correlation between RNA-Seq and qRT-PCR. (a) Expression estimates before bias correction (tail of arrows) and after correction (points of arrows) for the SRA012427 dataset compared to qRT-PCR values for the same transcripts. Red arrows show decrease in expression after correction and blue an increase. Note that we have zoomed in on lower-expression transcripts (the majority) for clarity. (b) Distribution of log-fold change in expression after bias correction.
Figure 5
Comparison with previous methods. A comparison of our method (
Cufflinks
) with
Genominator
[8] and
mseq
[10]. The _y_-axis shows the _R_2 value for the correlation between uncorrected (green) and bias corrected (orange) RNA-Seq expression estimates and qRT-PCR for the three methods. Correlation plots for these data can be found in Figure S3 of Additional file 1.
Figure 6
Variable technical replicates. Results of correlation tests showing improvement after bias correction for technical replicates. Fraction Explained Discrepancy was calculated by dividing the change in _R_2 after bias correction by the difference of the initial _R_2 from one (a perfect correlation). Note that when two RNA-Seq datasets are compared, the correction in the legend was applied to both. The pairwise correlations of the four SRA010153 replicates versus qRT-PCR and SRA008403, respectively, were averaged for the figure. Even though the same RH priming protocol was used in both labs, the bias differs slightly (see Figure S2 of Additional file 1) between the preps, which is why our correction method was able to improve the correlation.
Figure 7
Variable library preparations. Results of correlation tests showing improvement after bias correction of datasets generated using different library prep methods, all of which are strand-specific. The first four protocols are described in [11] and the final in [21]. All datasets were compared against a control that was generated using the standard Illumina RH protocol. The first four datasets used the control from [11] with the same yeast sample. The last dataset (NSR) was compared against the HBR dataset from SRA010153 since it is also consists of single-end reads.
Figure 8
Bias in different sequence technologies. Results of correlation tests showing improvement after bias correction of datasets generated using different sequencing technologies. The Illumina dataset is SRA012427 (_x_-axes) and the SOLiD data is SOLiD4_ HBR_PE_50x25 (_y_-axes). Both used the same MAQC HBR sample. Red axes and lines denote uncorrected FPKM values and blue corrected, while purple regression lines denote a comparison between corrected and uncorrected values. Both datasets are being corrected for different biases, which causes their expression estimates to become more correlated. Note that the plot is zoomed in on the lower abundance transcripts for clarity but captures over 98% of those in the experiment.
Figure 9
Positional bias correction effect on expression. This figure shows the effective length fold change due to positional bias correction for the SRA012427 dataset. So that the parameters would be consistent for all transcripts, we have limited the analysis to transcripts with length greater than 2,433, which is the largest of the 5 length bins we use for measuring positional bias. As expected, all single isoform genes are adjusted in the same way, whereas isoform deconvolution is altered based on the difference in relative position within a transcript for a read that maps to multiple isoforms.
Figure 10
Correlation of GC content with measured bias. Panel a shows an example from human and panel b from yeast. Since the log fold change in effective length does not capture the full bias measurement for multiple isoform genes, the plots are limited to those with single isoforms. GC content appears to be correlated with our sequence bias measurements in some datasets, but not in others. GC content alone is not always a good proxy for fragment end bias.
Similar articles
- Peregrine: A rapid and unbiased method to produce strand-specific RNA-Seq libraries from small quantities of starting material.
Langevin SA, Bent ZW, Solberg OD, Curtis DJ, Lane PD, Williams KP, Schoeniger JS, Sinha A, Lane TW, Branda SS. Langevin SA, et al. RNA Biol. 2013 Apr;10(4):502-15. doi: 10.4161/rna.24284. Epub 2013 Apr 1. RNA Biol. 2013. PMID: 23558773 Free PMC article. - Library preparation methods for next-generation sequencing: tone down the bias.
van Dijk EL, Jaszczyszyn Y, Thermes C. van Dijk EL, et al. Exp Cell Res. 2014 Mar 10;322(1):12-20. doi: 10.1016/j.yexcr.2014.01.008. Epub 2014 Jan 15. Exp Cell Res. 2014. PMID: 24440557 Review. - Optimization for Sequencing and Analysis of Degraded FFPE-RNA Samples.
Levin Y, Talsania K, Tran B, Shetty J, Zhao Y, Mehta M. Levin Y, et al. J Vis Exp. 2020 Jun 8;(160):10.3791/61060. doi: 10.3791/61060. J Vis Exp. 2020. PMID: 32568231 Free PMC article. - Mixture models reveal multiple positional bias types in RNA-Seq data and lead to accurate transcript concentration estimates.
Tuerk A, Wiktorin G, Güler S. Tuerk A, et al. PLoS Comput Biol. 2017 May 15;13(5):e1005515. doi: 10.1371/journal.pcbi.1005515. eCollection 2017 May. PLoS Comput Biol. 2017. PMID: 28505151 Free PMC article. - Preparation of Single-Cell RNA-Seq Libraries for Next Generation Sequencing.
Trombetta JJ, Gennert D, Lu D, Satija R, Shalek AK, Regev A. Trombetta JJ, et al. Curr Protoc Mol Biol. 2014 Jul 1;107:4.22.1-4.22.17. doi: 10.1002/0471142727.mb0422s107. Curr Protoc Mol Biol. 2014. PMID: 24984854 Free PMC article. Review.
Cited by
- CGF therapy: bridging androgenetic alopecia observations to psoriasis treatment via IL-17 pathway.
Xiao Q, Chu W, Guo J, Gao J, Yao W, Huang M, Lu Y, Xu Q, Xu N. Xiao Q, et al. Stem Cell Res Ther. 2024 Oct 8;15(1):353. doi: 10.1186/s13287-024-03959-y. Stem Cell Res Ther. 2024. PMID: 39380104 Free PMC article. - Mechanical force promotes tissue and molecular changes in adipose tissue regeneration post-transplantation.
Ye Y, Ma J, Guo BY, Li XJ, Hu KK, Tan MJ, Zhang L. Ye Y, et al. Front Cell Dev Biol. 2024 Sep 18;12:1472575. doi: 10.3389/fcell.2024.1472575. eCollection 2024. Front Cell Dev Biol. 2024. PMID: 39359720 Free PMC article. - Enhancing RNA-seq bias mitigation with the Gaussian self-benchmarking framework: towards unbiased sequencing data.
Su Q, Long Y, Gou D, Quan J, Lian Q. Su Q, et al. BMC Genomics. 2024 Sep 30;25(1):904. doi: 10.1186/s12864-024-10814-0. BMC Genomics. 2024. PMID: 39350040 Free PMC article. - Integrated miRNA and mRNA Transcriptome Analysis Reveals Regulatory Mechanisms in the Response of Winter Brassica rapa to Drought Stress.
Ma L, Xu Y, Tao X, Fahim AM, Zhang X, Han C, Yang G, Wang W, Pu Y, Liu L, Fan T, Wu J, Sun W. Ma L, et al. Int J Mol Sci. 2024 Sep 20;25(18):10098. doi: 10.3390/ijms251810098. Int J Mol Sci. 2024. PMID: 39337583 Free PMC article. - Strongly deleterious mutations influence reproductive output and longevity in an endangered population.
Hasselgren M, Dussex N, von Seth J, Angerbjörn A, Dalén L, Norén K. Hasselgren M, et al. Nat Commun. 2024 Sep 27;15(1):8378. doi: 10.1038/s41467-024-52741-4. Nat Commun. 2024. PMID: 39333094 Free PMC article.
References
- Paşaniuc B, Zaitlen N, Halperin E. In: Research in Computational Molecular Biology. Berger B, editor. Berlin/Heidelberg: Springer; 2010. Accurate estimation of expression levels of homologous genes in RNA-seq experiments. pp. 397–409. [Lecture Notes in Computer Science, vol 6044.] - PubMed
Publication types
MeSH terms
LinkOut - more resources
Full Text Sources
Other Literature Sources