Identify peaks in ChIP-seq data with biological replicates using a one-sided quasi-likelihood ratio test in quasi-Poisson or quasi-negative binomial models.
This package provides code to call peaks in ChIP-seq data with biological replicates using the BinQuasi algorithm of Goren, Liu, Wang, and Wang (2018) doi.org/10.1093/bioinformatics/bty227.
The BinQuasi package for R can be installed from Github using devtools following the code below.
devtools::install_github("emilygoren/BinQuasi", args = "--preclean", build_vignettes = TRUE)library(BinQuasi)
BinQuasi accepts sorted and indexed BAM files (note that it does not perform genome alignment of raw reads). If your BAM files are not indexed and sorted, we recommend using samtools.
Once installed, BinQuasi calls peaks with the function "BQ()." Below is code to run BinQuasi with all default settings, where the sorted and indexed BAM files are stored in the directory specified by "fpath" under the file names "C1.bam", " C2.bam" and "I1.bam", "I2.bam" for ChIP and input files, respectively.
fpath <- paste0(system.file(package = 'BinQuasi'), '/extdata/')results <- BQ(fpath,ChIP.files = c('C1.bam', 'C2.bam'),control.files = c('I1.bam', 'I2.bam'))head(results$peaks)
See the package documentation for information on changing the default settings.
The code below saves the called peaks in BED format in the file "BinQuasiPeaks.bed".
# Sort peaks by p-valueopeaks <- results$peaks[order(results$peaks$P.val),]# Name the peaks by rankopeaks$name <- paste0('BQ_Peak_', 1:nrow(opeaks))# Save as .bed file, setting the scores to be -log10(p-value)bedout <- data.frame(chrom = opeaks$chr,chromStart = opeaks$start,chromEnd = opeaks$end,name = opeaks$name,score = -log10(opeaks$P.val),strand = c(rep(".", nrow(opeaks))))head(bedout)write.table(bedout, file="BinQuasiPeaks.bed", quote = FALSE, sep = "\t", row.names = FALSE, col.names = FALSE)