This is a test version of Biostars. For the public version, visit https://www.biostars.org.
genomeCoverageBed for per nucleotide coverage

Hello,

I want to get per nucleotide coverage for my particular gene (geneX) which I added in the fasta file. There are other chromosomes present in the bam file, whose per nu coverage I am not interested in. When I use genomeCoverageBed tool to calculate per nucleotide coverage for my gene I also get coverage for other chromosomes which should not be the case if you look at my command below

samtools view -bh test.bam geneX | genomeCoverageBed - ibam stdin -g my.genome -d > out.txt

my my.genome looks like this

geneX 1400

This runs for a very long time because the program still looks for chr1, chr2 etc??

Anyway I can only get per nuc counts for geneX.

Thanks

Varun

bedtools
• 3,444 views • link
updated by Kevin Blighe Kevin Blighe • written by Varun Gupta Varun Gupta

1 answer

You only need to use BEDTools coverage and will need your aligned BAM and a simple BED file with the co-ordinates for your gene (in BED or GTF/GFF format); Specifically, you need the -d flag of bedtools coverage:

cat mygene.bed
chr5  10000 23456

[let's assume that your gene is on chr5, between positions 10000 and 23456]

bedtools coverage -ibam test.bam -b mygene.bed -d

Similar question here: Extracting coverage per base pair

Kevin

• 1 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin, Thanks for your reply. I looked at the bedtools manual for coverage function. If I want per nucleotide coverage for the bed file(as is my case), the command should be like this

bedtools coverage -a mygenes.bed -b test.bam -d

Then this is the output

chr5    10000   23456   10001   14
chr5    10000   23456   10002   16
chr5    10000   23456   10003   21
chr5    10000   23456   10004   22
chr5    10000   23456   10005   23
chr5    10000   23456   10006   24

Thanks

• 0 views • link
written by Varun Gupta Varun Gupta

Yes, is that what you needed? The 4th column is the exact base number, and the 5th is the read-depth at that position. There should be 13,456 lines in your file.

• 0 views • link
written by Kevin Blighe Kevin Blighe

yup, works fine, must say better than genomeCoverageBed because in genomeCoverageBed, if you have all the headers(having all the chromosomes) in the bam file and only a particular chr alignments you want, it still goes through all the chromosomes.

• 0 views • link
written by Varun Gupta Varun Gupta

Log in to answer this question.

More posts like this