Fix CRAM QS values where ambiguity codes are in the reference.

Commit 0aa2af3 changes QS encoding from verbatim copying out of
s->qual_blk and into calling codecs[DS_QS]->encode.  This was part of
the plan to be able to perform the data transforms on QS, included in
CRAM 4.0 prototype.

Unfortunately the ordering was wrong, with a copy of all quals
followed by encoding any explicit QS items produced by the feature
codes.  The correct ordering is the reversal of this.

Author: jkbonfield

io_lib/cram_encode.c

@@ -535,6 +535,7 @@ static int cram_encode_slice_read(cram_fd *fd,
     int32_t i32;
     int64_t i64;
     unsigned char uc;
+    int explicit_qual = 0;
 
     //fprintf(stderr, "Encode seq %d, %d/%d FN=%d, %s\n", rec, core->byte, core->bit, cr->nfeature, s->name_ds->str + cr->name);
 
@@ -609,11 +610,6 @@ static int cram_encode_slice_read(cram_fd *fd,
     /* Aux tags */
     r |= h->codecs[DS_TL]->encode(s, h->codecs[DS_TL], (char *)&cr->TL, 1);
 
-    // qual
-    r |= h->codecs[DS_QS]->encode(s, h->codecs[DS_QS],
-				  (char *)BLOCK_DATA(s->qual_blk) + cr->qual,
-				  cr->len);
-
     // features (diffs)
     if (!(cr->flags & BAM_FUNMAP)) {
 	int prev_pos = 0, j;
@@ -686,6 +682,7 @@ static int cram_encode_slice_read(cram_fd *fd,
 		uc  = f->B.qual;
 		r |= h->codecs[DS_QS]->encode(s, h->codecs[DS_QS],
 					      (char *)&uc, 1);
+		explicit_qual++;
 		break;
 
 	    case 'b':
@@ -700,6 +697,7 @@ static int cram_encode_slice_read(cram_fd *fd,
 		uc  = f->Q.qual;
 		r |= h->codecs[DS_QS]->encode(s, h->codecs[DS_QS],
 					      (char *)&uc, 1);
+		explicit_qual++;
 		break;
 
 	    case 'N':
@@ -736,6 +734,11 @@ static int cram_encode_slice_read(cram_fd *fd,
 	    r |= h->codecs[DS_BA]->encode(s, h->codecs[DS_BA], seq, cr->len);
     }
 
+    // qual
+    r |= h->codecs[DS_QS]->encode(s, h->codecs[DS_QS],
+				  (char *)BLOCK_DATA(s->qual_blk) + cr->qual
+				  + explicit_qual, cr->len);
+
     return r ? -1 : 0;
 }
 

tests/data/xx#MD.sam

@@ -1,7 +1,7 @@
 @SQ	SN:xx	LN:30
 @CO	All MD and NM should match the stored values
 a	0	xx	6	1	10M	*	0	0	AAAAATTTTT	*	co:Z:no fields
-a	0	xx	6	1	10M	*	0	0	AAAAGGTTTT	*
+a	0	xx	6	1	11M	*	0	0	AAAAGRTTTTT	ABCDEFGHIJK
 a	0	xx	6	1	10M	*	0	0	GAAAATTTTG	*
 i	0	xx	6	1	5M1I5M	*	0	0	AAAAAGTTTTT	*
 i	0	xx	6	1	5M3I5M	*	0	0	AAAAAGGGTTTTT	*
@@ -11,12 +11,12 @@ d	0	xx	6	1	5M10D5M	*	0	0	AAAAACCCCC	*
 d	0	xx	6	1	5M10N5M	*	0	0	AAAAACCCCC	*
 sid	0	xx	6	1	1S4M10D5I4M1S	*	0	0	AAAAAGGGGGCCCCC	*
 a	0	xx	6	1	10M	*	0	0	AAAAATTTTT	*	MD:Z:10	NM:i:0	co:Z:correct fields
-a	0	xx	6	1	10M	*	0	0	AAAAGGTTTT	*	MD:Z:4A0T4	NM:i:2
+a	0	xx	6	1	11M	*	0	0	AAAAGRTTTTT	ABCDEFGHIJK	MD:Z:4A0T4Y0	NM:i:3
 a	0	xx	6	1	10M	*	0	0	GAAAATTTTG	*	MD:Z:0A8T0	NM:i:2
 i	0	xx	6	1	5M1I5M	*	0	0	AAAAAGTTTTT	*	MD:Z:10	NM:i:1
 i	0	xx	6	1	5M3I5M	*	0	0	AAAAAGGGTTTTT	*	MD:Z:10	NM:i:3
 i	0	xx	6	1	10M2I	*	0	0	AAAAATTTTTCC	*	MD:Z:10	NM:i:2
 i	0	xx	6	1	10M2P2I	*	0	0	AAAAATTTTTCC	*	MD:Z:10	NM:i:2
-d	0	xx	6	1	5M10D5M	*	0	0	AAAAACCCCC	*	MD:Z:5^TTTTTTTTTT5	NM:i:10
+d	0	xx	6	1	5M10D5M	*	0	0	AAAAACCCCC	*	MD:Z:5^TTTTTYTTTT5	NM:i:10
 d	0	xx	6	1	5M10N5M	*	0	0	AAAAACCCCC	*	MD:Z:10	NM:i:0
-sid	0	xx	6	1	1S4M10D5I4M1S	*	0	0	AAAAAGGGGGCCCCC	*	MD:Z:4^ATTTTTTTTT0T3	NM:i:16
+sid	0	xx	6	1	1S4M10D5I4M1S	*	0	0	AAAAAGGGGGCCCCC	*	MD:Z:4^ATTTTTYTTT0T3	NM:i:16

tests/data/xx#MD2.sam

@@ -1,7 +1,8 @@
 @SQ	SN:xx	LN:30
 @CO	All MD and/or NM should differ to the stored values
 a	0	xx	6	1	10M	*	0	0	AAAAATTTTT	*	MD:Z:9	NM:i:0	co:Z:MD incorrect fields
-a	0	xx	6	1	10M	*	0	0	AAAAGGTTTT	*	MD:Z:4A0A4	NM:i:2
+a	0	xx	6	1	11M	*	0	0	AAAAGGTTTTT	*	MD:Z:4A0T4Y0	NM:i:2
+a	0	xx	6	1	11M	*	0	0	AAAAGGTTTTT	*	MD:Z:4A0T4N0	NM:i:3
 a	0	xx	6	1	10M	*	0	0	GAAAATTTTG	*	MD:Z:0G8T0	NM:i:2
 i	0	xx	6	1	5M1I5M	*	0	0	AAAAAGTTTTT	*	MD:Z:11	NM:i:1
 i	0	xx	6	1	5M3I5M	*	0	0	AAAAAGGGTTTTT	*	MD:Z:1A1	NM:i:3

tests/data/xx#rg.sam

@@ -1,5 +1,5 @@
 @HD	VN:1.4	SO:coordinate
-@SQ	SN:xx	LN:30	AS:?	SP:?	UR:?	M5:bbf4de6d8497a119dda6e074521643dc
+@SQ	SN:xx	LN:30	AS:?	SP:?	UR:?	M5:1224b81d8664d77635e1620d7f2c1523
 @RG	ID:x1	SM:x1
 @RG	ID:x2	SM:x2	LB:x	PG:foo:bar	PI:1111
 @PG	ID:emacs	PN:emacs	VN:23.1.1

tests/data/xx.fa

@@ -1,5 +1,5 @@
 >xx
-AAAAAAAAAATTTTTTTTTTCCCCCCCCCC
+AAAAAAAAAATTTTTYTTTTCCCCCCCCCC
 >yy
 AAAAAAAAAATTTTTTTTTT