remove local realignment (probabilistic realignment is still there)

diff --git a/bam_md.c b/bam_md.c
index 7023046b87748de170de994a2565ec0b7a8172a2..f65b845b5d594980d4b4dba36122672c71b19cc9 100644 (file)
--- a/bam_md.c
+++ b/bam_md.c
@@ -161,106 +161,6 @@ int bam_cap_mapQ(bam1_t *b, char *ref, int thres)
        return (int)(t + .499);
 }
 
-// local realignment
-
-#define MIN_REF_LEN 10
-#define MIN_BAND_WIDTH 11
-
-int bam_realn(bam1_t *b, const char *ref)
-{
-       int k, score, q[2], r[2], kk[2], kl[2], x, y, max, j, n_cigar, endx = -1;
-       uint32_t *cigar = bam1_cigar(b);
-       uint8_t *seq = bam1_seq(b);
-       bam1_core_t *c = &b->core;
-       q[0] = q[1] = r[0] = r[1] = kk[0] = kk[1] = kl[0] = kl[1] = -1;
-       // find the right boundary
-       for (k = 0, score = max = 0, x = c->pos, y = 0; k < c->n_cigar; ++k) {
-               int op = cigar[k]&0xf;
-               int ol = cigar[k]>>4;
-               if (op == BAM_CMATCH) {
-                       for (j = 0; j < ol; ++j) {
-                               int c1, c2, z = y + j;
-                               c1 = bam_nt16_nt4_table[bam1_seqi(seq, z)];
-                               if (ref[x+j] == 0) return -1;
-                               c2 = bam_nt16_nt4_table[(int)bam_nt16_table[(int)ref[x+j]]];
-                               if (c1 < 3 && c2 < 3)
-                                       score += c1 == c2? 5 : -4;
-                               if (score < 0) score = 0;
-                               if (score > max) max = score, q[1] = z, r[1] = x+j, kk[1] = k, kl[1] = j + 1;
-                       }
-                       x += ol; y += ol;
-               } else if (op == BAM_CINS) {
-                       score -= 4 - ol * 3;
-                       y += ol;
-                       if (score < 0) score = 0;
-               } else if (op == BAM_CDEL) {
-                       score -= 4 - ol * 3;
-                       x += ol;
-                       if (score < 0) score = 0;
-               } else if (op == BAM_CSOFT_CLIP) y += ol;
-               else if (op == BAM_CREF_SKIP) x += ol;
-       }
-       if (q[1] < 0) return -1; // no high scoring segments
-       endx = x - 1;
-       // find the left boundary
-       for (k = c->n_cigar - 1, score = max = 0, x = x-1, y = y-1; k >= 0; --k) {
-               int op = cigar[k]&0xf;
-               int ol = cigar[k]>>4;
-               if (op == BAM_CMATCH) {
-                       for (j = 0; j < ol; ++j) {
-                               int c1, c2, z = y - j;
-                               c1 = bam_nt16_nt4_table[bam1_seqi(seq, z)];
-                               if (ref[x+j] == 0) return -1;
-                               c2 = bam_nt16_nt4_table[(int)bam_nt16_table[(int)ref[x-j]]];
-                               if (c1 < 3 && c2 < 3)
-                                       score += c1 == c2? 5 : -4;
-                               if (score < 0) score = 0;
-                               if (score > max) max = score, q[0] = z, r[0] = x-j, kk[0] = k, kl[0] = j + 1;
-                       }
-                       x -= ol; y -= ol;
-               } else if (op == BAM_CINS) {
-                       score -= 4 - ol * 3;
-                       y -= ol;
-                       if (score < 0) score = 0;
-               } else if (op == BAM_CDEL) {
-                       score -= 4 - ol * 3;
-                       x -= ol;
-                       if (score < 0) score = 0;
-               } else if (op == BAM_CSOFT_CLIP) y -= ol;
-               else if (op == BAM_CREF_SKIP) x -= ol;
-       }
-       if (q[0] < 0 || q[1] - q[0] < 15) return -1; // the high-scoring segment is too short
-       // modify CIGAR
-       n_cigar = 0;
-       cigar = calloc(c->n_cigar + 4, 4);
-       if (q[0] != 0) cigar[n_cigar++] = (uint32_t)q[0]<<4 | BAM_CSOFT_CLIP;
-       if (r[0] != c->pos) cigar[n_cigar++] = (uint32_t)(r[0] - c->pos)<<4 | BAM_CREF_SKIP;
-       if (kk[0] == kk[1]) {
-               cigar[n_cigar++] = (uint32_t)(kl[0] + kl[1] - (bam1_cigar(b)[kk[0]]>>4))<<4 | BAM_CMATCH;
-       } else {
-               cigar[n_cigar++] = (uint32_t)kl[0]<<4 | BAM_CMATCH;
-               for (k = kk[0] + 1; k < kk[1]; ++k)
-                       cigar[n_cigar++] = bam1_cigar(b)[k];
-               cigar[n_cigar++] = (uint32_t)kl[1]<<4 | BAM_CMATCH; // FIXME: add ref_skip after this line
-       }
-       if (q[1] + 1 != c->l_qseq)
-               cigar[n_cigar++] = (uint32_t)(c->l_qseq - q[1] - 1)<<4 | BAM_CSOFT_CLIP;
-       // copy over the alignment
-       if (4 * (n_cigar - (int)c->n_cigar) + b->data_len > b->m_data) { // enlarge b->data
-               b->m_data = 4 * (n_cigar - (int)c->n_cigar) + b->data_len;
-               kroundup32(b->m_data);
-               b->data = realloc(b->data, b->m_data);
-       }
-       if (n_cigar != (int)c->n_cigar) { // move data
-               memmove(b->data + c->l_qname + 4 * n_cigar, bam1_seq(b), b->data_len - c->l_qname - 4 * c->n_cigar);
-               b->data_len += 4 * (n_cigar - (int)c->n_cigar);
-       }
-       memcpy(bam1_cigar(b), cigar, n_cigar * 4);
-       c->n_cigar = n_cigar;
-       free(cigar);
-       return 0;
-}
-
 int bam_prob_realn(bam1_t *b, const char *ref)
 {
        int k, i, bw, x, y, yb, ye, xb, xe;