From: Petr Danecek <pd3@sanger.ac.uk>
Date: Fri, 25 May 2012 07:41:01 +0000 (+0100)
Subject: Fix to prevent long stretches of Ns to be mistaken for indels; Always output the... 
X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=8c15f916dabce475febdf508a9cc0c708c5a7747

Fix to prevent long stretches of Ns to be mistaken for indels; Always output the most likely ALT allele
---

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 5142b3e..11cd371 100644
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -161,6 +161,11 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
 				if (j > max_rd_len) max_rd_len = j;
 			}
 		}
+        // To prevent long stretches of N's to be mistaken for indels (sometimes thousands of bases),
+        //  check the number of N's in the sequence. TODO: this may not be the best place and the best way of doing it
+        int nN=0; for (i=0; i<max_rd_len && ref[i]; i++) if ( ref[i]=='N' ) nN++;
+        if ( nN*2>i ) return -1;
+
 		ks_introsort(uint32_t, m, aux);
 		// squeeze out identical types
 		for (i = 1, n_types = 1; i < m; ++i)
diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index a380484..15735af 100644
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -193,24 +193,21 @@ void bcf_p1_destroy(bcf_p1aux_t *ma)
 
 static int cal_pdg(const bcf1_t *b, bcf_p1aux_t *ma)
 {
-	int i, j;
-	long *p, tmp;
-	p = alloca(b->n_alleles * sizeof(long));
-	memset(p, 0, sizeof(long) * b->n_alleles);
+	int i, j, imax=0;
 	for (j = 0; j < ma->n; ++j) {
 		const uint8_t *pi = ma->PL + j * ma->PL_len;
 		double *pdg = ma->pdg + j * 3;
 		pdg[0] = ma->q2p[pi[2]]; pdg[1] = ma->q2p[pi[1]]; pdg[2] = ma->q2p[pi[0]];
-		for (i = 0; i < b->n_alleles; ++i)
-			p[i] += (int)pi[(i+1)*(i+2)/2-1];
+        int ib,ia=0,n=(b->n_alleles+1)*b->n_alleles/2;
+        for (i=0; i<n;) {
+            for (ib=0; ib<=ia; ib++) {
+                if ( pi[i]==0 && ia>imax ) imax=ia;
+                i++;
+            }
+            ia++;
+        }
 	}
-	for (i = 0; i < b->n_alleles; ++i) p[i] = p[i]<<4 | i;
-	for (i = 1; i < b->n_alleles; ++i) // insertion sort
-		for (j = i; j > 0 && p[j] < p[j-1]; --j)
-			tmp = p[j], p[j] = p[j-1], p[j-1] = tmp;
-	for (i = b->n_alleles - 1; i >= 0; --i)
-		if ((p[i]&0xf) == 0) break;
-	return i;
+	return imax;
 }
 
 int bcf_p1_call_gt(const bcf_p1aux_t *ma, double f0, int k)