coverage calculation for exome data

[kjaja]

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]

Sure. Use BEDTools.

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

[kjaja]

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]

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

http://picard.sourceforge.net/command-line-overview.shtml#CalculateHsMetrics

[rbagnall]

I also use Bedtools. You need coverageBed with the -d option to show the coverage per base...

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

wc -l per_base_coverage.txt

Then count the number of lines (i.e. bases) having at least 50 fold coverage..

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]

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

[rbagnall]

Hi kjaja,

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

[cuiya]

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?

[woodydon]

Sorry I was wrong...

Woody