Statistically robust methylation calling for whole-transcriptome bisulfite sequencing reveals distinct methylation patterns for mouse RNAs

Abstract

Cytosine-5 RNA methylation plays an important role in several biologically and pathologically relevant processes. However, owing to methodological limitations, the transcriptome-wide distribution of this mark has remained largely unknown. We previously established RNA bisulfite sequencing as a method for the analysis of RNA cytosine-5 methylation patterns at single-base resolution. More recently, next-generation sequencing has provided opportunities to establish transcriptome-wide maps of this modification. Here, we present a computational approach that integrates tailored filtering and data-driven statistical modeling to eliminate many of the artifacts that are known to be associated with bisulfite sequencing. By using RNAs from mouse embryonic stem cells, we performed a comprehensive methylation analysis of mouse tRNAs, rRNAs, and mRNAs. Our approach identified all known methylation marks in tRNA and two previously unknown but evolutionary conserved marks in 28S rRNA. In addition, mRNAs were found to be very sparsely methylated or not methylated at all. Finally, the tRNA-specific activity of the DNMT2 methyltransferase could be resolved at single-base resolution, which provided important further validation. Our approach can be used to profile cytosine-5 RNA methylation patterns in many experimental contexts and will be important for understanding the function of cytosine-5 RNA methylation in RNA biology and in human disease.

Figure 1.

Schematic outline of whole-transcriptome bisulfite sequencing. The illustration shows key steps of library preparation (A) and data analysis (B).

Figure 2.

WTBS of mouse ES cell mRNA. (A) Frequency of bisulfite nonconversion (resp. methylation) ratios in mRNA. (B) Estimated Poisson rate λ_p/coverage as a function of coverage, as well as corresponding boxplot. (C) Bisulfite nonconversion ratio (x-axis) and −log of methylation P-values (y-axis) for 53,510 cytosines passing the 0.2 ratio threshold. Color indicates density (count of points per symbol area). (D) Cytosine position and methylation tracks for a representative mRNA transcript (NM_007984) in three replicates.

Figure 3.

Statistical analysis of WTBS data sets. (A) Bisulfite nonconversion ratio (x-axis) and −log of methylation P-values (y-axis) for 56,940 cytosines with ratio >λ_p/coverage that are common to the three wild-type replicates. Color indicates density (count of points per symbol area). (B) Number of cytosines and nonconversion ratios, as determined experimentally (left) or by simulation (right). (C) Statistical power stratified by coverage and by nonconversion ratio. (D) Histogram of coverage in a representative (Wt1L) sequencing data set. The red line indicates a coverage of 20×. (E) Logo plot for all cytosines with a ratio >0.1 and significant P-value (<0.05) in at least one sample.

Figure 4.

LC-MS/MS analysis of RNA samples from mouse ES cells that were subjected to multistep mRNA enrichment. (A) Basic outline of the mRNA enrichment protocol. (B) Relative amounts of mRNA and rRNA as determined by RNA-seq. The proportion of tRNA reads was <1% for all samples. (C) Modification analysis of m5C, hm5C, m26A, m6A, and m1A content, relative to A content.

Figure 5.

Methylation frequencies and reproducibility in mRNA, rRNA, and tRNA. (A) Proportion of nonconverted cytosines (coverage ≥ 20×). (B) Standard error (coverage ≥ 20× and nonconversion ratio ≥0.1). Values are based on triplicate wild-type data sets for mRNA and tRNA and on duplicate data sets for rRNA.

Figure 6.

Site-specific methylation analysis by whole-transcriptome bisulfite sequencing. Scatter plots show nonconverted cytosines for tRNA (A), rRNA (B), and mRNA (C) in wild-type and Dnmt2 knockout ES cells. Methylation ratios are specifically reduced for C38 of tRNA(Asp), tRNA(Gly), and tRNA(Val) in Dnmt2 knockouts. (D) Scatter plot for mRNA in wild-type and TET-deficient ES cells.