Register FAQ Members List Calendar Search Today's Posts Mark Forums Read

 Similar Threads Thread Thread Starter Forum Replies Last Post Newsbot! Literature Watch 1 04-11-2012 01:18 PM liu_xt005 Bioinformatics 12 10-06-2011 08:14 AM jjk Bioinformatics 6 04-13-2011 04:32 PM sheremey Bioinformatics 3 11-02-2010 11:24 AM dukevn Bioinformatics 15 08-15-2010 09:51 AM

 02-15-2012, 11:58 AM #1 bioinf newbie Member   Location: India Join Date: Feb 2012 Posts: 13 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.
 02-15-2012, 12:14 PM #2 twaddlac Member   Location: Pittsburgh, PA Join Date: Feb 2011 Posts: 49 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.
 02-15-2012, 01:26 PM #3 Richard Finney Senior Member   Location: bethesda Join Date: Feb 2009 Posts: 700 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.
 02-15-2012, 08:36 PM #4 bioinf newbie Member   Location: India Join Date: Feb 2012 Posts: 13 Thank you very much Richard and @twaddlac.
 01-30-2013, 03:02 AM #5 ppoudel Junior Member   Location: UK Join Date: Feb 2012 Posts: 5 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?
 01-30-2013, 03:34 AM #6 Danielbenitezr Junior Member   Location: Concepcion, Chile Join Date: Dec 2012 Posts: 9 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!

 Tags bedtools bam, depth, genome coverage