SAM format. Strange encoding in QUAL field -no ASCII

[tomjan]

Hello,

I have found strange encoding in QUAL field in sam format file which I'm obtaining with command "samtools view file.bam Chr10:73-73 > tmp.sam"

unfortunately the forum edytor does not allow paste this strange encoded field but I have attached the file with example line. Please look at last thirty qualities.

As you see some part of bases quality are encoded properly.

when I'm decoding this with perl according to the site http://maq.sourceforge.net/fastq.shtml $Q = ord($q) - 33;

I'm able to decode only properly ASCI encoded qualities.
When trying decode qualities encoded with 'strange encoding' the result is '-3'

C 29
A 24
G 28
C 32
A 30
G -3
G -3
G -3
C -3

do you have any suggestion how to fix this problem?

Thak You,

tj

[maubp]

I would guess your BAM file is corrupt/invalid. Is is possible to reduce this to a small test BAM file? e.g. If you use 'samtools view' to extract a small region as BAM, for instance using "samtools view -b file.bam Chr10:73-73 > tmp.bam", does the problem persist on converting that to SAM? If so, could you share the small broken BAM file?

[tomjan]

hello maubp, thank you for your help.

>I would guess your BAM file is corrupt/invalid.

I was also thinking that file is corupted. the more because it shows me "[bam_header_read] EOF marker is absent. The input is probably truncated." sentence when it is loaded by samtools.

But when load it to samtools tview I can see that bases with strange encoding are differently interpreted by samtools. I mean, when I press 'b' I can see that bases has different colours. I dont know if it mean that encoding is correct - I just share my obserwations


>Is is possible to reduce this to a small test BAM file? e.g. If you use 'samtools view' to
> extract a small region as BAM, for instance using "samtools view -b file.bamChr10:73-73
> tmp.bam", does the problem persist on converting that to SAM?

yes, the problem still persists

>If so, could you share the small broken BAM file?

Of Course, please find attached zip file

Thank you again, be happy and healthy

tj

[maubp]

After a minor hassle with the Mac decompression utility being to clever by half and unzipping and then ungzipping automatically I could index and view the tmp.bam file.

I don't see anything amiss:

