From: Heng Li <lh3@live.co.uk>
Date: Thu, 27 Jan 2011 19:50:02 +0000 (+0000)
Subject:  * fixed a rare memory issue in bam_md.c
X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=d1643d6e5e6ec3a7cb1ef8b62b9a6f9b21ff2527;p=samtools.git

 * fixed a rare memory issue in bam_md.c
 * fixed a bug in indel calling related to unmapped and refskip reads
---

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index c71280c..c157bc6 100644
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -318,8 +318,14 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
 			// align each read to ref2
 			for (i = 0; i < n_plp[s]; ++i, ++K) {
 				bam_pileup1_t *p = plp[s] + i;
-				int qbeg, qend, tbeg, tend, sc;
+				int qbeg, qend, tbeg, tend, sc, kk;
 				uint8_t *seq = bam1_seq(p->b);
+				uint32_t *cigar = bam1_cigar(p->b);
+				if (p->b->core.flag&4) continue; // unmapped reads
+				// FIXME: the following loop should be better moved outside; nonetheless, realignment should be much slower anyway.
+				for (kk = 0; kk < p->b->core.n_cigar; ++kk)
+					if ((cigar[kk]&BAM_CIGAR_MASK) == BAM_CREF_SKIP) break;
+				if (kk < p->b->core.n_cigar) continue;
 				// FIXME: the following skips soft clips, but using them may be more sensitive.
 				// determine the start and end of sequences for alignment
 				qbeg = tpos2qpos(&p->b->core, bam1_cigar(p->b), left,  0, &tbeg);
@@ -414,7 +420,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
 				indelQ2 = tmp > 111? 0 : (int)((1. - tmp/111.) * indelQ2 + .499);
 				// pick the smaller between indelQ1 and indelQ2
 				indelQ = indelQ1 < indelQ2? indelQ1 : indelQ2;
-				p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ;
+				p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ; // use 22 bits in total
 				sumq[sc[0]&0x3f] += indelQ < seqQ? indelQ : seqQ;
 //				fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), types[sc[0]&0x3f], indelQ);
 			}
diff --git a/bam_md.c b/bam_md.c
index a1a6d66..02d6c5b 100644
--- a/bam_md.c
+++ b/bam_md.c
@@ -223,8 +223,10 @@ int bam_prob_realn_core(bam1_t *b, const char *ref, int flag)
 		s = calloc(c->l_qseq, 1);
 		for (i = 0; i < c->l_qseq; ++i) s[i] = bam_nt16_nt4_table[bam1_seqi(seq, i)];
 		r = calloc(xe - xb, 1);
-		for (i = xb; i < xe; ++i)
+		for (i = xb; i < xe; ++i) {
+			if (ref[i] == 0) { xe = i; break; }
 			r[i-xb] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[i]]];
+		}
 		state = calloc(c->l_qseq, sizeof(int));
 		q = calloc(c->l_qseq, 1);
 		kpa_glocal(r, xe-xb, s, c->l_qseq, qual, &conf, state, q);