Include leading CIGAR P operators in 'samtools depad'

[samtools.git] / padding.c

diff --git a/padding.c b/padding.c
index 6c8b31ff3de5d01afa79bb94428bd4e711078116..e2afefb4840cd8acca35f07f200501b87990051a 100644 (file)
--- a/padding.c
+++ b/padding.c
@@ -38,8 +38,7 @@ static void unpad_seq(bam1_t *b, kstring_t *s)
                op = bam_cigar_op(cigar[k]);
                ol = bam_cigar_oplen(cigar[k]);
                if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) {
-                       for (i = 0; i < ol; ++i) s->s[s->l++] = bam1_seqi(seq, j);
-                       ++j;
+                       for (i = 0; i < ol; ++i, ++j) s->s[s->l++] = bam1_seqi(seq, j);
                } else if (op == BAM_CSOFT_CLIP) {
                        j += ol;
                } else if (op == BAM_CDEL || op == BAM_CPAD) {
@@ -93,6 +92,9 @@ int bam_pad2unpad(bamFile in, bamFile out)
                        }
                        unpad_seq(b, &q);
                        if (bam_cigar_op(cigar[0]) == BAM_CSOFT_CLIP) write_cigar(cigar2, n2, m2, cigar[0]);
+                       /* Include any pads if starts with an insert */
+                       for (k = 0; k+1 < b->core.pos && !r.s[b->core.pos - k - 1]; ++k);
+                       if (k) write_cigar(cigar2, n2, m2, bam_cigar_gen(k, BAM_CPAD));
                        /* Determine CIGAR operator for each base in the aligned read */
                        for (i = 0, k = b->core.pos; i < q.l; ++i, ++k)
                                q.s[i] = q.s[i]? (r.s[k]? BAM_CMATCH : BAM_CINS) : (r.s[k]? BAM_CDEL : BAM_CPAD);