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.