<p><em>Note: not yet tested, so there may be some additional fiddling with command line options needed</em></p>

<p>The per-base depth can be obtained from <em>samtools depth</em> (-a includes zero-coverage positions):</p>

<pre><code>samtools depth -a in1.bam &gt; depth_in1_both.tsv
</code></pre>

<p>To split this by forward and reverse, you can use an initial pipe through <em>samtools view</em> to exclude or include reverse-complement mappings:</p>

<pre><code>samtools view -b -F 0x10 in1.bam | samtools depth -a &gt; depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a &gt; depth_in1_rev.tsv
</code></pre>

<p>To filter for particular regions, this can be pre-piped through <em>bedtools intersect</em>:</p>

<pre><code>bedtools intersect -a in1.bam -b regions.bed | \
  samtools view -b -F 0x10 in1.bam | \
  samtools depth -a &gt; depth_in1_regions_fwd.tsv
</code></pre>

<p>Combining both, forward, and reverse is a bit trickier. I'd recommend doing it in a more manual fashion in R. But if you want to keep to the command line...</p>

<p>If you are absolutely sure that the 'depth' output will appear in the same order in all files (it might be worth using <code>samtools depth -aa</code>), then the lines can be <em>paste</em>d, with anonymous pipes, with unnecessary fields filtered out using <em>cut</em>:</p>

<pre><code>paste &lt;(samtools depth -aa in1.bam) \
      &lt;(samtools view -b -F 0x10 in1.bam | samtools depth -aa) \
      &lt;(samtools view -b -f 0x10 in1.bam | samtools depth -aa) | \
  cut -f 1,2,3,6,9 &gt; depth_in1_bfr.tsv
</code></pre>

<p>Or equivalently, if using intermediate files that have been generated already:</p>

<pre><code>paste depth_in1_both.tsv depth_in1_fwd.tsv depth_in1_rev.tsv |
  cut -f 1,2,3,6,9 &gt; depth_in1_bfr.tsv
</code></pre>

<p>It would also be possible to use <em>join</em> to combine the files (in cases where file output is <em>not</em> consistent), but the first two fields of each file would then need to be combined, because 'join' only joins on the first field... but if you need to get that complex, it'd be better to use something other than pipes on the command line.</p>

CC BY-SA 4.0
2019-05-17T23:49:00.650Z
2019-05-17T23:49:00.650Z
gringer
<pre><code>git clone https://github.com/lh3/htsbox
cd htsbox &amp;&amp; make
./htsbox pileup -cCf ref.fa aln.bam | less -S
</code></pre>

<p>This output a VCF containing positions covered by at least one read. Tag <code>ADF</code> and <code>ADR</code> give the forward/reverse depth for each allele. You can combine them to get the total counts on both strands.</p>

CC BY-SA 4.0
2019-05-18T02:38:35.590Z
2019-05-18T02:38:35.590Z
user172818
<p>Brent Pedersen (<a href="https://twitter/com/brentp" rel="nofollow noreferrer">@brentp</a>) claims <a href="https://github.com/brentp/mosdepth" rel="nofollow noreferrer">mosdepth</a> is 2x as fast as <code>samtools depth</code>. I'm not sure whether this provides strand information, but it might be the fastest tool available.</p>

<p>Another tool of possible interest (also written by Brent) is <a href="https://github.com/brentp/goleft/tree/master/indexcov#indexcov" rel="nofollow noreferrer">indexcov</a>. This doesn't give single bp precision, but is able to provide even quicker estimates of coverage at 16Kb resolution.</p>

<p>Installing these tools from source is non-trivial, and I struggled with installation in the past. However, they are both now available in bioconda, which should make installation simple.</p>

CC BY-SA 4.0
2019-05-21T15:52:52.357Z
2019-05-21T15:52:52.357Z
Daniel Standage
<p>I'm looking for a way to input a vcf or bed file (with specific base positions) and a bam file, and get the coverage at each base position (ie single base bins) using the bam file. I also want the strand information so ideal output would be something like:</p>

<pre><code>chr base_position total_coverage fwd_coverage rev_coverage
</code></pre>

<p>I've tried bamCoverage and bedtools Coverage but with bamCoverage, I can't figure out how to get single-base pair resolution and with bedtools Coverage, I can't figure out how to get the strand coverage information. </p>

<sam><bedtools><coverage>
How can I calculate coverage at single bases using a bam file?
# How can I calculate coverage at single bases using a bam file?

(Asked by: username_0 on 2019-05-17)

I'm looking for a way to input a vcf or bed file (with specific base positions) and a bam file, and get the coverage at each base position (ie single base bins) using the bam file. I also want the strand information so ideal output would be something like:
chr base_position total_coverage fwd_coverage rev_coverage

I've tried bamCoverage and bedtools Coverage but with bamCoverage, I can't figure out how to get single-base pair resolution and with bedtools Coverage, I can't figure out how to get the strand coverage information.

**Tags:** sam, bedtools, coverage

---

## 3 answers

### Answer by username_1 on 2019-05-17. Score: 5

Note: not yet tested, so there may be some additional fiddling with command line options needed
The per-base depth can be obtained from samtools depth (-a includes zero-coverage positions):
samtools depth -a in1.bam > depth_in1_both.tsv

To split this by forward and reverse, you can use an initial pipe through samtools view to exclude or include reverse-complement mappings:
samtools view -b -F 0x10 in1.bam | samtools depth -a > depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a > depth_in1_rev.tsv

To filter for particular regions, this can be pre-piped through bedtools intersect:
bedtools intersect -a in1.bam -b regions.bed | \
  samtools view -b -F 0x10 in1.bam | \
  samtools depth -a > depth_in1_regions_fwd.tsv

Combining both, forward, and reverse is a bit trickier. I'd recommend doing it in a more manual fashion in R. But if you want to keep to the command line...
If you are absolutely sure that the 'depth' output will appear in the same order in all files (it might be worth using samtools depth -aa), then the lines can be pasted, with anonymous pipes, with unnecessary fields filtered out using cut:
paste <(samtools depth -aa in1.bam) \
      <(samtools view -b -F 0x10 in1.bam | samtools depth -aa) \
      <(samtools view -b -f 0x10 in1.bam | samtools depth -aa) | \
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

Or equivalently, if using intermediate files that have been generated already:
paste depth_in1_both.tsv depth_in1_fwd.tsv depth_in1_rev.tsv |
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

It would also be possible to use join to combine the files (in cases where file output is not consistent), but the first two fields of each file would then need to be combined, because 'join' only joins on the first field... but if you need to get that complex, it'd be better to use something other than pipes on the command line.


---

### Answer by username_2 on 2019-05-18. Score: 4

git clone https://github.com/lh3/htsbox
cd htsbox && make
./htsbox pileup -cCf ref.fa aln.bam | less -S

This output a VCF containing positions covered by at least one read. Tag ADF and ADR give the forward/reverse depth for each allele. You can combine them to get the total counts on both strands.


---

### Answer by username_3 on 2019-05-21. Score: 1

Brent Pedersen ([@brentp](https://twitter/com/brentp)) claims [mosdepth](https://github.com/brentp/mosdepth) is 2x as fast as `samtools depth`. I'm not sure whether this provides strand information, but it might be the fastest tool available.

Another tool of possible interest (also written by Brent) is [indexcov](https://github.com/brentp/goleft/tree/master/indexcov#indexcov). This doesn't give single bp precision, but is able to provide even quicker estimates of coverage at 16Kb resolution.

Installing these tools from source is non-trivial, and I struggled with installation in the past. However, they are both now available in bioconda, which should make installation simple.


---

