---
title: "The *exomePeak2* User's Guide"
author: 
- name: Zhen Wei
  affiliation:
  - Department of Biological Sciences, Xi’an Jiaotong-Liverpool University, Suzhou, Jiangsu, 215123, China
  - Institute of Integrative Biology, University of Liverpool, L7 8TX, Liverpool, United Kingdom
  email: Zhen.Wei01@xjtlu.edu.cn
- name: Jia Meng
  affiliation:
  - Department of Biological Sciences, Xi’an Jiaotong-Liverpool University, Suzhou, Jiangsu, 215123, China
  - Institute of Integrative Biology, University of Liverpool, L7 8TX, Liverpool, United Kingdom
  email: Jia.Meng@xjtlu.edu.cn
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
graphics: yes
vignette: >
  %\VignetteIndexEntry{The exomePeak2 user's guide}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignettePackage{exomePeak2}
  %\VignetteEncoding{UTF-8}
---

```{r para, echo = FALSE, results='hide'}
knitr::opts_chunk$set(dev="png",fig.show="hold",
               fig.width=8,fig.height=4.5,fig.align="center",
               message=FALSE,collapse=TRUE)
```

# Introduction

**exomePeak2** provides bias-aware quantification and peak detection for Methylated RNA immunoprecipitation sequencing data (**MeRIP-Seq**). *MeRIP-Seq* is a commonly applied sequencing technology that can measure the location and abundance of RNA modification sites under given cell line conditions. However, quantification and peak calling in *MeRIP-Seq* are sensitive to PCR amplification biases, which generally present in next-generation sequencing (NGS) technologies. In addition, the count data generated by **RNA-Seq** exhibits significant biological variations between biological replicates. *exomePeak2* collectively address the challenges by introducing a series of robust data science tools tailored for MeRIP-Seq. Using *exomePeak2*, users can perform peak calling, modification site quantification and differential analysis through a straightforward single-step function. Alternatively, multi-step functions can be used to generate diagnostic plots and perform customized analyses.

# Installation

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

```{r, eval=FALSE}
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 that the alignment result be stored in **BAM** format. User can specify the corresponding **BAM** directories of IP and input samples at the arguments `bam_ip` and `bam_input` respectively. 

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

The genome sequence is necessary for the correction of GC content bias. If both the `genome` and the `bsgenome` arguments are missing ( `= NULL` ), exomPeak2 will perform peak calling without correcting PCR amplification bias.

```{r, eval = TRUE}
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_dir = GENE_ANNO_GTF,
           genome = "hg19",
           paired_end = FALSE)
```

exomePeak2 will export the significant peaks in the format of **BED** file and **CSV** table, and the output results will be automatically saved in a folder named `exomePeak2_output`.

If only IP and input samples are provided, the peak statistics are calculated from the ${\beta}_{i,1}$ terms in the following linear regression design:

$$log2(Q_{i,j}) = {\beta}_{i,0} + {\beta}_{i,1}1(\rho(j)=\text{IP}) + t_{i,j}$$

$Q_{i,j}$ is the expected value of read abundance of peak $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 the indicator variable of the sample $j$ being IP sample. $t_{i,j}$ is the peak specific offsets account for the sequencing depth variation and the GC content biases. $t{i,j}$ is estimated by the median regression method implemented in the package **cqn**[1]; the median regression is smoothed using the natural cubic splines with the knots defined in function `cqn::cqn()`, and the threshold for read count for regression analysis is set at 50. By default, the quantile normalization in **cqn** is turned off, so the peak specific offsets are calculated only based on the sequencing depth and GC content biases.

In the default setting, the fitted linear models are the regularized **GLM (Generalized Linear Model)** of NB implemented in **DESeq2**[2]. If there is no biological variation in one of the IP and input groups, Poisson GLM will be fitted for peak detection. Both the size factors used in **cqn** and **DESeq2** are estimated by the robust size factor estimator defined in the function `DESeq2::estimateSizeFactors()` with type = "ratio".

For peak calling, 2 rounds of GLM will be fitted. The first round of GLM fitting will identify the significantly modified sliding windows; the detection criteria is Wald test p-values < 1e-05 under the alternative hypothesis of ${\beta}_{i,1} > 0$. The second round of GLM fitting will re-calculate the peak statistics from the merged sliding windows (peaks). For the two rounds of GLM fitting, the correction offsets $t_{i,j}$ will be recalculated from the GC contents of the sliding windows and the peaks, respectively.

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.
- ***score***: the -log2 p value of the peak.
- ***strand***: the strand of the peak on genome.
- ***thickStart***: the start position of the peak.
- ***thickEnd***: the end position of the peak.
- ***itemRgb***: the column for the RGB encoded color in BED file visualization.
- ***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.
- ***ReadsCount.input***: the total read count of input samples.
- ***ReadsCount.IP***: the total read count of IP samples.
- ***log2FoldChange***: the estimate of IP over input log2 fold change (coefficient estimates of ${\beta}_{i,1}$ in GLM), the Bayesian estimation implemented in the bioconductor package `apeglm`[3] will be returned if the regularized NB GLMs are fitted using DESeq2.
- ***pvalue***: the Wald test p value on the ${\beta}_{i,1}$ coefficient in the GLM.
- ***padj***: the adjusted Wald test p value using Benjamini Hochberg approach.

# Differential Modification Analysis 

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

In the differential modification mode, exomePeak2 will first perform Peak calling on the exon regions using both the control and treated samples. Next, it will conduct the differential modification analysis on the significant peaks using GLMs of interactive designs.

```{r, eval = TRUE}
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_treated_input = TREATED_INPUT_BAM,
           bam_treated_ip = TREATED_IP_BAM,
           gff_dir = GENE_ANNO_GTF,
           genome = "hg19",
           paired_end = FALSE)
