samtools: BAQ differs between mpileup and calmd

[adeirossi]

Hello,

Though samtools is a wonderful tool, the other day I came across a behavior that left me confused. I have found some reads in an Illumina dataset where the BAQ calculation appears to give a different result depending on whether it is performed by samtools mpileup or samtools calmd. This contradicts my assumption that BAQ-adjusted base qualities would be the same regardless of whether mpileup or calmd is used. (I am running samtools version 0.1.17 on Linux.)

To demonstrate, I put one such read into its own BAM file and ran mpileup and calmd:

$ samtools view my.bam
SEQUENCER02_165:1:1207:8548:44865 99 chr3.fa 178936033 254 50M = 178936190 207 AAATGACAAAGAACAGCTCAAAGCAATTTCTACACGAGATCCTCTCTCTG CCCFFFFFHHHHHJJJJJJIJJJIJJJJJJJJJJJJJJJJIJIJJJJJJJ BC:Z:CGATGT XD:Z:50 SM:i:3 AS:i:406
## Here's the original base quality string.

$ samtools calmd -Ar my.bam hg19_sorted.bam.fa
@HD VN:1.0 SO:coordinate
...
SEQUENCER02_165:1:1207:8548:44865 99 chr3.fa 178936033 254 50M = 178936190 207 AAATGACAAAGAACAGCTCAAAGCAATTTCTACACGAGATCCTCTCTCTG ;;>FFFFFHHHHHJJJJJJIJJJIJJJJJJJJJJJJJJJJIJIJJJJJJE BC:Z:CGATGT XD:Z:50 SM:i:3 AS:i:406 ZQ:Z:HHE@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@E NM:i:0 MD:Z:50
## Here's the BAQ-adjusted base quality string, using calmd.

$ samtools mpileup -A -Q 0 -d 10000000 -f hg19_sorted.bam.fa my.bam
[mpileup] 1 samples in 1 input files
(mpileup) Max depth is above 1M. Potential memory hog!
chr3.fa 178936033 A 1 ^~. C
chr3.fa 178936034 A 1 . C
chr3.fa 178936035 A 1 . C
chr3.fa 178936036 T 1 . F
chr3.fa 178936037 G 1 . F
chr3.fa 178936038 A 1 . F
chr3.fa 178936039 C 1 . F
chr3.fa 178936040 A 1 . F
chr3.fa 178936041 A 1 . H
chr3.fa 178936042 A 1 . H
chr3.fa 178936043 G 1 . H
chr3.fa 178936044 A 1 . H
chr3.fa 178936045 A 1 . H
chr3.fa 178936046 C 1 . J
chr3.fa 178936047 A 1 . J
chr3.fa 178936048 G 1 . J
chr3.fa 178936049 C 1 . J
chr3.fa 178936050 T 1 . J
chr3.fa 178936051 C 1 . J
chr3.fa 178936052 A 1 . I
chr3.fa 178936053 A 1 . J
chr3.fa 178936054 A 1 . J
chr3.fa 178936055 G 1 . J
chr3.fa 178936056 C 1 . I
chr3.fa 178936057 A 1 . J
chr3.fa 178936058 A 1 . J
chr3.fa 178936059 T 1 . J
chr3.fa 178936060 T 1 . J
chr3.fa 178936061 T 1 . J
chr3.fa 178936062 C 1 . J
chr3.fa 178936063 T 1 . J
chr3.fa 178936064 A 1 . J
chr3.fa 178936065 C 1 . J
chr3.fa 178936066 A 1 . J
chr3.fa 178936067 C 1 . J
chr3.fa 178936068 G 1 . J
chr3.fa 178936069 A 1 . J
chr3.fa 178936070 G 1 . J
chr3.fa 178936071 A 1 . J
chr3.fa 178936072 T 1 . J
chr3.fa 178936073 C 1 . I
chr3.fa 178936074 C 1 . J
chr3.fa 178936075 T 1 . I
chr3.fa 178936076 C 1 . J
chr3.fa 178936077 T 1 . J
chr3.fa 178936078 C 1 . J
chr3.fa 178936079 T 1 . J
chr3.fa 178936080 C 1 . J
chr3.fa 178936081 T 1 . J
chr3.fa 178936082 G 1 .$ J
## Here's the BAQ-adjusted base quality string, using mpileup.
## Note that calmd and mpileup give different BAQ-adjusted base qualities at the start and end of the read.


None of the following modifications to the options seems to make any difference in this case:

• omitting the -A option with calmd
• using the default -Q value with mpileup
• using the -E option with mpileup

Why do the BAQ-adjusted base qualities differ between mpileup and calmd? Any thoughts or explanation would be greatly appreciated. Thanks in advance!

[kjohnsen]

I haven't conducted as thorough an investigation of the matter as adeirossi, but six years later I can confirm that in the pipeline we're developing, using calmd -bAr beforehand is not yielding the same results as performing the BAQ calculation in mpileup.