How can I calculate coverage at single bases using a bam file? 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.
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 &gt; 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 &gt; depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a &gt; 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 &gt; 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 &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


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 &gt; 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.