exomePeak2: README.md

2019-12-20

exomePeak2 provides bias awared quantification and peak detection on Methylated RNA immunoprecipitation sequencing data (MeRIP-Seq). MeRIP-Seq is a commonly applied sequencing technology to measure the transcriptome-wide location and abundance of RNA modification sites under a given cellular condition. However, the quantification and peak calling in MeRIP-Seq are sensitive to PCR amplification bias which is prevalent in next generation sequencing (NGS) techniques. In addition, the RNA-Seq based count data exhibits biological variation in small reads count. exomePeak2 collectively address these challanges by introducing a rich set of robust data science models tailored for MeRIP-Seq. With exomePeak2, users can perform peak calling, modification site quantification, and differential analysis with a straightforward one-step function. Alternatively, users could customize methods for their own analysis through multi-step functions and diagnostic plots.

To install exomePeak2 from Github, use the following codes.

For peak calling of MeRIP-Seq experiment, exomePeak2 demands the reads alignment results in BAM files. Users can specify the biological replicates of the IP and input samples by a character vector of the corresponding BAM directories at the arguments bam_ip and bam_input separately.

In the following example, the transcript annotation is provided using GFF files. Transcript annotation can also be provided by the TxDb object. exomePeak2 will automatically download the TxDb if the genome argument is filled with the corresponding UCSC genome name.

The genome sequence is required to conduct GC content bias correction. If the genome argument is missing ( = NULL ), exomPeak2 will perform peak calling without correcting the GC content bias.

exomePeak2 will export the modification peaks in formats of BED file and CSV table, the data will be saved automatically under a folder named by exomePeak2_output.

The modification peak statistics are derived from the βi, 1 terms in the following linear regression design.

log2(Qi, j) = βi, 0 + βi, 1I(ρ(j) = IP)+ti, j

Qi, j is the expected value of reads abundence of modification i under sample j. βi, 0 is the intercept coefficient, βi, 1 is the coefficient for IP/input log2 fold change, I(ρ(j)=IP) is the regression covariate that is the indicator variable for the sample j being IP sample. ti, j is the regression offset that account for the sequencing depth variation and the GC content biases.

Under the default settings, the linear models fitted are the regularized GLM (Generalized Linear Model) of NB developed by DESeq2. If one of the IP and input group has no biological replicates, Poisson GLMs will be fitted to the modification peaks.

Explaination over the columns of the exported table:

The code below could conduct differential modification analysis (Comparison of Two Conditions) on exon regions defined by the transcript annotation.

In differential modification mode, exomePeak2 will first perform Peak calling on exon regions using both the control and treated samples. Then, it will conduct the differential modification analysis on peaks reported from peak calling using an interactive GLM.

In differential modification mode, exomePeak2 will export the differential modification peaks in formats of BED file and CSV table, the data will also be saved automatically under a folder named by exomePeak2_output.

The peak statistics in differential modification setting are derived from the interactive coefficient βi, 3 in the following regression design of the NB GLM:

log2(Qi, j) = βi, 0 + βi, 1I(ρ(j)=IP)+βi, 2I(ρ(j) = Treatment)+βi, 3I(ρ(j) = IP&Treatment)+ti, j

Explaination for the additional table columns:

exomePeak2 supports the modification quantification and differential modification analysis on single based modification annotation. The modification sites with single based resolution can provide a more accurate mapping of modification locations compared with the peaks called directly from the MeRIP-seq datasets.

Some of the datasets in epitranscriptomics have single based resolution, e.x. Data generated by the m6A-CLIP-Seq or m6A-miCLIP-Seq techniques. Reads count on the single based modification sites could provide a more accurate and consistent quantification on MeRIP-Seq experiments with single based annotation.

exomePeak2 will automatically initiate the mode of single based modification quantification by providing a sigle based annotation file under the argument mod_annot.

The single based annotation information should be provided to the exomePeak2 function in the format of a GRanges object.