SEQanswers coverage calculation for exome data
 Register FAQ Members List Calendar Search Today's Posts Mark Forums Read

 Similar Threads Thread Thread Starter Forum Replies Last Post apratap General 16 03-30-2020 06:32 AM kjaja Bioinformatics 0 12-15-2012 06:30 PM lyz1030 Bioinformatics 1 10-22-2012 03:26 AM pravee1216 Bioinformatics 10 04-12-2012 12:20 PM marlei Bioinformatics 2 06-03-2010 10:44 AM

 02-20-2013, 11:48 AM #1 kjaja Member   Location: NY Join Date: Aug 2011 Posts: 55 coverage calculation for exome data Hi, I am looking to calculate what percentage of all target positions achieved a coverage of at least 50x? can this be done? Appreciate your help
 02-20-2013, 11:51 AM #2 swbarnes2 Senior Member   Location: San Diego Join Date: May 2008 Posts: 912 Sure. Use BEDTools. You'll also want a .bed file of target regions.
 02-20-2013, 04:39 PM #3 kjaja Member   Location: NY Join Date: Aug 2011 Posts: 55 Hi Below is a sample output for the coverage file from bedtools, the coverage ranges from 0 to 8 for the same region (column 5), so I want to see what percentage of the exome was covered by at least 10x. it seems that for each region, I have to find the max number of bases covered by at least 10 x. is there an easy way to do this? chr1 176160618 176161057 Spna1-41 0 8 439 0.0182232 chr1 176160618 176161057 Spna1-41 1 59 439 0.1343964 chr1 176160618 176161057 Spna1-41 2 92 439 0.2095672 chr1 176160618 176161057 Spna1-41 3 103 439 0.2346241 chr1 176160618 176161057 Spna1-41 4 55 439 0.1252847 chr1 176160618 176161057 Spna1-41 5 21 439 0.0478360 chr1 176160618 176161057 Spna1-41 6 49 439 0.1116173 chr1 176160618 176161057 Spna1-41 7 27 439 0.0615034 chr1 176160618 176161057 Spna1-41 8 25 439 0.0569476
 02-20-2013, 05:15 PM #4 nexgengirl Member   Location: Maryland Join Date: Apr 2010 Posts: 31 Picard's hybridization statistics tool could be used. Description at the link below: http://picard.sourceforge.net/comman...ulateHsMetrics
 02-20-2013, 06:06 PM #5 rbagnall Member   Location: Sydney, Australia Join Date: Jun 2010 Posts: 34 I also use Bedtools. You need coverageBed with the -d option to show the coverage per base... Code: `samtools view -b bamfile.bam | coverageBed -abam stdin -b intervals.bed -d > per_base_coverage.txt` The file can be large as it has one row per base of the exome, but you can firstly use it to count the total number of bases in your exome if you are unsure.. Code: `wc -l per_base_coverage.txt` Then count the number of lines (i.e. bases) having at least 50 fold coverage.. Code: `awk '{FS="\t"}{if(\$5 > "49") print \$0}' per_base_coverage.txt | wc -l` It works for me, but there are could be quicker ways. Last edited by rbagnall; 02-21-2013 at 02:55 PM.
 02-21-2013, 07:04 AM #6 kjaja Member   Location: NY Join Date: Aug 2011 Posts: 55 THANKS! but there seem to be syntax error in the last code awk command? I am not sure how to fix it?
 02-21-2013, 02:58 PM #7 rbagnall Member   Location: Sydney, Australia Join Date: Jun 2010 Posts: 34 Hi kjaja, Sorry, there was an extra close parenthesis ) in the awk command. I have corrected it.
 03-23-2013, 09:56 PM #8 cuiya Junior Member   Location: beijing Join Date: Mar 2013 Posts: 1 right? Hi,rbagnall your answer is useful. but the last code may be "awk '{FS="\t"}{if(\$5 > 49) print \$0}' per_base_coverage.txt | wc -l". right?
 03-20-2014, 12:34 AM #9 woodydon Member   Location: https://www.researchgate.net/profile/Woody_Lin Join Date: Jan 2010 Posts: 52 Sorry I was wrong... Woody Last edited by woodydon; 03-20-2014 at 12:49 AM.