[emboss-dev] Bug report and patch - BAM quality score reading

Peter biopython at maubp.freeserve.co.uk
Mon Aug 2 11:37:35 UTC 2010

Hi all,

Since I had several queries about what EMBOSS would do with SAM/BAM,
I decided to try it and see. I believe I have found a bug with reading quality
scores from BAM files in EMBOSS 6.3.1

$ seqret -version

I have been using a small pair of SAM and BAM files, originally downloaded
as a SAM file with reference FASTA sequence from the pysam project, which
I converted to BAM using samtools as they specify in their readme file


curl -O http://pysam.googlecode.com/hg/tests/ex1.fa
curl -O http://pysam.googlecode.com/hg/tests/ex1.sam.gz
gunzip ex1.sam.gz
samtools faidx ex1.fa
samtools import ex1.fa.fai ex1.sam ex1.bam

If we look at the first two reads in the SAM file, notice their quality strings:

$ head ex1.sam
EAS56_57:6:190:289:82	69	chr1	100	0	*	=	100	0	CTCAAGGTTGTTGCAAGGGGGTCTATGTGAACAAA	<<<7<<<;<<<<<<<<8;;<7;4<;<;;;;;94<;	MF:i:192
EAS56_57:6:190:289:82	137	chr1	100	73	35M	=	100	0	AGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCAC	<<<<<<;<<<<<<<<<<;<<;<<<<;8<6;9;;2;	MF:i:64	Aq:i:0	NM:i:0	UQ:i:0	H0:i:1	H1:i:0

Now let's ask EMBOSS seqret to convert from SAM to Sanger FASTQ,

$ seqret -sformat sam -osformat fastq-sanger ex1.sam -stdout -auto | head
@EAS56_57:6:190:289:82 chr1
@EAS56_57:6:190:289:82 chr1
@EAS51_64:3:190:727:308 chr1

The quality strings agree with the SAM file, good.

Now let's ask EMBOSS seqret to convert from BAM to Sanger FASTQ,

$ seqret -sformat bam -osformat fastq-sanger ex1.bam -stdout -auto | head

The quality strings differ, this is bad.

In the SAM file these two reads have quality strings starting the "<",
ASCII 60 meaning PHRED 60-33 = 27.

In the funny BAM to Sanger FASTQ conversion, EMBOSS has used
"]" which is ASCII 93, giving PHRED 93-33 = 60. i.e. 33 more than it
should be. I suspected that the EMBOSS code for reading BAM files
was wrongly applying a 33 offset to the quality scores. In BAM files
the scores are simply encoded directly as uint8_t without any offset.

Looking at the source code, file ajax/core/ajseqread.c we have:

        for(i=0; i < (ajuint) c->l_qseq; i++)
            ajFmtPrintAppS(&qualstr, " %02x", 33+d[dpos]);
            thys->Accuracy[i] = (float) (33 + d[dpos++]);

The creation of a quality string appears to be for debug only,
and here adding 33 to make it scores printable ASCII using
the Sanger FASTQ encoding makes sense. However, adding
the offset to the accuracy looks like an oversight. How about:

        for(i=0; i < (ajuint) c->l_qseq; i++)
            ajFmtPrintAppS(&qualstr, " %02x", 33+d[dpos]);
            thys->Accuracy[i] = (float) d[dpos++];

With this tiny change, I get the expected Sanger FASTQ output
from a BAM file using seqret.


Peter C.

More information about the emboss-dev mailing list