Answers.Body: <p>Just for reference, there is a not very in depth, but first-hand reasoning (Jim Kent) given here about why bedGraphToBigWig does not support streaming <a href="http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin" rel="nofollow noreferrer">http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin</a></p>
  
  <p>Edit: as noted wigToBigWig allows bedgraph streaming input via stdin but takes (possibly much) more memory than bedGraphToBigWig itself does</p>
  
  <pre><code>cat file.bedgraph | wigToBigWig stdin chrom.sizes output.bigwig
  </code></pre>
  
  <p>Note: You can also use <code>wiggletools bin 10</code> for example to bin rows and downsample</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-07T21:16:54.647Z
Answers.Text: Ah! "However if you do have enough memory, you can still use wigToBigWig on bedGraphs, and it will fit in a pipe."
Answers.UserDisplayName: Nathan S. Watson-Haigh
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-07T21:17:28.497Z
Answers.Text: That's sort of a weird decision. The zoom levels are written after the raw entries, so one can write all of the raw interval entries to disk and then do another pass over them when making the zoom levels. Then you don't need to store everything in memory. That's how libBigWig works, in fact.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-07T21:13:05.843Z
Answers.LastActivityDate: 2017-08-09T16:13:38.660Z
Answers.LastEditDate: 2017-08-09T16:13:38.660Z
Answers.UserDisplayName: Colin D
Answers.Body: <p>The <a href="http://bioconductor.org/packages/release/bioc/html/QuasR.html" rel="nofollow noreferrer">QuasR</a> package in R has a useful function called <code>qExportWig</code>. </p>
  
  <p>To use it you first need to create a sample sheet and a <code>qProject</code> object - this can be done starting with bam files or fastq files (in which case it also  performs alignments). You can then create wig or bigWig files without writing a bedGraph file to disk.</p>
  
  <p>An example would look something like this:</p>
  
  <pre><code>library(QuasR)
  
  # Create qProject object using qAlign
  # Sample file should contain tab-separated file names (BAM or fastq) and sample names
  myproj &lt;- qAlign("samples.txt", genome=mygenome, paired="no")
  
  # Export all samples in project as bigWig, scaling to mean number of alignments by default
  qExportWig(myproj, binsize=50, createBigWig=T)
  </code></pre>
  
  <p>It has a lot of flexible options, including scaling or log-transforming counts, shifting reads and changing the bin size.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-07T13:55:28.977Z
Answers.Text: It would be helpful to provide some example code to demonstrate how one can go from BAM to BigWig. Like the answer provided by @story
Answers.UserDisplayName: Nathan S. Watson-Haigh
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-07T12:07:52.207Z
Answers.LastActivityDate: 2017-06-08T06:56:00.910Z
Answers.LastEditDate: 2017-06-08T06:56:00.910Z
Answers.UserDisplayName: Sarah Carl
Answers.Body: <p>This can be done in R very easily from an indexed .bam file.</p>
  
  <p>Given single-end file for <strong>sample1</strong>.</p>
  
  <pre><code>library(GenomicAlignments)
  library(rtracklayer)
  
  ## read in BAM file (use readGAlignmentPairs for paired-end files)
  gr &lt;- readGAlignments('sample1.bam')
  
  ## convert to coverages
  gr.cov &lt;- coverage(gr)
  
  ## export as bigWig
  export.bw(gr.cov,'sample1.bigwig')
  </code></pre>
  
  <p>Be aware that this method doesn't include normalization steps (such as normalizing to total coverage). Most of these additional steps can be added if necessary. </p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T04:15:50.473Z
