From: Heng Li <lh3@live.co.uk>
Date: Tue, 26 Oct 2010 19:27:46 +0000 (+0000)
Subject: remove local realignment (probabilistic realignment is still there)
X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=8a994c1783c189bcf23777ca911c4b0eeeb66d1e;p=samtools.git

remove local realignment (probabilistic realignment is still there)
---

diff --git a/bam_md.c b/bam_md.c
index 7023046..f65b845 100644
--- 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;