Extract coverage from each nucleotide position of a genome in a bam file

Question

[MAPK]

Hi All, I have a bed file called my_region.bed as below:

Ch1 0   3951982 my_region       0       +
Ch2 0   3683506 my_region       0       +
Ch3 0   3351453 my_region       0       +
Ch4 0   2873318 my_region       0       +
Ch5 0   2822964 my_region       0       +
Ch6 0   2483831 my_region       0       +
Ch7 0   2434682 my_region       0       +
Ch8 0   2299506 my_region       0       +
Ch9 0   2122865 my_region       0       +
Ch10    0   2105496 my_region       0       +
Ch11    0   2052242 my_region       0       +
Ch12    0   1878461 my_region       0       +
Ch13    0   1845946 my_region       0       +
Ch14    0   1815632 my_region       0       +
Ch15    0   1765292 my_region       0       +
Ch16    0   1419421 my_region       0       +

I tried to extract the coverage for each nucleotide position on these chromosomes from my bam file using these commands: for positive strand,

bedtools coverage -a my_region.bed -b file.bam -bed -d -s | gzip > coverage_positive.tsv.gz

and for negative strand bedtools coverage -a my_region.bed -b file.bam -bed -d -S | gzip > coverage_negative.tsv.gz

However, these commands only give me positions expaning all chromosomes without their actual coverage (everything is zero in coverage column). Also, these commands work perfectly fine when I have only one chromosome on my bed file, but doesn't seem to work if I have all chromosomes as shown above. Could someone please let me know what I need to do instead of these commands. Thanks!

[GenoMax]

mosdepth (https://github.com/brentp/mosdepth ) is recommended for this application.

[MAPK]

Thanks, but couldn't make this work on individual base positions.

[GenoMax]

What does your bed file look like? Did you try making the bed regions 1 bp in length, if you only have one bp locations?

$ more test.bed
chr14   30602680        30602781
chr14   30602682        30602783

generates

chr14   30602680        30602781        219.50
chr14   30602682        30602783        228.55

[WouterDeCoster]

That's the default...

Answer 1

[Pierre Lindenbaum]

samtools depth -a -b in.bed in.bam

[MAPK]

Thanks, does this also separate based on strand? Is there a way to sort the output according to bed file's chromosome order. My alignments in bam file are in different chromosomal order (eg. Ch16, Ch15, Ch14...Ch1).

[Pierre Lindenbaum]

Thanks, does this also separate based on strand?

one can use an upstream samtools view to filter the read on strand

Is there a way to sort the output according to bed file's chromosome order.

run a loop with the bed file for each chromosome.