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 MACS2
. The 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
.
In image format, I want the leftmost example, where the three peaks come from three separate files instead of all being described in the same file:
Additionally, bedtools does have a poorly documented function called multiIntersectBed
that appears to run the cat
and 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 1.bed
and 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?