$ samtools index tmp.bam
$ samtools view -h tmp.bam
@SQ	SN:Chr10	LN:20055733
@SQ	SN:Chr11	LN:19490693
@SQ	SN:Chr12	LN:19992050
@SQ	SN:Chr13	LN:17705833
@SQ	SN:Chr14	LN:15238285
@SQ	SN:Chr15	LN:12988017
@SQ	SN:Chr16	LN:539312
@SQ	SN:Chr17	LN:10597884
@SQ	SN:Chr18	LN:10944493
@SQ	SN:Chr19	LN:10362646
@SQ	SN:Chr1	LN:196441523
@SQ	SN:Chr20	LN:14157587
@SQ	SN:Chr21	LN:6914634
@SQ	SN:Chr22	LN:4122003
@SQ	SN:Chr23	LN:5769937
@SQ	SN:Chr24	LN:6301418
@SQ	SN:Chr25	LN:2104337
@SQ	SN:Chr26	LN:5236184
@SQ	SN:Chr27	LN:5316963
@SQ	SN:Chr28	LN:4897722
@SQ	SN:Chr2	LN:148351000
@SQ	SN:Chr32	LN:2087
@SQ	SN:Chr38	LN:1256500
@SQ	SN:Chr39	LN:79566926
@SQ	SN:Chr3	LN:110372186
@SQ	SN:Chr40	LN:1004462
@SQ	SN:Chr41	LN:825617
@SQ	SN:Chr4	LN:90259632
@SQ	SN:Chr5	LN:59208505
@SQ	SN:Chr6	LN:35063351
@SQ	SN:Chr7	LN:36220071
@SQ	SN:Chr8	LN:28840261
@SQ	SN:Chr9	LN:23535872
@SQ	SN:chrZ_061609	LN:83952565
EBRI093151_0051:5:10:12239:4897#0	99	Chr10	10	0	101M	=	412	503	AGAAGCCTCATTTAGACGTGCCTACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGC	FFFFEFFFFFFFFFFFFFDFFFFFFEFEFFFFCFFFEDDFCEDFDFDFFEAEFFF@FDDFDFED>DEDDEDEADDDDDDDDDDD>DDDAD?DDDDDDDD:B
EBRI093151_0054:6:3:1763:12550#0	163	Chr10	14	29	101M	=	470	557	GCCTCATTTAGACGTGCCTACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGAAACTCAAGCTCTGCACCACAAACG	BBBBBABBABBBBB@AABBBBABBBBBBBBB:A?BB;B9B?:B>@AA;BA??B7?BBB?B?=BAA*1;:;)''74:*;@9@5?<>A;B?
EBRI093151_0054:6:113:12058:5288#0	99	Chr10	32	23	101M	=	512	581	TACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGA	BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCBCBABBB?BBBBBA@BAA@BBB>BBBBBA=@@:4BB?<B<<@@@<A?@B@
EBRI093151_0051:4:55:2998:9540#0	99	Chr10	33	16	101M	=	248	316	ACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGAT	FFFFFFFFFFDFFBFDEAEEEFFFFFFFFCFFEFFCEEEDDFFEEEFEADDFDFDEEDFFE@FCDDD>ACDFADD?CCECDB<?@047:9@?BB+B@@@
EBRI093151_0054:6:5:6793:14536#0	73	Chr10	34	0	101M	=	34	0	CAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGCCTTTGGACG	BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB?BBBBABBBBBBBBBBB?BA@@BBB=>B@9B@B@7B;:881>2>?
EBRI093151_0051:5:20:5981:18486#0	89	Chr10	35	0	11M1I89M	=	35	0	CACAGAGGGGGTTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATG	><2@>@>?><<>>4=9=:4>>9A<><<?3DECEFEDFFFFFFFFFFFFFFFFFFFFFF
EBRI093151_0051:4:68:7839:4415#0	99	Chr10	38	15	3S98M	=	235	298	CCAACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGG	FFFFFFFDFFDDFFFFFFFFFFFEFFDFFFFFFFFFDFEFFFFFFFFFFFFFFEFDFFFFFEEDAFDDDFCFFECEDDFA?@9BBBBBABB9@@BABAAAD
EBRI093151_0054:6:39:4855:15851#0	99	Chr10	38	15	27S74M	=	475	538	ACGATTGTATAAAACATTTCTACTCCAACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCC	BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB@BBBCBBBBBCBBBBBCBBBBBBBB??@BBAABBAC?BABA@AB8?<5-@?><BCAAC<BA@A<@
EBRI093151_0051:4:50:11602:8115#0	99	Chr10	39	9	101M	=	252	314	CAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGC	FFFFEFFFFFFFFFFFFFFFFFDFFFEEFFFFFFEDDEDCDEEFDDDD9DADFEFDDDDDDDFDDA@BB4@BB@CB<B69;4=8+55:5+>5>8;?@>>;8
EBRI093151_0051:5:20:4520:12245#0	163	Chr10	40	11	101M	=	316	377	AGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCA	FFF@FFFFFFFFCFFFCFFFFFFFFFFFFFFCFFFFFEFFFFFFFDFFDFDFFDFEFFDDDDFDEDDFDD?CDBD>B@BABE@BBB@=@@BB;9<?<+>
EBRI093151_0051:4:112:7929:1711#0	163	Chr10	45	23	101M	=	268	324	GTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGACGGCAGCAGGGCC	EE@EECEC@ED4DD>EEEFDEFBADFCF@EECDEDFCFEFFCDFEDFEEBFFFF>DF9DC@D?ADEACCD<@AA4,;>+<,=>9=B@@B4<DBBB>>;@78
EBRI093151_0051:4:17:10201:5574#0	163	Chr10	47	8	101M	=	259	313	TGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAACGGCAAGGCTTTGGACGGCAGCAGGGCCGG	A>4CC@>CCCCCDD9@@CC><C>>@DCDDDD<9D@@<B>4CC9CA9><@BD<9DCC999C66,A3B3,3B>5949=9B<B?<>
EBRI093151_0051:5:109:5445:11427#0	137	Chr10	49	0	101M	=	49	0	GAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAA	FFEFFFFFFFFFFFFFFEE>DCFFFFFAFFFFDFFFFFEFDEFEEDFEDECFFC@D@DDCDADDDDBDDBC?CAADDB?DBBB?3?9+?;>@==@
EBRI093151_0051:5:67:9954:9956#0	99	Chr10	51	20	101M	=	464	514	GTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGG	FFFFFFFFFFFFFFFFFFFFFFFFFFFDFFFFFFFEFFFFFFFFFFDFFCFFFFFFFFFEFCFF>BDABBEACEDC@BDCDC<>DBBDABC?B1@BC;@@A
EBRI093151_0051:4:83:6307:17773#0	99	Chr10	56	9	101M	=	264	309	GCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGGCATAA	FFFDFEFFFFFFFEFFFDFFFFFFFFFFFEFFEEFFCF@FFFEEFDB9DDEEDEEFEFECAECA>>>94A99@@@DCBCE?AA4@@?>5@>>
EBRI093151_0051:4:89:9204:12518#0	99	Chr10	56	10	101M	=	252	297	GCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAATGCATAA	EFEFFFEFFFEFFFFFFFFFFFFFFFFFFFE@FFFFFFFFFCECFCEFFEFFFEFFFFDFDFFD>ED>E9@ABBB@CDFAD9A?AB9BBB?,>>>>>B4>>
EBRI093151_0051:4:118:9782:16360#0	163	Chr10	61	14	101M	=	279	319	TCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGGCATAAAATGC	DD@>D><4CCC>><CDD9D4DDD@A>C4CACDCDC@E@EEEC@C9EEE>EADCC4,@?B4EEAEDCEECDBDCDDD,4D<@4A@?
EBRI093151_0051:4:82:14603:9124#0	163	Chr10	62	15	101M	=	282	321	CCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGGCATAAAATGAT	FFFFFFFFFFFFFFFFFFFFFFDFFFFCDECEFFFEFFDEFFFEFFFFFDDDDDADDDDCDDAC>D@DD?>BCCADD@DDDD9DDB;B>??@>@
EBRI093151_0051:5:107:14690:3406#0	99	Chr10	63	0	101M	=	581	619	CAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGGCATAAAATTCTC	FFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFFFEFFFFEFDDFFFFAFFFFFDDFFFFFECFCDDBFCDCCDEDCBDA>9>>>??+@<A>@9@=>>5?<@6
EBRI093151_0051:4:44:8923:12230#0	73	Chr10	71	0	101M	=	71	0	CTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGGCATAAAATGCTCTGTTTCGT	FFFFFFFFFFEFFFFFFFFFEFFFFFEFFFFFFFDFFFFFFFFFEFDFDDFFDFDCDDEDCDDDDDDDDDD@DADBBAAABC<B9@@>
EBRI093151_0054:6:31:8979:9911#0	73	Chr10	71	0	101M	=	71	0	CTATCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCATTGGACGGCAGCAGGGCCGGAAGGCATAAAATGCTCTGTTTCGT	>>%>>BBBBBBBCBBBBBBABBBBABBCB@BB@BBBCB@BACCC@A@AAACA>@@>)>?>,4174>9=A?
EBRI093151_0051:5:33:14948:6712#0	73	Chr10	72	0	101M	=	72	0	TTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGATGGCAGCAGGGCCGGAAGGCATAAAATGCTCTGTTTCGTG	FFFFFFFFFFFFFFFCFFFFFFFFEFFFEEEFFFFFFFDFFFFCDDFDDFDDDDDDDDFADBDDDDDDD@DDADDD;DDD>>@=@DDD9DAAA>AAAAA>@

