coverage precentage of the hg19 using bowtie2

[papori]

Hi all,
i am using bowtie2 for mapping WGS reads to hg19.
I am trying to find out what is the coverage precentage of the hg19 using my reads.
For example:
if the genome siae is 300 bases and my reads cover 30 bases, so my genome coverage is 10%.
i am looking for the 10%...
Can i know it from bowtie2 output?
Any ideas?

Thanks,
Pap

[GenoMax]

While it probably won't give you the gross % value you are looking for "genomecov" from BEDTools would be a good starting point to experiment: http://bedtools.readthedocs.org/en/l...genomecov.html

[jgibbons1]

If I understand you correctly, you'd like to know the average coverage of your read set. Is that right?

You may consider using the samtools depth function. This calculates depth for every base of the reference. You can then take the sum of the depth value and divide by non 'N' base pairs in hg19.

This will give you a quick and dirty estimate.

samtools depth you.file.bam | awk '{sum=sum+$3} END {print sum}'

(this will give you the sum of depth values)

[jgibbons1]

If you just want to know how many reads are covered, you can use the same tool, as it does not report '0' depth values. You can simply subtract "covered" bases from total bases.

[papori]

If I understand you correctly, you'd like to know the average coverage of your read set. Is that right?

You may consider using the samtools depth function. This calculates depth for every base of the reference. You can then take the sum of the depth value and divide by non 'N' base pairs in hg19.

This will give you a quick and dirty estimate.

samtools depth you.file.bam | awk '{sum=sum+$3} END {print sum}'

(this will give you the sum of depth values)

I want to know how many bases from the reference has hit from my reads.
I dont care from the depth.
Is it clearer?
Thanks

[papori]

While it probably won't give you the gross % value you are looking for "genomecov" from BEDTools would be a good starting point to experiment: http://bedtools.readthedocs.org/en/l...genomecov.html

This function is great, but the output is huge!!! more than 2TB per bam file..
This is to much for my storage
Any solutions?

[dpryan]

This function is great, but the output is huge!!! more than 2TB per bam file..
This is to much for my storage
Any solutions?

Just don't store the output then. If you just want to know how many bases have a read aligning to them, just use the "-d" option and then pipe the output to awk:

bedtools genomecov ...stuff.. | awk '{if($3==0) { nocov+=1;}else{cov+=1}}END{printf("%i have hits and %i do not\n",cov,nocov)}'

That will occupy no space and be faster since you don't have to write anything to disk.

[papori]

Just don't store the output then. If you just want to know how many bases have a read aligning to them, just use the "-d" option and then pipe the output to awk:

bedtools genomecov ...stuff.. | awk '{if($3==0) { nocov+=1;}else{cov+=1}}END{printf("%i have hits and %i do not\n",cov,nocov)}'

That will occupy no space and be faster since you don't have to write anything to disk.

Is it report 1 for base which have also more than 1 hits?
I guess yes is the answer, just want to be sure..
Many Thanks!

[dpryan]

It does. A base with a million reads mapping to it counts the same as another base with only a single read covering it.

You'll find awk and the other standard unix tools extremely useful.

[papori]

It does. A base with a million reads mapping to it counts the same as another base with only a single read covering it.

You'll find awk and the other standard unix tools extremely useful.

Thanks!
Very Help

[SeqMlife]

Hi papori,
I am currently doing something similar to you.
Could you share one of the example final output from the code?
like cov= ,ncov=?
The result I got is extremely huge.
879370384 hits and 1984311365 no hits. Not sure if I was doing it right.
Thanks for your reply!