Samtools/vcftools related question

[Pradhaun]

Hello all,

Can anybody help me with a script parsing the vcftools output file. I have a bed file with start and end regions for which I need the depth of each region from bam file. Therefore, I used samtools mpileup, bcftools and vcftools. Now the output file of vcftools gives me depth in single nucleotide level. Though it is useful to know the depth per base, it would be more appropriate for my analysis to know the total, maximum and avg number of reads aligned to each BED region.

My bed file:
chr_1 5325 5500
chr_1 10909 11000
chr_1 12300 12598
.......
So, for the first bed region, the VCFtool output gives depth per base as shown below and it produce the same for other bed regions (it is a very huge file):
chr_1 5325 1
chr_1 5326 1
chr_1 5327 2
chr_1 5328 2
chr_1 5329 3
.....
chr_1 5499 1
chr_1 5500 1

And I would appreciate if anybody can help me with any script which outputs the total, maximum and average number of reads aligned to each bed region.

Thanks!!

[lindenb]

Don't reinvent the wheel: try something like gatk https://www.broadinstitute.org/gatk/...OfCoverage.php

[Pradhaun]

Thanks lindenb