That gave 22 reads, and by eye looks like a normal SAM file. Does it match what you get? Are any of these reads giving you a problem.

[tomjan]

please look at the length of QUAL string
there is 101 bases in each read, so should be 101 ASCI signs in QUAL string. because each one code quality amount for one base, is that correct?

if yes you will find that QUAL strings has uncorrect length.

for example in last one

'FFFFFFFFFFFFFFFCFFFFFFFFEFFFEEEFFFFFFFDFFFFCDDFDDFDDDDDDDDFADBDDDDDDD@DDADDD;DDD>>@=@DDD9DAAA>AAAAA>@"

there is 101 signs, everything is ok


but in one before last read

'>>%>>BBBBBBBCBBBBBBABBBBABBCB@BB@BBBCB@BACCC@A@AAACA>@@>)>?>,4174>9=A?'

you can see only about seventy signs. I can see (also working on mac) the rest thirty - please find attahed screen, look at the bottom. they appears as such 'strange rectangles' on the end.


tj

[maubp]

please look at the length of QUAL string
there is 101 bases in each read, so should be 101 ASCI signs in QUAL string. because each one code quality amount for one base, is that correct?

Yes, I see the problem now. Curious - I'll try to look into this tomorrow.

[maubp]

There does seem to be a problem in the BAM file itself, which samtools is not handling nicely.

