SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
PubMed: What can you do with 0.1x genome coverage? A case study based on a genome sur Newsbot! Literature Watch 1 04-11-2012 01:18 PM
How to obtain coverage % of the genome? liu_xt005 Bioinformatics 12 10-06-2011 08:14 AM
Genome coverage jjk Bioinformatics 6 04-13-2011 04:32 PM
"nucleotide coverage" to genome feature coverage sheremey Bioinformatics 3 11-02-2010 11:24 AM
genome (%) coverage - what people usually get? dukevn Bioinformatics 15 08-15-2010 09:51 AM

Reply
 
Thread Tools
Old 02-15-2012, 11:58 AM   #1
bioinf newbie
Member
 
Location: India

Join Date: Feb 2012
Posts: 13
Default Genome coverage

How do you calculate the genome coverage for a sequenced genome?
I have the BAM file for my genome, and I have calculated the read depth at each position using samtools depth command.
Is there a way I can calculate genome coverage using this?

Last edited by bioinf newbie; 02-15-2012 at 12:01 PM.
bioinf newbie is offline   Reply With Quote
Old 02-15-2012, 12:14 PM   #2
twaddlac
Member
 
Location: Pittsburgh, PA

Join Date: Feb 2011
Posts: 49
Default

If I'm not mistaken, the output of 'samtools depth' provides the index of the reference (second column) and the depth of that index (third column). So every reference index that has a depth greater than 0 reads (or whatever your cutoff value is) you can say that the base is covered. A simple perl script could do the trick or you could pipe it through the command line.
twaddlac is offline   Reply With Quote
Old 02-15-2012, 01:26 PM   #3
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 700
Default calculate coverage from bam file

samtools depth whole.bam | awk '{sum+=$3;cnt++}END{print sum/cnt" "sum}'

It outputs 2 numbers: average coverage for a coveraged base (sum/cnt) AND total base coverage (sum).

Note, some bases may not be covered.

Divide (total base coverage) by (genome size).

example : /samtools depth TCGA-AG-4008-01A-01D-1115-02_IlluminaHiSeq-DNASeq_whole.bam | awk '{sum+=$3;cnt++}END{print sum/cnt" "sum}
3.86809 10573174269

10573174269/3100000000 = 3.41070137709677419354 (note 3.1 billion is aproxmate size of human genome)

so 3.4 coverage
____________
Sanity check ....
A way to check this, and an faster way to calculate it is ... run samtools flagstat and get the mapped reads. Multiply this by read size then divide by genome size.

Example:
Get number of mapped reads ...
samtools flagstat TCGA-AG-4008-01A-01D-1115-02_IlluminaHiSeq-DNASeq_whole.bam | grep mapped
213309151 + 0 mapped (87.77%:nan%)
208629507 + 0 with itself and mate mapped
3785276 + 0 with mate mapped to a different chr

213309151 is what we're after = number of mapped reads

Find out read length
samtools view TCGA-AG-4008-01A-01D-1115-02_IlluminaHiSeq-DNASeq_whole.bam | awk '{print length($10)}' | head -1
50
(read lentgth is 50, beware there can be mixed read lengths)

So ...
(213309151*50)/3100000000 = 3.44047017741935483870

Again, 3.1 billion is the genome size of the fearsome creature "homo sapiens".

check. close enough. we good.
Richard Finney is offline   Reply With Quote
Old 02-15-2012, 08:36 PM   #4
bioinf newbie
Member
 
Location: India

Join Date: Feb 2012
Posts: 13
Default

Thank you very much Richard and @twaddlac.
bioinf newbie is offline   Reply With Quote
Old 01-30-2013, 03:02 AM   #5
ppoudel
Junior Member
 
Location: UK

Join Date: Feb 2012
Posts: 5
Default

Hi,

I have a bedgraph file, which goes like this

chr start end coverage

1 213 214 890
2 214 215 900
2 340 345 900
4 34090 34809 34

Now I would like to calculate the average coverage, is there any better way to do it using awk or any other tools?
ppoudel is offline   Reply With Quote
Old 01-30-2013, 03:34 AM   #6
Danielbenitezr
Junior Member
 
Location: Concepcion, Chile

Join Date: Dec 2012
Posts: 9
Default

Hey:
I have a question. All these applies only for genome assemblies?
Can you calculate the coverage for an "assembled" Transcriptome with the same command line that Richard Finney posted here?

I'm trying to assess transcriptome assemblies, looking for the best algorithm, at least, for my data. I used Trinity and MIRA. MIRA gives an output with a lot of stats, but Trinity doesn´t. So, in order to compare them I need to retrieve stats from Trinity outputs. Someone knows how to do that?

Also, I don't understand why Samtools needs a mapping information. Can you assess Quality or stats of the indexed reference with samtools? I mean, with the *.fasta.fai file.


Thanks!
Danielbenitezr is offline   Reply With Quote
Reply

Tags
bedtools bam, depth, genome coverage

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:49 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO