I am trying to extract the coverage and the average quality score for each position of a reference assembly in bam/sam format. I have managed to get the coverage using BEDtools
genomeCoverageBed -ibam mybamfile.bam -g my_genome -d > my_coverage.txt
but am at a loss on how to get some measure of the quality of the base calls at each position. I was thinking that I could use the bcftools to get a variant call formatted file
samtools mpileup -uf ref.fa mybamfile.bam | bcftools view -bvcg - > var.raw.bcf bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf
but this only provides the sites for which there are SNPs. Any advice greatly appreciated.