SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 02-20-2013, 11:48 AM   #1
kjaja
Member
 
Location: NY

Join Date: Aug 2011
Posts: 55
Default 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
kjaja is offline   Reply With Quote
Old 02-20-2013, 11:51 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Sure. Use BEDTools.

You'll also want a .bed file of target regions.
swbarnes2 is offline   Reply With Quote
Old 02-20-2013, 04:39 PM   #3
kjaja
Member
 
Location: NY

Join Date: Aug 2011
Posts: 55
Default

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
kjaja is offline   Reply With Quote
Old 02-20-2013, 05:15 PM   #4
nexgengirl
Member
 
Location: Maryland

Join Date: Apr 2010
Posts: 31
Default

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

http://picard.sourceforge.net/comman...ulateHsMetrics
nexgengirl is offline   Reply With Quote
Old 02-20-2013, 06:06 PM   #5
rbagnall
Member
 
Location: Sydney, Australia

Join Date: Jun 2010
Posts: 34
Default

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.
rbagnall is offline   Reply With Quote
Old 02-21-2013, 07:04 AM   #6
kjaja
Member
 
Location: NY

Join Date: Aug 2011
Posts: 55
Default

THANKS! but there seem to be syntax error in the last code awk command? I am not sure how to fix it?
kjaja is offline   Reply With Quote
Old 02-21-2013, 02:58 PM   #7
rbagnall
Member
 
Location: Sydney, Australia

Join Date: Jun 2010
Posts: 34
Default

Hi kjaja,

Sorry, there was an extra close parenthesis ) in the awk command. I have corrected it.
rbagnall is offline   Reply With Quote
Old 03-23-2013, 09:56 PM   #8
cuiya
Junior Member
 
Location: beijing

Join Date: Mar 2013
Posts: 1
Default 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?
cuiya is offline   Reply With Quote
Old 03-20-2014, 12:34 AM   #9
woodydon
Member
 
Location: https://www.researchgate.net/profile/Woody_Lin

Join Date: Jan 2010
Posts: 52
Default

Sorry I was wrong...

Woody

Last edited by woodydon; 03-20-2014 at 12:49 AM.
woodydon is offline   Reply With Quote
Reply

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 09:37 PM.


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