Fix CRAM embed_ref=2 with seqs overlapping ref end.

cram/cram_encode.c

@@ -1761,7 +1761,6 @@ static int cram_generate_reference(cram_container *c, cram_slice *s, int r1) {
     c->ref_start = ref_start+1;
     c->ref_end   = ref_end+1;
     c->ref_free  = 1;
-
     return 0;
 
  err:
@@ -1997,6 +1996,9 @@ int cram_encode_container(cram_fd *fd, cram_container *c) {
                 fd->no_ref_counter -= (fd->no_ref_counter > 0);
                 pthread_mutex_unlock(&fd->ref_lock);
             }
+
+            if (c->ref_end > fd->refs->ref_id[c->ref_id]->LN_length)
+                c->ref_end = fd->refs->ref_id[c->ref_id]->LN_length;
         }
 
         // Iterate through records creating the cram blocks for some

cram/cram_io.c

@@ -2787,10 +2787,12 @@ static int refs_from_header(cram_fd *fd) {
 
         /* Initialise likely filename if known */
         if ((ty = sam_hrecs_find_type_id(h->hrecs, "SQ", "SN", h->hrecs->ref[i].name))) {
-            if ((tag = sam_hrecs_find_key(ty, "M5", NULL))) {
+            if ((tag = sam_hrecs_find_key(ty, "M5", NULL)))
                 r->ref_id[j]->fn = string_dup(r->pool, tag->str+3);
-                //fprintf(stderr, "Tagging @SQ %s / %s\n", r->ref_id[h]->name, r->ref_id[h]->fn);
-            }
+
+            if ((tag = sam_hrecs_find_key(ty, "LN", NULL)))
+                // LN tag used when constructing consensus reference
+                r->ref_id[j]->LN_length = strtoll(tag->str+3, NULL, 0);
         }
 
         k = kh_put(refs, r->h_meta, r->ref_id[j]->name, &n);

cram/cram_structs.h

@@ -671,7 +671,8 @@ struct cram_slice {
 typedef struct ref_entry {
     char *name;
     char *fn;
-    int64_t length;
+    int64_t length;        // if 0 this indicates we haven't loaded it yet
+    int64_t LN_length;     // @SQ LN length, used to trim consensus ref
     int64_t offset;
     int bases_per_line;
     int line_length;

test/sam.c

@@ -2254,7 +2254,7 @@ static void test_bam_set1_write_and_read_back(void)
     w_bam = bam_init1();
     VERIFY(w_bam != NULL, "failed to initialize BAM struct.");
     r = bam_set1(w_bam, strlen(qname), qname,
-                 BAM_FREVERSE, 0, 1000, 42,
+                 BAM_FPAIRED | BAM_FREVERSE, 0, 1000, 42,
                  sizeof(cigar) / 4, cigar, 0, 2000, 3000,
                  strlen(seq), seq, qual, 64);
     VERIFY(r >= 0, "call to bam_set1() failed.");

test/test.pl

@@ -865,6 +865,14 @@ sub test_view
     testv $opts, "./test_view $tv_args -C -p $ercram $ersam";
     testv $opts, "./test_view $tv_args -p $ersam2 $ercram";
     testv $opts, "./compare_sam.pl $ersam $ersam2";
+
+    $ersam = "c1#bounds.sam";
+    $ercram = "c1#bounds_er.tmp.cram";
+    $ersam2 = "${ercram}.sam";
+    testv $opts, "./test_view $tv_args -C -p $ercram $ersam";
+    testv $opts, "./test_view $tv_args -p $ersam2 $ercram";
+    testv $opts, "./compare_sam.pl $ersam $ersam2";
+
     if ($test_view_failures == 0) {
         passed($opts, "embed_ref=2 tests");
     } else {