Extract coverage from each nucleotide position of a genome in a bam file
1
1
Entering edit mode
5.7 years ago
MAPK ★ 2.1k

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!

bam bed coverage • 2.0k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode

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
ADD REPLY
0
Entering edit mode

That's the default...

ADD REPLY
2
Entering edit mode
5.7 years ago
samtools depth -a -b in.bed in.bam
ADD COMMENT
0
Entering edit mode

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).

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 1936 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6