Recently, our lab ran a 192 sample experiment through our ATACseq pipeline. In doing so, HOMER's
mergePeaks told us that our 512 GB RAM server had too little memory to process all of the samples at once. So here we are trying to replace it with an alternative.
It seems that
bedTools merge would be the way to go to replace it. However, going into this step, each sample has its own bed file of peaks called by
mergeBed only takes in 1 input file that must be sorted by chromosome and peak start position before we run it. I can run a
cat call on all of the samples followed by a
sort, but those individual sample bed files are already sorted. Is there a way to run
mergeBed with multiple input files at once? I tried
mergeBed -i 1.bed -i 2.bed > out.bed, but that is equivalent to
mergeBed -i 1.bed > out.bed.
Additionally, bedtools does have a poorly documented function called
multiIntersectBedthat appears to run the
sort together, but is both slower and doesn't merge bookended sections--eg files look like this:
chr1 4658186 4658322 chr1 4658322 4658775 chr1 4658775 4658777 chr1 4671049 4671351
When they should look like this:
chr1 4658186 4658777 # This was 3 lines before chr1 4671049 4671351
In summary, I have sorted bed files
2.bed from the same experiment, and I want to merge all bookended and overlapping peak sections between them. How do I do so in 1 step that doesn't involve HOMER's
mergePeaks? Thank you!
EDIT: Here's the speed of my current solutions at n=192:
HOMER: Unstable `multiIntersectBed | mergeBed` real 3m53.230s user 6m28.166s sys 0m14.912s `cat | sort | mergeBed` real 0m24.860s user 0m30.939s sys 0m2.946s `cat > tmp1; sort > tmp2; mergeBed > out; rm tmp?` real 0m11.824s user 0m20.710s sys 0m2.021s `sort -m | mergeBed` real 0m17.905s user 0m25.112s sys 0m1.816s `sort -m --batch-size=192 | mergeBed` real 0m14.687s user 0m25.589s sys 0m0.747s `sort -m > tmp ; mergeBed -i tmp ; rm tmp` real 0m17.472s user 0m15.251s sys 0m2.212s
So... Followup questions, why does
sort -m take longer than
sort | cat especially when the sys time is indeed smaller? And why does piping take more time than creating temporary files?