Answers.Text: I had hoped for a command line tool but I'm not totally against going into `R` if needs be. I need to test more extensively with regards to speed so I need to have more of a play. If tests show that it's slow for the number and sizes of BAM files I have, or a better command line answer materialises, I will accept this answer.
Answers.UserDisplayName: Nathan S. Watson-Haigh
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-07T13:32:48.943Z
Answers.LastActivityDate: 2017-06-07T13:32:48.943Z
Answers.UserDisplayName: story
Body: <p>I have a pipeline for generating a BigWig file from a BAM file:</p>
  
  <pre><code>BAM -&gt; BedGraph -&gt; BigWig
  </code></pre>
  
  <p>Which uses <code>bedtools genomecov</code> for the <code>BAM -&gt; BedGraph</code> part and <code>bedGraphToBigWig</code> for the <code>BedGraph -&gt; BigWig</code> part.</p>
  
  <p>The use of <code>bedGraphToBigWig</code> to create the <code>BigWig</code> file requires a BedGraph file to reside on disk in uncompressed form as it performs seeks. This is problematic for large genomes and variable coverage BAM files when there are more step changes/lines in the <code>BedGraph</code> file. My <code>BedGraph</code> files are in the order of 50 Gbytes in size and all that IO for 10-20 BAM files seems unnecessary.</p>
  
  <p><strong>Are there any tools capable generate <code>BigWig</code> without having to use an uncompressed <code>BedGraph</code> file on disk? I'd like this conversion to hapen as quickly as possible.</strong></p>
  
  <p>I have tried the following tools, but they still create/use a <code>BedGraph</code> intermediary file:</p>
  
  <ul>
  <li><code>deepTools</code></li>
  </ul>
  
  <h1>Some Benchmarks, Ignoring IO</h1>
  
  <p>Here are some timings I get for creating a <code>BigWig</code> file from a <code>BAM</code> file using 3 different pipelines. All files reside on a <code>tmpfs</code> i.e. in memory.</p>
  
  <h2>BEDTools and Kent Utils</h2>
  
  <p>This is the approach taken by most.</p>
  
  <pre><code>time $(bedtools genomecov -bg -ibam test.bam -split -scale 1.0 &gt; test.bedgraph \
    &amp;&amp; bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
    &amp;&amp; rm test.bedgraph)
  
  real    1m20.015s
  user    0m56.608s
  sys     0m27.271s
  </code></pre>
  
  <h2>SAMtools and Kent Utils</h2>
  
  <p>Replacing <code>bedtools genomecov</code> with <code>samtools depth</code> and a custom <code>awk</code> script (<a href="https://gist.github.com/nathanhaigh/c8322002a42e40c479d3d766aab32547" rel="noreferrer"><code>depth2bedgraph.awk</code></a>) to output <code>bedgraph</code> format has a significant performance improvement:</p>
  
  <pre><code>time $(samtools depth -Q 1 --reference test.fasta test.bam \
    | mawk -f depth2bedgraph.awk \
    &gt; test.bedgraph \
    &amp;&amp; bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
    &amp;&amp; rm test.bedgraph)
  
  real    0m28.765s
  user    0m44.999s
  sys     0m1.166s
  </code></pre>
  
  <p>Although it is has less features, we used <code>mawk</code> here as it's faster than <code>gawk</code> (we don't need those extra features here).</p>
  
  <h3>Parallelising with xargs</h3>
  
  <p>If you want to have a <code>BigWig</code> file per chromosome, you can easily parallelise this across chromosome/reference sequences. We can use <code>xargs</code> to run 5 parallel BAM->BedGraph->BigWig pipelines, each using the <code>tmpfs</code> mounted <code>/dev/shm</code> for the intermediary <code>BedGraph</code> files.</p>
  
  <pre><code>cut -f1 test.fasta.chrom.sizes \
    | xargs -I{} -P 5  bash -c 'mkdir /dev/shm/${1} \
    &amp;&amp; samtools depth -Q 1 --reference test.fasta -r "${1}" test.bam \
    | mawk -f scripts/depth2bedgraph.awk \
    &gt; "/dev/shm/${1}/test.bam.bedgraph" \
    &amp;&amp; mkdir "./${1}" \
    &amp;&amp; bedGraphToBigWig \
      "/dev/shm/${1}/test.bam.bedgraph" \
      test.fasta.chrom.sizes \
      "./${1}/test.bam.bw" \
    &amp;&amp; rm "/dev/shm/${1}/test.bam.bedgraph"' -- {}
  </code></pre>
  
  <h2>deepTools</h2>
  
  <p>Let's see how <code>deepTools</code> performs.</p>
  
  <pre><code>time bamCoverage --numberOfProcessors max \
    --minMappingQuality 1 \
    --bam test.bam --binSize 1 --skipNonCoveredRegions \
    --outFileName deeptools.bw
  
  real    0m40.077s
  user    3m56.032s
  sys     0m9.276s
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T09:41:33.003Z
Comments.Text: It's been on my to do list to remove the "make a single bedGraph file" intermediate step from deepTools. Given how that's all implemented there's no good way to avoid at least some intermediate file creation (without blowing up the memory at least).
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T10:25:52.487Z
Comments.Text: I needed to generate the intermediate bedgraph files when I was creating BigWig files for GBrowse in 2011; it's disappointing to see that it hasn't changed much in the many years since then. The file seeking is a killer -- it excludes most of the shell tricks that can be used to work around programs that don't explicitly support the common UNIX pipeline (e.g. ```<()```, FIFO buffers, ```/dev/stdin```).
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T10:33:13.367Z
Comments.Text: @DevonRyan when you say blow up memory, will that be equivilent to the size of the intermediary bedgraph file? I have a 720Gbyte fat memory compute node...so for me, speed is a bigger issue than memory. Perhaps an option to state the max memory usage, at which point you start going to disk?
Comments.UserDisplayName: Nathan S. Watson-Haigh
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T11:08:33.857Z
Comments.Text: @NathanS.Watson-Haigh Yeah, it'd be equivalent more or less to the file size. deepTools itself was never engineered to make this easier to implement and the generally crappy nature of multithreading in python doesn't help that :(
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T11:17:47.417Z
Comments.Text: @DevonRyan how about using the system's temp dir as the default location for the intermediary bedgraph file and expose this as a command line argument. That way I could use my raided SSD's for more convenient and faster IO. I was doing this with my bedtools and bedGraphToBigWig pipeline.
Comments.UserDisplayName: Nathan S. Watson-Haigh
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T11:28:03.147Z
Comments.Text: The system temp dir (or shared memory) is already what's used for all of the tmp files. Honestly, the current implementation is a hold-over from before I wrote pyBigWig (or even joined the group) and deepTools was using bedGraphToBigwig in a system call. It's mostly a matter of me finding time to change how bamCoverage does this final step.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T11:48:31.130Z
Comments.Text: @DevonRyan is there any reason a named pipe as suggested by my answer wouldn't work for this? I have never used these tools so I can't be sure.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T12:07:42.703Z
Comments.Text: @terdon: Sorry, I was referring specifically to deepTools, which was mentioned in the post. One can certainly do what you've described, it'll just be slower than what deeptools can do (but better with IO, so maybe it's all a wash at that point!).
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T15:12:52.817Z
Comments.Text: @terdon I don't see your named pipe based answer
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T15:14:46.817Z
Comments.Text: @bli no I deleted it since it won't work here because the question requires the ability to seek the file. FIFOs won't work with seeking AFAIK, so the answer wasn't useful. I hadn't realized the OP required a seekable file.
Comments.UserDisplayName: terdon
Tags: <sam><file-formats><format-conversion>
Title: BAM to BigWig without intermediary BedGraph
ThreadText: # BAM to BigWig without intermediary BedGraph
  
  (Asked by: username_0 on 2017-06-07)
  
  I have a pipeline for generating a BigWig file from a BAM file:
  
  ```
  BAM -> BedGraph -> BigWig
  
  ```
  
  Which uses `bedtools genomecov` for the `BAM -> BedGraph` part and `bedGraphToBigWig` for the `BedGraph -> BigWig` part.
  
  The use of `bedGraphToBigWig` to create the `BigWig` file requires a BedGraph file to reside on disk in uncompressed form as it performs seeks. This is problematic for large genomes and variable coverage BAM files when there are more step changes/lines in the `BedGraph` file. My `BedGraph` files are in the order of 50 Gbytes in size and all that IO for 10-20 BAM files seems unnecessary.
  
  **Are there any tools capable generate `BigWig` without having to use an uncompressed `BedGraph` file on disk? I'd like this conversion to hapen as quickly as possible.**
  
  I have tried the following tools, but they still create/use a `BedGraph` intermediary file:
  
  * `deepTools`
  
  # Some Benchmarks, Ignoring IO
  
  Here are some timings I get for creating a `BigWig` file from a `BAM` file using 3 different pipelines. All files reside on a `tmpfs` i.e. in memory.
  
  ## BEDTools and Kent Utils
  
  This is the approach taken by most.
  
  ```
  time $(bedtools genomecov -bg -ibam test.bam -split -scale 1.0 > test.bedgraph \
    && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
    && rm test.bedgraph)
  
  real    1m20.015s
  user    0m56.608s
  sys     0m27.271s
  
  ```
  
  ## SAMtools and Kent Utils
  
  Replacing `bedtools genomecov` with `samtools depth` and a custom `awk` script ([`depth2bedgraph.awk`](https://gist.github.com/nathanhaigh/c8322002a42e40c479d3d766aab32547)) to output `bedgraph` format has a significant performance improvement:
  
  ```
  time $(samtools depth -Q 1 --reference test.fasta test.bam \
    | mawk -f depth2bedgraph.awk \
    > test.bedgraph \
    && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
    && rm test.bedgraph)
  
  real    0m28.765s
  user    0m44.999s
  sys     0m1.166s
  
  ```
  
  Although it is has less features, we used `mawk` here as it's faster than `gawk` (we don't need those extra features here).
  
  ### Parallelising with xargs
  
  If you want to have a `BigWig` file per chromosome, you can easily parallelise this across chromosome/reference sequences. We can use `xargs` to run 5 parallel BAM->BedGraph->BigWig pipelines, each using the `tmpfs` mounted `/dev/shm` for the intermediary `BedGraph` files.
  
  ```
  cut -f1 test.fasta.chrom.sizes \
    | xargs -I{} -P 5  bash -c 'mkdir /dev/shm/${1} \
    && samtools depth -Q 1 --reference test.fasta -r "${1}" test.bam \
    | mawk -f scripts/depth2bedgraph.awk \
    > "/dev/shm/${1}/test.bam.bedgraph" \
    && mkdir "./${1}" \
    && bedGraphToBigWig \
      "/dev/shm/${1}/test.bam.bedgraph" \
      test.fasta.chrom.sizes \
      "./${1}/test.bam.bw" \
    && rm "/dev/shm/${1}/test.bam.bedgraph"' -- {}
  
  ```
  
  ## deepTools
  
  Let's see how `deepTools` performs.
  
  ```
  time bamCoverage --numberOfProcessors max \
    --minMappingQuality 1 \
    --bam test.bam --binSize 1 --skipNonCoveredRegions \
    --outFileName deeptools.bw
  
  real    0m40.077s
  user    3m56.032s
  sys     0m9.276s
  
  ```
  
  **Tags:** sam, file-formats, format-conversion
  
  > **Comments:**
  > **username_1**: It's been on my to do list to remove the "make a single bedGraph file" intermediate step from deepTools. Given how that's all implemented there's no good way to avoid at least some intermediate file creation (without blowing up the memory at least).  
  > **username_2**: I needed to generate the intermediate bedgraph files when I was creating BigWig files for GBrowse in 2011; it's disappointing to see that it hasn't changed much in the many years since then. The file seeking is a killer -- it excludes most of the shell tricks that can be used to work around programs that don't explicitly support the common UNIX pipeline (e.g. ```<()```, FIFO buffers, ```/dev/stdin```).  
  > **username_0**: @username_1 when you say blow up memory, will that be equivilent to the size of the intermediary bedgraph file? I have a 720Gbyte fat memory compute node...so for me, speed is a bigger issue than memory. Perhaps an option to state the max memory usage, at which point you start going to disk?  
  > **username_1**: @NathanS.username_0 Yeah, it'd be equivalent more or less to the file size. deepTools itself was never engineered to make this easier to implement and the generally crappy nature of multithreading in python doesn't help that :(  
  > **username_0**: @username_1 how about using the system's temp dir as the default location for the intermediary bedgraph file and expose this as a command line argument. That way I could use my raided SSD's for more convenient and faster IO. I was doing this with my bedtools and bedGraphToBigWig pipeline.  
  > **username_1**: The system temp dir (or shared memory) is already what's used for all of the tmp files. Honestly, the current implementation is a hold-over from before I wrote pyBigWig (or even joined the group) and deepTools was using bedGraphToBigwig in a system call. It's mostly a matter of me finding time to change how bamCoverage does this final step.  
  > **username_3**: @username_1 is there any reason a named pipe as suggested by my answer wouldn't work for this? I have never used these tools so I can't be sure.  
  > **username_1**: @username_3: Sorry, I was referring specifically to deepTools, which was mentioned in the post. One can certainly do what you've described, it'll just be slower than what deeptools can do (but better with IO, so maybe it's all a wash at that point!).  
  > **username_4**: @username_3 I don't see your named pipe based answer  
  > **username_3**: @username_4 no I deleted it since it won't work here because the question requires the ability to seek the file. FIFOs won't work with seeking AFAIK, so the answer wasn't useful. I hadn't realized the OP required a seekable file.  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_7 on 2017-06-07. Score: 11 | Accepted answer
  
  This can be done in R very easily from an indexed .bam file.
  
  Given single-end file for **sample1**.
  
  ```
  library(GenomicAlignments)
  library(rtracklayer)
  
  ## read in BAM file (use readGAlignmentPairs for paired-end files)
  gr <- readGAlignments('sample1.bam')
  
  ## convert to coverages
  gr.cov <- coverage(gr)
  
  ## export as bigWig
  export.bw(gr.cov,'sample1.bigwig')
  
  ```
  
  Be aware that this method doesn't include normalization steps (such as normalizing to total coverage). Most of these additional steps can be added if necessary.
  
  > **Comments:**
  > **username_0**: I had hoped for a command line tool but I'm not totally against going into `R` if needs be. I need to test more extensively with regards to speed so I need to have more of a play. If tests show that it's slow for the number and sizes of BAM files I have, or a better command line answer materialises, I will accept this answer.  
  
  ---
  
  ### Answer by username_5 on 2017-08-07. Score: 4
  
  Just for reference, there is a not very in depth, but first-hand reasoning (Jim Kent) given here about why bedGraphToBigWig does not support streaming <http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin>
  
  Edit: as noted wigToBigWig allows bedgraph streaming input via stdin but takes (possibly much) more memory than bedGraphToBigWig itself does
  
  ```
  cat file.bedgraph | wigToBigWig stdin chrom.sizes output.bigwig
  
  ```
  
  Note: You can also use `wiggletools bin 10` for example to bin rows and downsample
  
  > **Comments:**
  > **username_0**: Ah! "However if you do have enough memory, you can still use wigToBigWig on bedGraphs, and it will fit in a pipe."  
  > **username_1**: That's sort of a weird decision. The zoom levels are written after the raw entries, so one can write all of the raw interval entries to disk and then do another pass over them when making the zoom levels. Then you don't need to store everything in memory. That's how libBigWig works, in fact.  
  
  ---
  
  ### Answer by username_6 on 2017-06-07. Score: 2
  
  The [QuasR](http://bioconductor.org/packages/release/bioc/html/QuasR.html) package in R has a useful function called `qExportWig`.
  
  To use it you first need to create a sample sheet and a `qProject` object - this can be done starting with bam files or fastq files (in which case it also performs alignments). You can then create wig or bigWig files without writing a bedGraph file to disk.
  
  An example would look something like this:
  
  ```
  library(QuasR)
  
  # Create qProject object using qAlign
  # Sample file should contain tab-separated file names (BAM or fastq) and sample names
  myproj <- qAlign("samples.txt", genome=mygenome, paired="no")
  
  # Export all samples in project as bigWig, scaling to mean number of alignments by default
  qExportWig(myproj, binsize=50, createBigWig=T)
  
  ```
  
  It has a lot of flexible options, including scaling or log-transforming counts, shifting reads and changing the bin size.
  
  > **Comments:**
  > **username_0**: It would be helpful to provide some example code to demonstrate how one can go from BAM to BigWig. Like the answer provided by @username_7  
  
  ---
  
  