```

In the differential modification mode, exomePeak2 will export the differentially modified peaks in the format of **BED** file and **CSV** table, and the outputs will also be saved in a folder named `exomePeak2_output`.

The peak statistics calculated in the differential modification setting are based on the interactive coefficient ${\beta}_{i,3}$ in the following regression design:

$$log2(Q_{i,j}) = {\beta}_{i,0} + {\beta}_{i,1}1(\rho(j)=\text{IP}) + {\beta}_{i,2}1(\rho(j)=\text{Treatment}) + {\beta}_{i,3}1(\rho(j)=\text{IP&Treatment}) + t_{i,j}$$

Annotations on the additional columns of the output table for differential analysis:

- ***ModLog2FC_control***: the IP/input log2 fold enrichment in the control condition.
- ***ModLog2FC_treated***: the IP/input log2 fold enrichment in the treatment condition.
- ***DiffModLog2FC***: the differential log2 fold change estimates (coefficient estimates of the interactive term ${\beta}_{i,3}$); similarly, the Bayesian estimation by `apeglm` will be applied given regularized NB GLMs of DESeq2.
- ***pvalue***: the Wald test p value on the interactive coefficient of the GLM.
- ***padj***: the adjusted Wald test p value using the BH approach.

By the default setup of *exomePeak2()*, the sequencing depth for the interactive GLM will be estimated from the background features, which are the disjoint regions of the detected peaks on exons.

# Reference based Quantification and Differencial Analysis

exomePeak2 supports modification quantification and differential modification analysis based on single nucleotides modification annotation. Some data sets in epitranscriptomics have single-base resolution, e.x. Data generated by the *m6A-CLIP-Seq* or *m6A-miCLIP-Seq* technologies. Compared with the peaks called from MeRIP-Seq, sites identified by single-base resolution techniques often have higher precision. 

By eliminating of the technical variation introduced by the peak lengths, the read counting at the single-base modification sites can provide more accurate and reliable quantification for the *MeRIP-Seq* experiments. exomePeak2 will automatically initiate the mode of single-base modification quantification/differential analysis when user provide a reference object in the argument `mod_annot`. 

The GLM designs used in the reference based mode are identical with the quantification and differential analysis under the peak calling mode. The difference is that only one round of GLM fitting is performed in the reference mode by treating the the single-base modification site as the detected peaks.

The single-base annotation information should be provided to the exomePeak2 function in *GRanges* object, the GRanges object can be imported directly from the **BED** files of single-base annotation files using `dplyr::import()`.

```{r, eval = TRUE}
f2 = system.file("extdata", "mod_annot.rds", package="exomePeak2")

MOD_ANNO_GRANGE <- readRDS(f2)

exomePeak2(bam_ip = IP_BAM,
           bam_input = INPUT_BAM,
           gff_dir = GENE_ANNO_GTF,
           genome = "hg19",
           paired_end = FALSE,
           mod_annot = MOD_ANNO_GRANGE)
```

Under the reference based mode, exomePeak2 will export the analysis results in the format of **BED** and **CSV** files, while the order of the rows in the files will match the corresponding ranges of the single-base `GRanges`.

# Multi-step Functions

The exomePeak2 package can achieve peak calling and peak statistics calculation with multiple functions.

**1. Check the BAM files of MeRIP-seq experiments**

The *scanMeripBAM()* is used to organize the *BAM* files used in a MeRIP-Seq analysis, it also specify the important BAM reading parameters such as the paired end library type and the FLAG filters.
```{r, eval = TRUE}
MeRIP_Seq_Alignment <- scanMeripBAM(
                         bam_ip = IP_BAM,
                         bam_input = INPUT_BAM,
                         paired_end = FALSE
                        )
