genomeCoverageBed for per nucleotide coverage

Question: genomeCoverageBed for per nucleotide coverage

15 months ago by

United States

Varun Gupta • 1.1k wrote:

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 • 771 views

ADD COMMENT • link •

modified 15 months ago by Kevin Blighe ♦ 46k • written 15 months ago by Varun Gupta • 1.1k

15 months ago by

Kevin Blighe ♦ 46k

Kevin Blighe ♦ 46k wrote:

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

ADD COMMENT • link written 15 months ago by Kevin Blighe ♦ 46k

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

ADD REPLY • link written 15 months ago by Varun Gupta • 1.1k

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.

ADD REPLY • link written 15 months ago by Kevin Blighe ♦ 46k

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.

ADD REPLY • link written 15 months ago by Varun Gupta • 1.1k

Similar posts • Search »