How can I calculate coverage at single bases using a bam file?

Question

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.

Answer 1

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 2

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 3

Brent Pedersen (@brentp) claims 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. 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.