exomePeak2 {exomePeak2} | R Documentation |
exomePeak2
conducts peak calling and peak statistics calculation from BAM files of a MeRIP-seq experiment.
The function integrates the following steps of a standard MeRIP-seq data analysis pipeline.
Check and index the BAM files with scanMeripBAM
.
Call modification peaks on exons with exomePeakCalling
.
Calculate offset factors of GC content biases with normalizeGC
.
Calculate (differential) modification statistics with the generalized linear model (GLM) using glmM
or glmDM
.
Export the peaks/sites statistics with user defined format by exportResults
.
See the help pages of the corresponding functions for the complete documentation.
exomePeak2( bam_ip = NULL, bam_input = NULL, bam_treated_ip = NULL, bam_treated_input = NULL, txdb = NULL, bsgenome = NULL, genome = NA, gff_dir = NULL, mod_annot = NULL, paired_end = FALSE, library_type = c("unstranded", "1st_strand", "2nd_strand"), fragment_length = 100, binding_length = 25, step_length = binding_length, min_peak_width = fragment_length/2, max_peak_width = Inf, pc_count_cutoff = 0, bg_count_cutoff = 5, p_cutoff = 1e-05, p_adj_cutoff = NULL, log2FC_cutoff = 0, parallel = 1, background_method = c("Gaussian_mixture", "m6Aseq_prior", "manual", "all"), manual_background = NULL, correct_GC_bg = FALSE, qtnorm = FALSE, glm_type = c("DESeq2", "Poisson", "NB"), LFC_shrinkage = c("apeglm", "ashr", "Gaussian", "none"), export_results = TRUE, export_format = c("CSV", "BED", "RDS"), table_style = c("bed", "granges"), save_plot_GC = TRUE, save_plot_analysis = FALSE, save_plot_name = "", save_dir = "exomePeak2_output", peak_calling_mode = c("exon", "full_tx", "whole_genome") )
bam_ip |
a |
bam_input |
a |
bam_treated_ip |
a |
bam_treated_input |
a If the bam files do not contain treatment group, user should only fill the arguments of |
txdb |
a If the |
bsgenome |
a If the |
genome |
a By default, the argument = NA, it should be provided when the |
gff_dir |
optional, a |
mod_annot |
a If user provides the single based RNA modification annotation, exomePeak2 will perform reads count and (differential) modification quantification on the provided annotation. The single base annotation will be flanked by length = floor(fragment_length - binding_length/2) to account for the fragment length of the sequencing library. |
paired_end |
a |
library_type |
a
|
fragment_length |
a positive integer number for the expected fragment length in nucleotides; default |
binding_length |
a positive integer number for the expected binding length of the anti-modification antibody in IP samples; default |
step_length |
a positive integer number for the shift distances of the sliding window; default |
min_peak_width |
a |
max_peak_width |
a |
pc_count_cutoff |
a |
bg_count_cutoff |
a |
p_cutoff |
a |
p_adj_cutoff |
a |
log2FC_cutoff |
a |
parallel |
a |
background_method |
a In order to accurately account for the technical variations, it is often important to estimate the sequencing depth and GC content linear effects on windows without modification signals. The following methods are supported in
|
manual_background |
a |
correct_GC_bg |
a If |
qtnorm |
a If |
glm_type |
a
By default, the DESeq2 GLMs are fitted on the data set with > 1 biological replicates for both the IP and input samples, the Poisson GLM will be fitted otherwise. |
LFC_shrinkage |
a see |
export_results |
a |
export_format |
a |
table_style |
a |
save_plot_GC |
a |
save_plot_analysis |
a |
save_plot_name |
a |
save_dir |
a |
peak_calling_mode |
a
P.S. The full transcript mode and the whole genome mode demand big memory size (> 4GB) for large genomes. |
exomePeak2
call RNA modification peaks and calculate peak statistics from BAM files of a MeRIP-seq experiment.
The transcript annotation (from either the TxDb
object or the GFF file) should be provided to perform analysis on exons.
The BSgenome
object is also required to perform the GC content bias adjustment.
If the bsgenome
and the genome
arguments are not provided (= NULL
), the downstream analysis will proceed without GC content bias corrections.
If the BAM files in treated samples are provided at the arguments bam_treated_ip
and bam_treated_input
, the statistics of differential modification analysis on peaks/sites will be reported.
Under the default setting, exomePeak2
will save the results of (differential) modification analysis under a folder named 'exomePeak2_output'
.
The results generated include a BED file and a CSV table that stores the locations and statistics of the (differential) modified peaks/sites.
When performing differential analysis using function exomePeak2(), the sequencing depth for the interactive GLM will be estimated from the background features, which by default are the disjoint regions of the detected peaks on exons.
a SummarizedExomePeak
object.
exomePeakCalling
, glmM
, glmDM
, normalizeGC
, exportResults
, plotLfcGC
#Specify File Directories 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) # Peak Calling sep <- exomePeak2(bam_ip = IP_BAM, bam_input = INPUT_BAM, gff_dir = GENE_ANNO_GTF, genome = "hg19", paired_end = FALSE) sep # Differential Modification Analysis on Modification Peaks (Comparison of Two Conditions) 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) sep <- 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) sep # Modification Level Quantification on Single Based Modification Annotation f2 = system.file("extdata", "mod_annot.rds", package="exomePeak2") MOD_ANNO_GRANGE <- readRDS(f2) sep <- exomePeak2(bam_ip = IP_BAM, bam_input = INPUT_BAM, gff_dir = GENE_ANNO_GTF, genome = "hg19", paired_end = FALSE, mod_annot = MOD_ANNO_GRANGE) sep # Differential Modification Analysis on Single Based Modification Annotation sep <- 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, mod_annot = MOD_ANNO_GRANGE) sep