Question: genomeCoverageBed for per nucleotide coverage
2
gravatar for Varun Gupta
15 months ago by
Varun Gupta1.1k
United States
Varun Gupta1.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 COMMENTlink modified 15 months ago by Kevin Blighe46k • written 15 months ago by Varun Gupta1.1k
3
gravatar for Kevin Blighe
15 months ago by
Kevin Blighe46k
Kevin Blighe46k 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 COMMENTlink written 15 months ago by Kevin Blighe46k

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 REPLYlink written 15 months ago by Varun Gupta1.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 REPLYlink written 15 months ago by Kevin Blighe46k
1

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 REPLYlink written 15 months ago by Varun Gupta1.1k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1582 users visited in the last hour