coverage calculation for exome data

kjaja · 02-20-2013, 11:48 AM

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

swbarnes2 · 02-20-2013, 11:51 AM

Sure. Use BEDTools.

You'll also want a .bed file of target regions.

kjaja · 02-20-2013, 04:39 PM

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

nexgengirl · 02-20-2013, 05:15 PM

Picard's hybridization statistics tool could be used. Description at the link below:

http://picard.sourceforge.net/comman...ulateHsMetrics

rbagnall · 02-20-2013, 06:06 PM

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.

kjaja · 02-21-2013, 07:04 AM

THANKS! but there seem to be syntax error in the last code awk command? I am not sure how to fix it?

rbagnall · 02-21-2013, 02:58 PM

Hi kjaja,

Sorry, there was an extra close parenthesis ) in the awk command. I have corrected it.

cuiya · 03-23-2013, 09:56 PM

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?

woodydon · 03-20-2014, 12:34 AM

Sorry I was wrong...

Woody

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Size of human transcriptome/exome for coverage calculation	apratap	General	16	03-30-2020 06:32 AM
power calculation for exome seq	kjaja	Bioinformatics	0	12-15-2012 06:30 PM
extra high coverage regions in exome sequencing data	lyz1030	Bioinformatics	1	10-22-2012 03:26 AM
Estimate on-target coverage for exome data	pravee1216	Bioinformatics	10	04-12-2012 12:20 PM
Coverage calculation (on GA data)	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, 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.

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.

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-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.