The first read is OK, but the second is not - EBRI093151_0054:6:3:1763:12550#0 - this has mostly PHRED 33, at lowest PHRED 6, but the last dozen quality scores are 253 (unbelievably high), which are too high to represent in the Sanger FASTQ encoding used in SAM.

Look at the first two reads as SAM, we get:

$ samtools view tmp.bam | head -n 2
EBRI093151_0051:5:10:12239:4897#0	99	Chr10	10	0	101M	=	412	503	AGAAGCCTCATTTAGACGTGCCTACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGC	FFFFEFFFFFFFFFFFFFDFFFFFFEFEFFFFCFFFEDDFCEDFDFDFFEAEFFF@FDDFDFED>DEDDEDEADDDDDDDDDDD>DDDAD?DDDDDDDD:B
EBRI093151_0054:6:3:1763:12550#0	163	Chr10	14	29	101M	=	470	557	GCCTCATTTAGACGTGCCTACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGAAACTCAAGCTCTGCACCACAAACG	BBBBBABBABBBBB@AABBBBABBBBBBBBB:A?BB;B9B?:B>@AA;BA??B7?BBB?B?=BAA*1;:;)''74:*;@9@5?<>A;B?

Now with hexdump,

$ samtools view tmp.bam | head -n 2 | hexdump -C
00000000  45 42 52 49 30 39 33 31  35 31 5f 30 30 35 31 3a  |EBRI093151_0051:|
00000010  35 3a 31 30 3a 31 32 32  33 39 3a 34 38 39 37 23  |5:10:12239:4897#|
00000020  30 09 39 39 09 43 68 72  31 30 09 31 30 09 30 09  |0.99.Chr10.10.0.|
00000030  31 30 31 4d 09 3d 09 34  31 32 09 35 30 33 09 41  |101M.=.412.503.A|
00000040  47 41 41 47 43 43 54 43  41 54 54 54 41 47 41 43  |GAAGCCTCATTTAGAC|
00000050  47 54 47 43 43 54 41 43  41 41 43 41 43 41 47 54  |GTGCCTACAACACAGT|
00000060  47 47 47 54 54 47 47 41  47 54 41 47 41 47 43 41  |GGGTTGGAGTAGAGCA|
00000070  54 43 54 43 43 41 41 41  47 47 43 43 43 54 54 54  |TCTCCAAAGGCCCTTT|
00000080  43 43 41 41 54 43 43 41  41 43 41 54 47 41 47 54  |CCAATCCAACATGAGT|
00000090  41 41 43 54 43 41 41 47  43 54 43 54 47 43 41 43  |AACTCAAGCTCTGCAC|
000000a0  43 41 47 43 09 46 46 46  46 45 46 46 46 46 46 46  |CAGC.FFFFEFFFFFF|
000000b0  46 46 46 46 46 46 46 44  46 46 46 46 46 46 45 46  |FFFFFFFDFFFFFFEF|
000000c0  45 46 46 46 46 43 46 46  46 45 44 44 46 43 45 44  |EFFFFCFFFEDDFCED|
000000d0  46 44 46 44 46 46 45 41  45 46 46 46 40 46 44 44  |FDFDFFEAEFFF@FDD|
000000e0  46 44 46 45 44 3e 44 45  44 44 45 44 45 41 44 44  |FDFED>DEDDEDEADD|
000000f0  44 44 44 44 44 44 44 44  44 3e 44 44 44 41 44 3f  |DDDDDDDDD>DDDAD?|
00000100  44 44 44 44 44 44 44 44  3a 42 0a 45 42 52 49 30  |DDDDDDDD:B.EBRI0|
00000110  39 33 31 35 31 5f 30 30  35 34 3a 36 3a 33 3a 31  |93151_0054:6:3:1|
00000120  37 36 33 3a 31 32 35 35  30 23 30 09 31 36 33 09  |763:12550#0.163.|
00000130  43 68 72 31 30 09 31 34  09 32 39 09 31 30 31 4d  |Chr10.14.29.101M|
00000140  09 3d 09 34 37 30 09 35  35 37 09 47 43 43 54 43  |.=.470.557.GCCTC|
00000150  41 54 54 54 41 47 41 43  47 54 47 43 43 54 41 43  |ATTTAGACGTGCCTAC|
00000160  41 41 43 41 43 41 47 54  47 47 47 54 54 47 47 41  |AACACAGTGGGTTGGA|
00000170  47 54 41 47 41 47 43 41  54 43 54 43 43 41 41 41  |GTAGAGCATCTCCAAA|
00000180  47 47 43 43 43 54 54 54  43 43 41 41 54 43 43 41  |GGCCCTTTCCAATCCA|
00000190  41 43 41 54 47 41 47 41  41 41 43 54 43 41 41 47  |ACATGAGAAACTCAAG|
000001a0  43 54 43 54 47 43 41 43  43 41 43 41 41 41 43 47  |CTCTGCACCACAAACG|
000001b0  09 42 42 42 42 42 41 42  42 41 42 42 42 42 42 40  |.BBBBBABBABBBBB@|
000001c0  41 41 42 42 42 42 41 42  42 42 42 42 42 42 42 42  |AABBBBABBBBBBBBB|
000001d0  3a 41 3f 42 42 3b 42 39  42 3f 3a 42 3e 40 41 41  |:A?BB;B9B?:B>@AA|
000001e0  3b 42 41 3f 3f 42 37 3f  42 42 42 3f 42 3f 3d 42  |;BA??B7?BBB?B?=B|
000001f0  41 41 2a 31 3b 3a 3b 29  27 27 37 34 3a 2a 3b 40  |AA*1;:;)''74:*;@|
00000200  39 40 35 3f 3c 3e 41 3b  42 3f 1e 1e 1e 1e 1e 1e  |9@5?<>A;B?......|
00000210  1e 1e 1e 1e 1e 1e 0a                              |.......|
00000217

