]> git.donarmstrong.com Git - samtools.git/commitdiff
Update BIN values in 'samtools depad'
authorpeterjc <p.j.a.cock@googlemail.com>
Fri, 19 Apr 2013 11:08:42 +0000 (12:08 +0100)
committerpeterjc <p.j.a.cock@googlemail.com>
Fri, 19 Apr 2013 11:08:46 +0000 (12:08 +0100)
By its nature, 'samtools depad' edits mapped reads, altering both
the mapping position and CIGAR string. This means the mapped start
and end can change, and therefore the BIN can change.

Prior to this fix the pre-editing BIN (from the padded reference)
was being written to the BAM output, leading to subtle problems
downstream when trying to use the depadded BAM file via the BAI
index. For example, in Picard v1.89 this could trigger an
java.lang.ArrayIndexOutOfBoundsException where an old BIN was
outside the length of the unpadded reference sequence.

padding.c

index a8da562d0872e8b0930bb2bdd3dcc9abb58bab7b..476e393166983b0d00434789ad8596f0403a8f34 100644 (file)
--- a/padding.c
+++ b/padding.c
@@ -195,6 +195,7 @@ int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai)
                        write_cigar(cigar2, n2, m2, bam_cigar_gen(b->core.l_qseq, BAM_CMATCH));
                        replace_cigar(b, n2, cigar2);
                        posmap = update_posmap(posmap, r);
+                       b->core.bin = bam_reg2bin(0, bam_calend(&b->core, bam1_cigar(b)));
                } else if (b->core.n_cigar > 0) {
                        int i, k, op;
                        if (b->core.tid < 0) {
@@ -272,6 +273,7 @@ int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai)
                        n2 = k;
                        replace_cigar(b, n2, cigar2);
                        b->core.pos = posmap[b->core.pos];
+                       b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&b->core.pos, bam1_cigar(b)));
                        if (b->core.mtid < 0 || b->core.mpos < 0) {
                                /* Nice case, no mate to worry about*/
                                // fprintf(stderr, "[depad] Read '%s' mate not mapped\n", bam1_qname(b));