```

For MeRIP-seq experiment with interactive design (contain control and treatment groups), we can add the treated input and IP files with the arguments `bam_treated_input` and `bam_treated_ip`:
```{r, eval = TRUE}
MeRIP_Seq_Alignment <- scanMeripBAM(
    bam_ip = IP_BAM,
    bam_input = INPUT_BAM,
    bam_treated_input = TREATED_INPUT_BAM,
    bam_treated_ip = TREATED_IP_BAM,
    paired_end = FALSE
  ) 
```


**2. Perform peak calling on exons**

After checking the MeRIP-seq BAMs, peak calling can be performed with *exomePeakCalling()*. The function will first generate sliding windows on exons of the annotation file, and it will then count the read ans it GLMs to detect the significantly modified peaks. If either `genome` or `bsgenome` arguments are provided, the GC content bias correction will be performed while peak calling.

```{r, eval = TRUE}
SummarizedExomePeaks <- exomePeakCalling(merip_bams = MeRIP_Seq_Alignment,
                                         gff_dir = GENE_ANNO_GTF,
                                         genome = "hg19") 
```

Alternatively, user can provide single-base RNA modification annotation at `mod_annot`, so that the quantification will be performed on the reference sites.

```{r, eval = TRUE}
SummarizedExomePeaks <- exomePeakCalling(merip_bams = MeRIP_Seq_Alignment,
                                         gff_dir = GENE_ANNO_GTF,
                                         genome = "hg19",
                                         mod_annot = MOD_ANNO_GRANGE) 
```

**3. Estimate correction factors**

Before performing modification level quantification or differential analysis on the peaks or sites, the normalization methods used will be controlled by the functions **estimateSeqDepth()** and **normalizeGC()**. User can combine multiple normalization approaches, such as adding the conditional quantile normalization introduced in [2].
```{r, eval = TRUE}
SummarizedExomePeaks <- estimateSeqDepth(SummarizedExomePeaks)
SummarizedExomePeaks <- normalizeGC(SummarizedExomePeaks)
```

P.S. The GC content bias factors can still be estimated by *normalizeGC()* if no genome sequence information are provided in the previous steps, but genome information need to be provided at *normalizeGC()* via the `bsgenome` argument.

In addition, the sequencing depth and GC content offset can be calculated according to 3 types of feature ranges: "All" (default), "Background" and "Modification". The impact of different estimation scopes are highly data dependent, and reasonable changes of these may significantly improve the biological outcomes of the downstream analysis. For example, in differential analysis, it is recommended to use `estimateSeqDepth(SummarizedExomePeaks, from='Background')` so that the direction of differential modification directions can be more consistent with the expectation of the perturbed protein factors. 

**4. Report the GLM statistics**

Using the normalization factors calculated above, the GLM statistics of the IP/input design can be calculated with the function *glmM()*:
```{r, eval = FALSE}
SummarizedExomePeaks <- glmM(SummarizedExomePeaks) 
```

If the treated IP and input BAM are provided, *glmDM()* function can be used to conduct differential modification analysis on modification Peaks with interactive GLM:
```{r, eval = TRUE, warning=FALSE}
SummarizedExomePeaks <- glmDM(SummarizedExomePeaks)
```

All of the GLM statistics and normalization factors calculated will be contained in the `SummarizedExomePeaks` object, which can be saved (ex. with `saveRDS()`) to store the intermediate results in a given pipeline of analysis. 

**5. Scatter plot between GC content and log2 fold changes (LFCs).**

A scatter plot can be produced with *plotLfcGC()* to visualize the trend of LFCs with GC. It is recommended to plot it twice before and after conducting the GC content normalization using *normalizeGC()*.

```{r, eval = TRUE, fig.align='center', fig.height = 2.8, fig.width = 5}
plotLfcGC(SummarizedExomePeaks, point_size = 1, xlim = c(0.4, 0.85)) 
```

**6. Bar plot of sequencing depth estimates**

The library size factors can be visualized and compared using *plotSizeFactors()*. The label will allow comparison of 3 types of feature ranges for size factor calculations, including the background, modification and both regions (All).

```{r, eval = TRUE}
plotSizeFactors(SummarizedExomePeaks)
```

**7. Export the modification peaks and their statistics**

The analysis results can be saved on the system in formats of **CSV**, **BED**, or **RDS**.
```{r, eval = TRUE, fig.width=10, fig.height=7}  
exportResults(SummarizedExomePeaks, format = "BED") 
```

Please note that *exportResults()* can also specify the filtering thresholds and style of the result table.

# Contact

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

**Zhen Wei** :  <zhen.wei01@xjtlu.edu.cn>

# References
1. KD Hansen, RA Irizarry, and Z Wu, Removing technical variability in RNA-seq data using conditional quantile normalization. Biostatistics 2012 vol. 13(2) pp. 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:550. https://doi.org/10.1186/s13059-014-0550-8

3. Zhu A, Ibrahim JG, Love MI (2018). “Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences.” Bioinformatics. doi: 10.1093/bioinformatics/bty895.

# Session Info
```{r}
sessionInfo()
```


