Calculating number of reads per target base

[Lilach]

Hi,
I have results of whole exome NGS (SureSelect, enrichment by hybridization) as fastq files and I want to get the the list of target sequences that are not covered by at least x reads (for example, 10 reads).
Do anyone know how to do that?

Even if there is a software that can give me for each base in the target bases how many reads cover it (consensus+variant) - it will be excellent.
I have the list of targets as positions (i.e. chrn:1234568-4564567)

I can do alignment with BWA, MAQ, etc.

I will be gratefull for any hint!
Lilach

[Heisman]

That's simple, align the data, remove duplicate reads, index the data with samtools index, and then run the GATK DepthOfCoverage program while specifying your target intervals. You can have it output the coverage for every single base; large file, but then you can easily use an awk command to keep only the bases with under 10x coverage, which will hopefully be a much smaller file, and then you can take those chromosomes and bases, add one more column with the base changed by one position (basically, make it a bed file) and then use the "merge" feature in bedtools to combine everything. Then you'll have all of the intervals (and individual bases) where there is less than 10x coverage.

[Lilach]

Thank you very much!

will try that