I have trouble to understand how bedtools computing coverage. I have vcf file (generated from Illumina somatic caller) - I was created bed file from my vcf.
Then I try to compute coverage with bedtools and I have for the same coordinate different coverage. Could you help me to explain why? Probably some reads where filter out during vcf calling algorithm.
coverage in vcf:
chr13 32890572 DP=45048 chr13 32893197 DP=1494 chr13 32899359 DP=12809 chr13 32900176 DP=57850 chr13 32900177 DP=61728
and in my bedtools output:
chr13 32890572 32890573 47857 1 1 1.0000000 chr13 32893197 32893198 1686 1 1 1.0000000 chr13 32899359 32899360 15673 1 1 1.0000000 chr13 32900176 32900177 65461 1 1 1.0000000 chr13 32900177 32900178 65461 1 1 1.0000000
Note, that my vcf coverage is before filtering.
Can anybody now why is this different? And if you can explain (mathematically) how bedtools computing coverage? Is it sum of all coverage in the interval and then divide length of interval? Of it is just sum of all coverage in interval?
Thank you so much for clarification.
Its important to know, how the coverage was calculated exactly, which sub-tool and which parameters were used. Just write the whole command used to generate the file.
Thank you for response. VCF was generate automatically from Illumina MiSeq Reporter - Illumina Somatic Variant Caller. And I am use coordinates of my variants to generate bed file to calculate bedtools coverage.
My code was:
I am not able to find how bedtools calculating coverage - mathematically meaning.
Did you read this, it should answer your question.