Background

MeRIP-Seq (methylated RNA immunoprecipitation followed by sequencing) is a method used to study the epigenetic regulation of gene expression by identifying methylated RNA molecules in a sample. It involves immunoprecipitation of methylated RNA using an antibody specific for methylated nucleotides, followed by high-throughput sequencing of the immunoprecipitated RNA. The resulting data provide a snapshot of the methylation status of RNA molecules, allowing researchers to investigate the role of methylation in the regulation of gene expression and other biological processes. After being collected from two or more biological conditions, MeRIP-Seq data is typically analyzed using computational tools to identify differentially methylated regions (DMRs). DMR detection can uncover the functional significance of m6A methylation and identify potential therapeutic targets for disease such as cancer.

To establish the statistical rigor of MeRIP-Seq experiments, it is important to carefully consider sample size during the study design process in order to ensure that the experiment is adequately powered to detect differentially methylated RNA (DMR) regions. However, there is no such tool available, so we developed the R package magpie, a simulation-based tool for performing power calculations on MeRIP-Seq data. magpie has two main functions:

• Power calculation: Given a MeRIP-Seq dataset, various experimental scenarios, and a selected test method, magpie can calculate the statistical power of the experiment and output the results.

• Results preservation and visualization: magpie can save the results of the power calculation as an Excel file, and it can also produce basic line plots that allow the user to visualize the results.

Installation

From GitHub:

install.packages("devtools") # if you have not installed "devtools" package
library(devtools)
install_github("https://github.com/dxd429/magpie",
    build_vignettes = TRUE
)

From Bioconductor:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("magpie")

To view the package vignette in HTML format, run the following lines in R:

library(magpie)
vignette("magpie")

Power evaluation with powerEval()

With .BAM files and an annotation file, powerEval() enables users to specify number of simulations (nsim), sample sizes (N.reps), sequencing depths (depth_factor), FDR thresholds (thres), and testing methods (Test_method).

Users should always indicate the percentage of the whole genome covered by their dataset, whether it’s complete or partial (bam_factor). Here, bam_factor allows for simulating whole-genome data using a smaller subset of it, under the assumption that DMR signals are relatively evenly distributed across chromosomes.

We demonstrate the usage of powerEval() with the example dataset from R package magpieData:

library(magpieData)
library(magpie)
### Get the example data
BAM_path <- getBAMpath()
### Call powerEval()
power.test <- powerEval(
    Input.file = c("Ctrl1.chr15.input.bam", "Ctrl2.chr15.input.bam", "Case1.chr15.input.bam", "Case2.chr15.input.bam"),
    IP.file = c("Ctrl1.chr15.ip.bam", "Ctrl2.chr15.ip.bam", "Case1.chr15.ip.bam", "Case2.chr15.ip.bam"),
    BamDir = BAM_path,
    annoDir = paste0(BAM_path, "/hg18_chr15.sqlite"),
    variable = rep(c("Ctrl", "Trt"), each = 2),
    bam_factor = 0.03,
    nsim = 10,
    N.reps = c(2, 3, 5, 7),
    depth_factor = c(1, 2),
    thres = c(0.01, 0.05, 0.1),
    Test_method = "TRESS" ## TRESS or exomePeak2
)

To use your own data, replace Input.file and IP.file with the names of your .BAM files, and set BamDir and annoDir to the file paths for your data. Make sure that the order of the variable corresponds to the order of your data files.

Power evaluation with quickPower()

Another power calculation function is quickPower(). Unlike powerEval() which often takes a while to run, quickPower() produces results in seconds, by directly extracting results from three pre-evaluated datasets on GEO:

• GSE46705: Human HeLa cell line: Two replicates of wild type (WT) and two replicates of knockdown (KD) of complex METTL3.

• GSE55575: Mouse embryonic fibroblasts: Two replicates of wild type (WT) and four replicates of knockdown (KD) of WTAP.

• GSE94613: Human leukemia cell line: Four replicates of wild type (WT) and eight replicates of knockdown (KD) of complex METTL3.

Here, we use quickPower() to get power evaluation results of GSE46705, tested by TRESS:

library(magpie)
power.test <- quickPower(dataset = "GSE46705")

Save results to .xlsx files

magpie provides two functions: writeToxlsx() and writeToxlsx_strata(), for saving results to a formatted .xlsx file. writeToxlsx() allows you to save power measurements for different sample sizes, FDR thresholds, and sequencing depths, while writeToxlsx_strata() writes out the results stratified by input expression level and for different sample sizes. Note that we only evaluate power under the original sequencing depth (depth_factor = 1) and FDR threshold of 0.05 (thres = 1).

Here we save the output from quickPower() to corresponding .xlsx files:

### write out .xlsx
writeToxlsx(power.test, file = "test_TRESS.xlsx")

### write out stratified results
writeToxlsx_strata(power.test, file = "test_strata_TRESS.xlsx")

The generated .xlsx files are formatted as follows:

Generate figures

As mentioned before, magpie also provides four plotting functions: plotRes(), plotAll(), plotStrata(), and plotAll_Strata(), for figure generating. In general, plotRes() and plotStrata() produce individual plots, while plotAll() and `plotAll_Strata() produce four plots in a 2 x 2 panel.

Again, we demonstrate these four functions with the previous output from `quickPower()``:

### plot FDR under sequencing depth 1x
plotRes(power.test, depth_factor = 1, value_option = "FDR")

### plot all in a panel under sequencing depth 1x
plotAll(power.test, depth_factor = 1)

### plot a FDR strata result
plotStrata(power.test, value_option = "FDR")

### plot all strata results in a panel
plotAll_Strata(power.test)