How Do You Get The Quality Score And Coverage For Every Single Position Of A Reference Assembly

Question

[Joseph Hughes]

Hi,

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.

Joseph

Answer 1

[Pierre Lindenbaum]

Don't use use bcftools with the "calling options" ( -v , -c -g ). You should get an output for all the positions

[Joseph Hughes]

Hi Pierre, I have tried that, there are more positions but not all of them. I am going to try setting -d10000000 in the mpileup command to see if that changes anything. Thanks