magpie Package User’s Guide
24 October 2023
Abstract
This package, magpie, aims to perform statistical power analysis for differential RNA methylation calling, using MeRIP-Seq data. It takes real MeRIP-Seq data as input for parameter estimation, allows for options of setting various sample sizes, sequencing depths, and testing methods, and calculates FDR, FDC, power, and precision as evaluation metrics. It also offers functions to save results into .xlsx files and produce basic line plots.
Package
magpie 1.2.0
Contents
1 Introduction
1.1 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.
1.2 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")2 Quickstart
To get started with magpie, users are expected to provide MeRIP-Seq data for parameter estimation. This includes paired input and IP .BAM files for each experimental condition, along with the matching annotation file in the *.sqlite format. Since .BAM files are generally large, it is also encouraged to use data of only one or a few chromosomes.
magpie offers the function quickPower() for users who need quick power calculation results without their own simulations. This function takes few seconds to run and uses pre-calculated results from three publicly available MeRIP-seq datasets. The details of these datasets can be found in 4.2). To use this function and save and plot the results:
library(magpie)
power.test <- quickPower(dataset = "GSE46705") # Options are 'GSE46705', 'GSE55575', and 'GSE94613'.Save results into .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")Produce basic line plots for visulization:
### 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)
3 Input data
In order to best mimic the characteristics of real data, magpie expects users to provide real MeRIP-Seq datasets, if not directly extracting our pre-calculated results.
magpie requires paired input and IP .BAM files for each replicate of both conditions. Example file names are as follows: “Ctrl1.input.bam” & Ctrl1.ip.bam“,”Ctrl2.input.bam" & Ctrl2.ip.bam“,”Case1.input.bam" & Case1.ip.bam“,”Case2.input.bam" & Case2.ip.bam", for a 2 controls vs 2 cases DMR calling study. Note that names of provided .BAM files should be ordered carefully so that samples from different conditions will not be mistreated.
For illustration purpose, we include a sample dataset (GSE46705) from a study investigating how METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation in our experimental data package magpieData on GitHub, which can be installed with:
install.packages("devtools") # if you have not installed "devtools" package
library(devtools)
install_github("https://github.com/dxd429/magpieData",
build_vignettes = TRUE
)The data package contains four .BAM files of two wild type (WT) replicates, four .BAM files of two knockdown of complex METTL3 replicates only on chromosome 15, and one genome annotation file.
In terms of the annotation file, it should match the reference genome version when obtaining the .BAM files, and be saved in the format of *.sqlite, which can be created with R function makeTxDbFromUCSC() from Bioconductor package GenomicFeatures:
## Use "makeTxDbFromUCSC" function to create an annotation file of hg18
library(GenomicFeatures)
hg18 <- makeTxDbFromUCSC(genome = "hg18", tablename = "knownGene")
saveDb(hg18, file = "hg18.sqlite")4 Power calculation
magpie offers two functions for power evaluation: powerEval() and quickPower(). If users prefer to quickly examine the power evaluation results using the built-in datasets, they can use the quickPower() function. Otherwise, users can use the powerEval() function to perform power evaluation on their own data and customize the simulation settings. The output of these two functions is a list of power measurements under various experimental settings, such as FDR, Precision, and Statistical power.
4.1 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.
4.2 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")5 Results preservation and visualization
Magpie computes power evaluation metrics for each experimental scenario defined in the function arguments, the results it generates include:
FDR: The ratio of number of false positives to the number of positive discoveries.
FDC: The ratio of number of false positives to the number of true positives.
Power: Statistical power.
Precision: The ratio of number of true positives to the number of positive discoveries.
Once users have obtained a list of power measurements using either powerEval() or quickPower(), they can use functions in magpie to save the results as a formatted .xlsx file and generate basic line plots.
5.1 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:
| FDR | N.rep | 0.01 | 0.05 | 0.1 |
|---|---|---|---|---|
| 2 | 0.36 | 0.48 | 0.57 | |
| 3 | 0.14 | 0.27 | 0.38 | |
| 5 | 0.06 | 0.13 | 0.21 | |
| 7 | 0.04 | 0.11 | 0.17 |
| FDR | N.rep | (0, 27.68] | (27.68, 54.3] | (54.3, 92.64] | (92.64, Inf] |
|---|---|---|---|---|---|
| 2 | 0.41 | 0.46 | 0.47 | 0.58 | |
| 3 | 0.23 | 0.30 | 0.27 | 0.28 | |
| 5 | 0.15 | 0.16 | 0.11 | 0.12 | |
| 7 | 0.12 | 0.14 | 0.08 | 0.11 |
5.2 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)
Session info
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] magpie_1.2.0 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.1.3 bitops_1.0-7
## [3] biomaRt_2.58.0 rlang_1.1.1
## [5] magrittr_2.0.3 matrixStats_1.0.0
## [7] compiler_4.3.1 RSQLite_2.3.1
## [9] GenomicFeatures_1.54.0 systemfonts_1.0.5
## [11] reshape2_1.4.4 png_0.1-8
## [13] vctrs_0.6.4 rvest_1.0.3
## [15] stringr_1.5.0 pkgconfig_2.0.3
## [17] crayon_1.5.2 fastmap_1.1.1
## [19] magick_2.8.1 dbplyr_2.3.4
## [21] XVector_0.42.0 utf8_1.2.4
## [23] Rsamtools_2.18.0 rmarkdown_2.25
## [25] bit_4.0.5 xfun_0.40
## [27] zlibbioc_1.48.0 cachem_1.0.8
## [29] GenomeInfoDb_1.38.0 jsonlite_1.8.7
## [31] progress_1.2.2 blob_1.2.4
## [33] highr_0.10 DelayedArray_0.28.0
## [35] BiocParallel_1.36.0 parallel_4.3.1
## [37] prettyunits_1.2.0 R6_2.5.1
## [39] bslib_0.5.1 stringi_1.7.12
## [41] RColorBrewer_1.1-3 rtracklayer_1.62.0
## [43] GenomicRanges_1.54.0 jquerylib_0.1.4
## [45] Rcpp_1.0.11 bookdown_0.36
## [47] SummarizedExperiment_1.32.0 knitr_1.44
## [49] IRanges_2.36.0 Matrix_1.6-1.1
## [51] tidyselect_1.2.0 rstudioapi_0.15.0
## [53] abind_1.4-5 yaml_2.3.7
## [55] codetools_0.2-19 curl_5.1.0
## [57] plyr_1.8.9 lattice_0.22-5
## [59] tibble_3.2.1 Biobase_2.62.0
## [61] KEGGREST_1.42.0 evaluate_0.22
## [63] BiocFileCache_2.10.0 zip_2.3.0
## [65] xml2_1.3.5 Biostrings_2.70.0
## [67] pillar_1.9.0 BiocManager_1.30.22
## [69] TRESS_1.8.0 filelock_1.0.2
## [71] MatrixGenerics_1.14.0 stats4_4.3.1
## [73] generics_0.1.3 RCurl_1.98-1.12
## [75] S4Vectors_0.40.0 hms_1.1.3
## [77] ggplot2_3.4.4 aod_1.3.2
## [79] munsell_0.5.0 scales_1.2.1
## [81] glue_1.6.2 tools_4.3.1
## [83] BiocIO_1.12.0 webshot_0.5.5
## [85] openxlsx_4.2.5.2 locfit_1.5-9.8
## [87] GenomicAlignments_1.38.0 XML_3.99-0.14
## [89] grid_4.3.1 AnnotationDbi_1.64.0
## [91] colorspace_2.1-0 GenomeInfoDbData_1.2.11
## [93] restfulr_0.0.15 cli_3.6.1
## [95] rappdirs_0.3.3 kableExtra_1.3.4
## [97] fansi_1.0.5 viridisLite_0.4.2
## [99] S4Arrays_1.2.0 svglite_2.1.2
## [101] dplyr_1.1.3 gtable_0.3.4
## [103] DESeq2_1.42.0 sass_0.4.7
## [105] digest_0.6.33 BiocGenerics_0.48.0
## [107] SparseArray_1.2.0 rjson_0.2.21
## [109] memoise_2.0.1 htmltools_0.5.6.1
## [111] lifecycle_1.0.3 httr_1.4.7
## [113] mime_0.12 bit64_4.0.5