Revert one of my earlier changes - Heng was right, CIGAR P not sensible in a padded...

[samtools.git] / bcftools / prob1.c

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 667a2732555dae4df733f12859c646eea6d63488..83bd8e2881f2f84688ce287b9a2ebe5e785cdace 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -194,23 +194,49 @@ 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 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]];
-               for (i = 0; i < b->n_alleles; ++i)
-                       p[i] += (int)pi[(i+1)*(i+2)/2-1];
-       }
-       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;
+        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)
@@ -488,6 +514,7 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r
                        break;
                }
        }
+       if (i == b->n_gi) return -1; // no PL
        if (b->n_alleles < 2) return -1; // FIXME: find a better solution
        // 
        rst->rank0 = cal_pdg(b, ma);
@@ -496,6 +523,13 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r
        for (k = 0, sum = 0.; k < ma->M; ++k)
                sum += ma->afs1[k];
        rst->p_var = (double)sum;
+       { // compute the allele count
+               double max = -1;
+               rst->ac = -1;
+               for (k = 0; k <= ma->M; ++k)
+                       if (max < ma->z[k]) max = ma->z[k], rst->ac = k;
+               rst->ac = ma->M - rst->ac;
+       }
        // calculate f_flat and f_em
        for (k = 0, sum = 0.; k <= ma->M; ++k)
                sum += (long double)ma->z[k];
@@ -508,16 +542,27 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r
        { // estimate equal-tail credible interval (95% level)
                int l, h;
                double p;
-               for (i = 0, p = 0.; i < ma->M; ++i)
+               for (i = 0, p = 0.; i <= ma->M; ++i)
                        if (p + ma->afs1[i] > 0.025) break;
                        else p += ma->afs1[i];
                l = i;
-               for (i = ma->M-1, p = 0.; i >= 0; --i)
+               for (i = ma->M, p = 0.; i >= 0; --i)
                        if (p + ma->afs1[i] > 0.025) break;
                        else p += ma->afs1[i];
                h = i;
                rst->cil = (double)(ma->M - h) / ma->M; rst->cih = (double)(ma->M - l) / ma->M;
        }
+       if (ma->n1 > 0) { // compute LRT
+               double max0, max1, max2;
+               for (k = 0, max0 = -1; k <= ma->M; ++k)
+                       if (max0 < ma->z[k]) max0 = ma->z[k];
+               for (k = 0, max1 = -1; k <= ma->n1 * 2; ++k)
+                       if (max1 < ma->z1[k]) max1 = ma->z1[k];
+               for (k = 0, max2 = -1; k <= ma->M - ma->n1 * 2; ++k)
+                       if (max2 < ma->z2[k]) max2 = ma->z2[k];
+               rst->lrt = log(max1 * max2 / max0);
+               rst->lrt = rst->lrt < 0? 1 : kf_gammaq(.5, rst->lrt);
+       } else rst->lrt = -1.0;
        rst->cmp[0] = rst->cmp[1] = rst->cmp[2] = rst->p_chi2 = -1.0;
        if (do_contrast && rst->p_var > 0.5) // skip contrast2() if the locus is a strong non-variant
                rst->p_chi2 = contrast2(ma, rst->cmp);