r3fang / SnapATAC

Analysis Pipeline for Single Cell ATAC-seq
GNU General Public License v3.0
296 stars 124 forks source link

runMACS on filtered regions #161

Open dinhdiep opened 4 years ago

dinhdiep commented 4 years ago

Hi Rongxin,

In the code below I attempted to remove the blacklist regions and some unplaced and random contigs prior to peak calling using runMACS. However, after calling runMACS the resulting peaks file still contain peaks from the random and unplaced contigs. Does runMACS subset reads from regions available in the BMAT? Or is there a step that I am missing?

Thank you,

Dinh

`black_list <- read.table("hg38.blacklist.bed") black_list.gr = GRanges( black_list[,1], IRanges(black_list[,2], black_list[,3]) )

idy1 = queryHits(findOverlaps(x.sp@feature, black_list.gr)) idy2 = grep("chrM|random|chrUn", x.sp@feature) idy = unique(c(idy1, idy2))

plotBinCoverage( obj=x.sp, col="grey", border="grey", breaks=20, xlim=c(-6, 6) )

x.sp.rmsk = x.sp[,-idy, mat="bmat"]

Use all the cells to call peaks

bulk.peaks = runMACS(obj=x.sp.rmsk, output.prefix = "MACS_peaks", path.to.snaptools="~/.local/bin/snaptools", path.to.macs="~/.local/bin/macs2", tmp.folder="./", num.cores=16, gsize="hs", buffer.size=500, macs.options="--nomodel --shift 75 --ext 150 --qvalue 1e-1 -B --SPMR --call-summits", keep.minimal = TRUE)

bulk.peaks.gr = GRanges( bulk.peaks[,1], IRanges(bulk.peaks[,2], bulk.peaks[,3]) ) `