You can see that after the "...A;B?" you can see in the SAM output, it actually continues with a dozen non-printable ASCII 0x1e characters (character 30).

Here the dozen quality scores of PHRED 253 would become 253+33 = 286, blindly wrapping at 256 (ignoring the high bits) gives ASCII 30 as observed. I will report this as a bug in samtools.

Do you mind that tmp.bam file being shared and used as a BAM test case? e.g. in samtools, picard, or other BAM handling tools?

It would be worth finding out how that BAM file was created. It could be a bug in the tool used, perhaps invalid input data, or more simply a data corruption (e.g. bad RAM, bad disk, failed network transfer).

[tomjan]

You're the man!
Thank you a lot.

if -3 =253 i can handle it in perlscript

>Do you mind that tmp.bam file being shared and used as a BAM test case? e.g.
> in samtools, picard, or other BAM handling tools?

because of my english I dont understand it. are you asking my if you may share tmp.bam?
if yes, the answer is ' yes of course'. (I have just shared it by posting to forum)

>It would be worth finding out how that BAM file was created.

I will ask the 'bam provider' about sequencing hardware and software.

tj

[maubp]

You're the man!
Thank you a lot.

if -3 =253 i can handle it in perlscript

In this specific case, yes.

>Do you mind that tmp.bam file being shared and used as a BAM test case? e.g.
> in samtools, picard, or other BAM handling tools?

because of my english I dont understand it. are you asking my if you may share tmp.bam?
if yes, the answer is ' yes of course'. (I have just shared it by posting to forum)

Thank you.

[maubp]

I will report this as a bug in samtools.

See this thread, http://sourceforge.net/mailarchive/m...sg_id=29178587