This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get average quality score per position in an aligned bam file

I have a bam file of nanopore reads aligned to a reference. I would like to calculate the average quality score of all bases mapping to every position in the reference. I thought that it should be easy to find a function for this in samtools or some other package, but apparently not.

I am aware of samtooks mpileup, but it only outputs a string of ASCII characters for each position. Is there not a package somewhere that can simply take my bam and give me the average quality score at each position without me having to write a whole script for it?

Thank you!

nanopore samtools bam phred33 quality
• 185 views • link
updated by GenoMax GenoMax • written by cedric.blais • cedric.blais

Not sure why you are looking to extract this information but assuming AI is not hallucinating this could be done by the the following two ways:

pysamstats --type baseq your_alignment.bam --genome /path/to/reference.fa > base_quality_per_pos.txt
  • Using samtools mpileup and awk :

samtools mpileup -f reference.fa your_alignment.bam | awk '{
    if ($5 ~ /^[A-Za-z]+$/) { 
        split($6, q, ""); 
        sum=0; 
        for(i=1; i<=length(q); i++) sum+=(ord[q[i]]-33); 
        print $1, $2, sum/length(q); 
    }
}'
• 0 views • link
written by GenoMax GenoMax

I'll give those a try, thank you!

As for why, I suspect there are systemic differences in the read quality supporting specific parts of the genome I'm working on, and I'd like to be able to precisely map/quantify them.

• 0 views • link
written by cedric.blais • cedric.blais

1 answer

Two adds to GenoMax's answer.

  1. The awk one-liner as written returns a constant negative value at every position: base awk has no ord() function, so ord[q[i]] is an undefined array element (0), and each term becomes 0 - 33 = -33. You have to build the ASCII-to-int map yourself in a BEGIN block:
    samtools mpileup -B -f reference.fa aln.bam | awk 'BEGIN{for(i=33;i<=126;i++)ord[sprintf("%c",i)]=i} $5 ~ /^[A-Za-z]+$/ {split($6,q,""); s=0; for(i=1;i<=length(q);i++) s+=ord[q[i]]-33; print $1,$2,s/length(q)}'
  1. Two mpileup gotchas that matter for a per-position quality map:

  • By default mpileup applies BAQ, which recomputes base qualities near indels -- usually the opposite of what you want when you are trying to see the raw quality supporting each position. Pass -B (--no-BAQ), as above, to keep the original Phred values.

  • The read-base column ($5) carries markers (^, $, +N, -N, *) for read starts/ends and indels, so the /^[A-Za-z]+$/ filter silently skips any position containing them, and indel-adjacent positions drop out. If you need those too, strip the markers rather than skip the line.

On pysamstats: --type baseq reports rms_baseq (root-mean-square), not the arithmetic mean. For Phred values the RMS sits a bit above the mean, so if you specifically want the mean, the mpileup route above is the one to use.

• 0 views • link
updated by GenoMax GenoMax • written by SeqBench • SeqBench

Log in to answer this question.

More posts like this