Error-rate of a subset of the genome

Question

[skanterakis]

Does anyone know of an efficient way to calculate error-rate (% bases that mis-match to the reference) from a bam, in a subset of the genome defined by a bed file? By efficient I was hoping for a c/htslib based method. Thank you!

[GenoMax]

igvtools will allow you to get a count of bases present at each site from a BAM file. You can subset your BAM before doing that.

[skanterakis]

as far as I can see igvtools count --bases counts the occurrence of each base per region. Are you suggesting to then compare to the reference base with another command?

[h.mon]

Do you want to calculate the sequencing error rate, or naturally occurring polymorphisms (SNPs, indels, and so on)?

[skanterakis]

The sequencing error rate. I don't mind polymorphisms being included in the calculation though.

Answer 1

[skanterakis]

looks like the 5th column of samtools mpileup -f hg38.fa the.bam could be used or this: https://github.com/genome/bam-readcount