alternative to HOMER's mergePeaks

Question

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 multiIntersectBedthat 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?

Answer 1

Via BEDOPS, you could combine bedops --merge with process substitutions sorting N files:

$ bedops --merge <(sort-bed 1.bed) <(sort-bed 2.bed) ... <(sort-bed N.bed) > answer.bed

The peak memory usage of a typical bedops --merge run should be around 10MB, if I remember correctly, as a result of using sorted inputs.

The memory usage of each sort-bed process can take up the size of each file. You can specify --max-mem and --tmpdir options to instead use disk space to sort files. You may want to do this if the sum of file sizes is equal to or greater than ~512GB, or however much system memory is available (typically less).

If your inputs are already sorted, or if your ATACseq pipeline is set up to create sorted files, then sort-bed does not need to be run again, and you can just use the merge operation on its own:

$ bedops --merge 1.bed 2.bed ... N.bed > answer.bed

The memory usage is then only what is needed for merging, as discussed above.

Using bedops --merge will merge adjacent ("bookended") intervals, so long as they are in sort-order, e.g.:

$ echo -e "chr1\t4658186\t4658322\nchr1\t4658322\t4658775" | bedops --merge -
chr1    4658186 4658775

Notes:

1. The file answer.bed will be sorted per sort-bed for downstream operations, such as bedmap association of merged regions to peaks. Additional sorting is not needed.
2. Using sort-bed will often be faster than Unix sort on BED files.
3. The limit of N on most platforms will be 1021 files. Most OSes have a limit of 1024 file descriptors, and three are reserved for standard IO. This soft limit can be dialed upwards on some hosts where the user has admin access.

References:

1. https://bedops.readthedocs.io/en/latest/content/summary.html#bedops
2. https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#merge-m-merge
3. https://bedops.readthedocs.io/en/latest/content/reference/file-management/sorting/sort-bed.html