Answers.Body: <p>Via BEDOPS, you could combine <code>bedops --merge</code> with process substitutions sorting <code>N</code> files:</p>
  <pre><code>$ bedops --merge &lt;(sort-bed 1.bed) &lt;(sort-bed 2.bed) ... &lt;(sort-bed N.bed) &gt; answer.bed
  </code></pre>
  <p>The peak memory usage of a typical <code>bedops --merge</code> run should be around 10MB, if I remember correctly, as a result of using sorted inputs.</p>
  <p>The memory usage of each <code>sort-bed</code> process can take up the size of each file. You can specify <code>--max-mem</code> and <code>--tmpdir</code> 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).</p>
  <p>If your inputs are already sorted, or if your ATACseq pipeline is set up to create sorted files, then <code>sort-bed</code> does not need to be run again, and you can just use the merge operation on its own:</p>
  <pre><code>$ bedops --merge 1.bed 2.bed ... N.bed &gt; answer.bed
  </code></pre>
  <p>The memory usage is then only what is needed for merging, as discussed above.</p>
  <p>Using <code>bedops --merge</code> will merge adjacent (&quot;bookended&quot;) intervals, so long as they are in sort-order, <em>e.g.</em>:</p>
  <pre><code>$ echo -e &quot;chr1\t4658186\t4658322\nchr1\t4658322\t4658775&quot; | bedops --merge -
  chr1    4658186 4658775
  </code></pre>
  <p>Notes:</p>
  <ol>
  <li>The file <code>answer.bed</code> will be sorted per <code>sort-bed</code> for downstream operations, such as <code>bedmap</code> association of merged regions to peaks. Additional sorting is not needed.</li>
  <li>Using <code>sort-bed</code> will often be faster than Unix <code>sort</code> on BED files.</li>
  <li>The limit of <code>N</code> 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 <a href="https://serverfault.com/questions/48717/practical-maximum-open-file-descriptors-ulimit-n-for-a-high-volume-system">dialed upwards</a> on some hosts where the user has admin access.</li>
  </ol>
  <p>References:</p>
  <ol>
  <li><a href="https://bedops.readthedocs.io/en/latest/content/summary.html#bedops" rel="nofollow noreferrer">https://bedops.readthedocs.io/en/latest/content/summary.html#bedops</a></li>
  <li><a href="https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#merge-m-merge" rel="nofollow noreferrer">https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#merge-m-merge</a></li>
  <li><a href="https://bedops.readthedocs.io/en/latest/content/reference/file-management/sorting/sort-bed.html" rel="nofollow noreferrer">https://bedops.readthedocs.io/en/latest/content/reference/file-management/sorting/sort-bed.html</a></li>
  </ol>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-21T23:27:20.747Z
Answers.LastActivityDate: 2021-06-22T01:13:10.837Z
Answers.LastEditDate: 2021-06-22T01:13:10.837Z
Answers.UserDisplayName: Alex Reynolds
Body: <p>Recently, our lab ran a 192 sample experiment through our ATACseq pipeline. In doing so, HOMER's <code>mergePeaks</code> 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.</p>
  <p>It seems that <a href="https://bedtools.readthedocs.io/en/latest/content/tools/merge.html" rel="nofollow noreferrer"><code>bedTools merge</code></a> 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 <code>MACS2</code>. The <code>mergeBed</code> only takes in 1 input file that must be sorted by chromosome and peak start position before we run it. I can run a <code>cat</code> call on all of the samples followed by a <code>sort</code>, but those individual sample bed files are already sorted. Is there a way to run <code>mergeBed</code> with multiple input files at once? I tried <code>mergeBed -i 1.bed -i 2.bed &gt; out.bed</code>, but that is equivalent to <code>mergeBed -i 1.bed &gt; out.bed</code>.</p>
  <p>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: <a href="https://i.sstatic.net/HVrE5.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/HVrE5.png" alt="mergeBed" /></a></p>
  <p>Additionally, bedtools does have a poorly documented function called <code>multiIntersectBed</code>that appears to run the <code>cat</code> and <code>sort</code> together, but is both slower and doesn't merge bookended sections--eg files look like this:</p>
  <pre><code>chr1 4658186 4658322
  chr1 4658322 4658775
  chr1 4658775 4658777
  chr1 4671049 4671351
  </code></pre>
  <p>When they should look like this:</p>
  <pre><code>chr1 4658186 4658777 # This was 3 lines before
  chr1 4671049 4671351
  </code></pre>
  <p>In summary, I have sorted bed files <code>1.bed</code> and <code>2.bed</code> 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 <code>mergePeaks</code>? Thank you!</p>
  <p>EDIT: Here's the speed of my current solutions at n=192:</p>
  <pre><code>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 &gt; tmp1; sort &gt; tmp2; mergeBed &gt; 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 &gt; tmp ; mergeBed -i tmp ; rm tmp`
  
  real    0m17.472s
  user    0m15.251s
  sys     0m2.212s
  </code></pre>
  <p>So... Followup questions, why does <code>sort -m</code> take longer than <code>sort | cat</code> especially when the sys time is indeed smaller? And why does piping take more time than creating temporary files?</p>
  
Tags: <bedtools><homer>
Title: alternative to HOMER's mergePeaks
ThreadText: # alternative to HOMER's mergePeaks
  
  (Asked by: username_0 on 2021-06-17)
  
  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`](https://bedtools.readthedocs.io/en/latest/content/tools/merge.html) 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: [![mergeBed](https://i.sstatic.net/HVrE5.png)](https://i.sstatic.net/HVrE5.png)
  
  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?
  
  **Tags:** bedtools, homer
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-06-21. Score: 1 | Accepted answer
  
  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](https://serverfault.com/questions/48717/practical-maximum-open-file-descriptors-ulimit-n-for-a-high-volume-system) 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>
  
  
  ---
  
  