# The exomePeak2 Guide

#### Zhen Wei, Jia Meng

# Introduction

exomePeak2 provides technical independent peak detection and differential methylation analysis for Methylated RNA immunoprecipitation sequencing data (MeRIP-Seq). MeRIP-Seq is the primary sequencing technology for epi-transcriptomic assessment. The peak calling processes in MeRIP-Seq is sensitive to GC content biases, which are generally present in NGS-based assays. Besides, the antibody pull-down efficiency do often vary across different IP replicates, introducing another layer of unwanted technical variation. exomePeak2 addresses these challenges by introducing a robust set of computation tools tailored for MeRIP-Seq. With exomePeak2, users can perform peak calling and differential analysis through a straightforward single-step function.

# Installation

To install exomePeak2 from bioconductor, start R (version >“3.6”) and enter:

    if(!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("exomePeak2")

For order versions of R, please refer to the appropriate [Bioconductor release](https://www.bioconductor.org/about/release-announcements/}.

# Peak Calling

For the peak calling of the MeRIP-Seq experiment, exomePeak2 requires alignment results to be stored in BAM format. User need to specify the corresponding BAM directories of IP and input samples in the arguments `bam_ip` and `bam_input` of the main function `exomePeak2()`, respectively.

The following example demonstrates the peak calling from BAM and GFF files. In addition to GFF files, transcript annotation can also be provided by the bioconductor TxDb object. TxDb will be automatically downloaded if the corresponding UCSC genome name is filled in the `genome` argument.

Note that the genome sequences are necessary for the correction of GC content biases. If the `genome` argument is missing ( `= NULL` ), exomPeak2 will perform peak calling without correcting the PCR amplification bias.

    library(exomePeak2)

    set.seed(1)

    GENE_ANNO_GTF = system.file("extdata", "example.gtf", package="exomePeak2")

    f1 = system.file("extdata", "IP1.bam", package="exomePeak2")
    f2 = system.file("extdata", "IP2.bam", package="exomePeak2")
    f3 = system.file("extdata", "IP3.bam", package="exomePeak2")
    f4 = system.file("extdata", "IP4.bam", package="exomePeak2")
    IP_BAM = c(f1,f2,f3,f4)

    f1 = system.file("extdata", "Input1.bam", package="exomePeak2")
    f2 = system.file("extdata", "Input2.bam", package="exomePeak2")
    f3 = system.file("extdata", "Input3.bam", package="exomePeak2")
    INPUT_BAM = c(f1,f2,f3)

    exomePeak2(bam_ip = IP_BAM,
               bam_input = INPUT_BAM,
               gff = GENE_ANNO_GTF,
               genome = "hg19")
    ## GRangesList object of length 29:
    ## $`1`
    ## GRanges object with 1 range and 0 metadata columns:
    ##     seqnames              ranges strand
    ##        <Rle>           <IRanges>  <Rle>
    ##   1     chr7 150742912-150742936      +
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
    ##
    ## $`2`
    ## GRanges object with 1 range and 0 metadata columns:
    ##     seqnames              ranges strand
    ##        <Rle>           <IRanges>  <Rle>
    ##   2     chr7 150742987-150743011      +
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
    ##
    ## $`3`
    ## GRanges object with 1 range and 0 metadata columns:
    ##     seqnames              ranges strand
    ##        <Rle>           <IRanges>  <Rle>
    ##   3     chr7 150743062-150743136      +
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
    ##
    ## ...
    ## <26 more elements>

Besides the GRangesList object of peaks returned by the main function, exomePeak2 will export significant peaks in BED and CSV format; the files will be automatically saved in a folder named `exomePeak2_output`.

In the peak calling mode, the peak statistics are calculated from the \({\beta}_{i,1}\) terms in the following Poisson GLMs:

\[K_{i,j} \sim Poisson(\lambda_{i,j})\]

\[log2(\lambda_{i,j}) = {\beta}_{i,0} + {\beta}_{i,1}1(\rho(j)=\text{IP}) + t_{i,j}\]

\(\lambda_{i,j}\) is the expected value of read abundance for sliding window (bin) \(i\) under sample \(j\). \({\beta}_{i,0}\) is the intercept coefficient; \({\beta}_{i,1}\) is the coefficient of IP/input log2 fold change; \(1(\rho(j)=\text{IP})\) is the regression covariate, which is an indicator (binary) variable of the sample \(j\) being IP sample. \(t_{i,j}\) is the bin-specific offset accounting for the sequencing depth variation, IP efficiency variation, and the GC content bias.

\[t_{i,j} = log2(s_j) + log2(f(GC_i))\]

By default, \(t_{i,j}\) is further decomposed into the sample-specific size factor \(s_j\) and the GC content bias offset \(f(GC_i)\). \(\hat s_j\) is estimated using the median read count on the background bins, which are classified by the GMM (K = 2) from the bin’s POI values: (IP count + 1)/(input count + 1). \(\hat f(GC_i)\) is estimated using cubic splines expanded Poisson GLM on bin’s read count with bin’s fragment GC content \(GC_i\) as covariate. Knots for the cubic splines used are set at `c(0, 0.4, 0.5, 0.6, 1)`. For identifiability, the fitted values \(\hat f(GC_i)\) are centered at 1 by dividing its sample mean.

For peak calling, Poisson GLMs using the above design are fitted to identify the significantly modified sliding windows; the default detection criteria is the Wald test p-values < 1e-10 under the alternative hypothesis of \({\beta}_{i,1} > 0\).

Annotations on the columns of the output table:

- chr: the chromosome name of the peak.
- chromStart: the start of the peak on the chromosome.
- chromEnd: the end of the peak on the chromosome.
- name: the unique ID of the modification peak.
- strand: the strand of the peak on genome.
- blockCount: the block (exon) number within the peak.
- blockSizes: the widths of the blocks.
- blockStarts: the start positions of the blocks.
- geneID: the gene ID of the peak; multiple gene IDs will be returned when peak is compatible with > 1 genes.
- RPM.input: the peak RPM (Reads per million mapped reads) calculated from the pooled input samples.
- RPM.IP: the peak RPM calculated from the pooled IP samples.
- log2FC: the estimate of IP over input log2 fold change (coefficient estimates of \({\beta}_{i,1}\)).
- pvalue: the p-value of the peak, which is calculated by the Wald test over the coefficient of \({\beta}_{i,1}\).
- fdr: the adjusted Wald test p-value using Benjamini Hochberg approach.
- score: the -log10 p value of the peak.

# Differential Methylation Analysis

The following example performs differential methylation analysis (comparison of two biological conditions) on exon regions defined by transcript annotation.

    f1 = system.file("extdata", "treated_IP1.bam", package="exomePeak2")
    TREATED_IP_BAM = c(f1)
    f1 = system.file("extdata", "treated_Input1.bam", package="exomePeak2")
    TREATED_INPUT_BAM = c(f1)

    exomePeak2(bam_ip = IP_BAM,
               bam_input = INPUT_BAM,
               bam_input_treated = TREATED_INPUT_BAM,
               bam_ip_treated = TREATED_IP_BAM,
               gff = GENE_ANNO_GTF,
               genome = "hg19")
    ## GRangesList object of length 2:
    ## $`1`
    ## GRanges object with 1 range and 0 metadata columns:
    ##     seqnames              ranges strand
    ##        <Rle>           <IRanges>  <Rle>
    ##   1     chr7 150910255-150910629      -
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
    ##
    ## $`2`
    ## GRanges object with 1 range and 0 metadata columns:
    ##     seqnames              ranges strand
    ##        <Rle>           <IRanges>  <Rle>
    ##   2     chr7 150934486-150935360      +
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

In the differential methylation mode, exomePeak2 will first perform peak calling on the exonic regions by pulling the control and treated samples. Next, it will conduct the differential methylation analysis over the read counts of significant peaks using an interactive Poisson GLM design.

The peak statistics calculated in the differential methylation setting are based on the interactive coefficient \({\beta}_{i,3}\) of the following regression equation:

\[log2(\lambda_{i,j}) = {\beta}_{i,0} + {\beta}_{i,1}1(\rho(j)=\text{IP}) + {\beta}_{i,2}1(\rho(j)=\text{Treated}) + {\beta}_{i,3}1(\rho(j)=\text{IP&Treated}) + t_{i,j}\]

In the output of differential methylation analysis, The diff.log2FC and pvalue are calculated from the coefficient estimates of the interactive term \({\beta}_{i,3}\). The differential log2FC here is interpreted as log( IP over input FC in treated group / IP over input FC in control group). The offset terms \(t_{i,j}\) used in the interactive GLM are calculated using the same method as the peak calling.

# Contact

If you encounter any problems during use, please contact the maintainer of exomePeak2:

Zhen Wei : zhen.wei01@xjtlu.edu.cn

# References

1. Hansen, K. D., Irizarry, R. A., & Wu, Z. (2012). Removing technical variability in RNA-seq data using conditional quantile normalization. Biostatistics, 13(2), 204-216.

2. Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology, 15(12), 1-21.

3. Benjamini, Y., & Speed, T. P. (2012). Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic acids research, 40(10), e72-e72.

4. Meng, J., Lu, Z., Liu, H., Zhang, L., Zhang, S., Chen, Y., … & Huang, Y. (2014). A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/Bioconductor package. Methods, 69(3), 274-281.

# Session Info

    sessionInfo()
    ## R version 4.2.1 (2022-06-23 ucrt)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.66.0
    ##  [3] rtracklayer_1.58.0                Biostrings_2.66.0
    ##  [5] XVector_0.38.0                    exomePeak2_1.10.0
    ##  [7] SummarizedExperiment_1.28.0       Biobase_2.58.0
    ##  [9] GenomicRanges_1.50.0              GenomeInfoDb_1.34.0
    ## [11] IRanges_2.32.0                    S4Vectors_0.36.0
    ## [13] BiocGenerics_0.44.0               MatrixGenerics_1.10.0
    ## [15] matrixStats_0.62.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] bitops_1.0-7             bit64_4.0.5              filelock_1.0.2
    ##  [4] RColorBrewer_1.1-3       progress_1.2.2           httr_1.4.4
    ##  [7] tools_4.2.1              bslib_0.4.0              utf8_1.2.2
    ## [10] R6_2.5.1                 DBI_1.1.3                colorspace_2.0-3
    ## [13] tidyselect_1.2.0         prettyunits_1.1.1        DESeq2_1.38.0
    ## [16] curl_4.3.3               bit_4.0.4                compiler_4.2.1
    ## [19] textshaping_0.3.6        cli_3.4.1                xml2_1.3.3
    ## [22] DelayedArray_0.24.0      labeling_0.4.2           sass_0.4.2
    ## [25] scales_1.2.1             genefilter_1.80.0        speedglm_0.3-4
    ## [28] rappdirs_0.3.3           systemfonts_1.0.4        stringr_1.4.1
    ## [31] digest_0.6.30            Rsamtools_2.14.0         rmarkdown_2.17
    ## [34] pkgconfig_2.0.3          htmltools_0.5.3          dbplyr_2.2.1
    ## [37] fastmap_1.1.0            rlang_1.0.6              RSQLite_2.2.18
    ## [40] farver_2.1.1             jquerylib_0.1.4          BiocIO_1.8.0
    ## [43] generics_0.1.3           jsonlite_1.8.3           mclust_6.0.0
    ## [46] BiocParallel_1.32.0      dplyr_1.0.10             RCurl_1.98-1.9
    ## [49] magrittr_2.0.3           GenomeInfoDbData_1.2.9   Matrix_1.5-1
    ## [52] Rcpp_1.0.9               munsell_0.5.0            fansi_1.0.3
    ## [55] lifecycle_1.0.3          stringi_1.7.8            yaml_2.3.6
    ## [58] MASS_7.3-58.1            zlibbioc_1.44.0          BiocFileCache_2.6.0
    ## [61] grid_4.2.1               blob_1.2.3               parallel_4.2.1
    ## [64] crayon_1.5.2             lattice_0.20-45          splines_4.2.1
    ## [67] GenomicFeatures_1.50.0   annotate_1.76.0          hms_1.1.2
    ## [70] KEGGREST_1.38.0          locfit_1.5-9.6           knitr_1.40
    ## [73] pillar_1.8.1             rjson_0.2.21             geneplotter_1.76.0
    ## [76] codetools_0.2-18         biomaRt_2.54.0           XML_3.99-0.12
    ## [79] glue_1.6.2               evaluate_0.17            png_0.1-7
    ## [82] vctrs_0.5.0              gtable_0.3.1             assertthat_0.2.1
    ## [85] cachem_1.0.6             ggplot2_3.3.6            xfun_0.34
    ## [88] xtable_1.8-4             restfulr_0.0.15          ragg_1.2.4
    ## [91] survival_3.4-0           tibble_3.1.8             GenomicAlignments_1.34.0
    ## [94] AnnotationDbi_1.60.0     memoise_2.0.1            ellipsis_0.3.2