From eb4c2ebca7bd29d3b850baf9b57c342742d06788 Mon Sep 17 00:00:00 2001
From: Petr Danecek <pd3@sanger.ac.uk>
Date: Wed, 30 May 2012 10:20:51 +0100
Subject: [PATCH] Fixed errors introduced by the
 8c15f916dabce475febdf508a9cc0c708c5a7747 commit.

---
 bam2bcf_indel.c  |  4 ++--
 bcftools/prob1.c | 51 +++++++++++++++++++++++++++++++++++++-----------
 2 files changed, 42 insertions(+), 13 deletions(-)

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 11cd371..ab9e83c 100644
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -162,8 +162,8 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
 			}
 		}
         // 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++;
+        //  check the number of N's in the sequence and skip places where half or more reference bases are Ns.
+        int nN=0; for (i=pos; i-pos<max_rd_len && ref[i]; i++) if ( ref[i]=='N' ) nN++;
         if ( nN*2>i ) return -1;
 
 		ks_introsort(uint32_t, m, aux);
diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 15735af..83bd8e2 100644
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -193,21 +193,50 @@ void bcf_p1_destroy(bcf_p1aux_t *ma)
 
 static int cal_pdg(const bcf1_t *b, bcf_p1aux_t *ma)
 {
-	int i, j, imax=0;
+	int i, j;
+    int n = (b->n_alleles+1)*b->n_alleles/2;
+    double *lk = alloca(n * sizeof(long));
+    memset(lk, 0, sizeof(double) * n);
 	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]];
-        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++;
-        }
-	}
-	return imax;
+        for (i=0; i<n; i++) lk[i] += pi[i];
+    }
+
+    double norm=lk[0]; 
+    for (i=1; i<n; i++) if (lk[i]<norm) norm=lk[i];
+    #if DBG
+    for (i=0,j=0; i<b->n_alleles; i++)
+    {
+        int k; for (k=0; k<=i; k++) printf("%.0f\t", lk[j++]);
+        printf("\n");
+    }
+    #endif
+    for (i=0; i<n; i++) lk[i] = pow(10,-0.1*(lk[i]-norm));
+
+    // Find out the most likely alleles. In contrast to the original version,
+    //  ALT alleles may not be printed when they are more likely than REF but
+    //  significantly less likely than the most likely ALT. The only criterion
+    //  is the LK ratio now. To obtain behaviour similar to the original one,
+    //  use the pref variant below.
+    double pmax=0; //,pref=0;
+    n = ma->is_indel ? b->n_alleles : b->n_alleles-1;
+    for (i=0; i<n; i++)
+    {
+        double pr=0;
+        int k=i*(i+1)/2;
+        for (j=0; j<=i; j++) { pr+=lk[k]; k++; }
+        for (j=i+1; j<b->n_alleles; j++) { k=j*(j+1)/2+i; pr+=lk[k]; }
+        #if DBG
+        printf("%d\t%e\n", i,pr);
+        #endif
+        if (pmax<pr) pmax=pr;
+        // if (i==0) pref=pr;
+        // if (pr<pref && pr/pmax < 1e-4) break;
+        if (pr/pmax < 1e-4) break;   // Assuming the alleles are sorted by the lk
+    }
+	return i-1;
 }
 
 int bcf_p1_call_gt(const bcf_p1aux_t *ma, double f0, int k)
-- 
2.39.2