VCF annotations with and without -m and more robust rm_info

diff --git a/bcftools/bcf.c b/bcftools/bcf.c
index 87ff6c7435abc65758017bf3f6307df803e0606c..6695c74226059ea55033792dabab3faba4dbda05 100644 (file)
--- a/bcftools/bcf.c
+++ b/bcftools/bcf.c
@@ -323,7 +323,7 @@ int bcf_append_info(bcf1_t *b, const char *info, int l)
        return 0;
 }
 
-int remove_tag(char *str, char *tag, char delim)
+int remove_tag(char *str, const char *tag, char delim)
 {
     char *tmp = str, *p;
     int len_diff = 0, ori_len = strlen(str);
@@ -336,11 +336,11 @@ int remove_tag(char *str, char *tag, char delim)
         }
         char *q=p+1;
         while ( *q && *q!=delim ) q++;
-        if ( p==str && *q ) q++;
+        if ( p==str && *q ) q++;        // the tag is first, don't move the delim char
         len_diff += q-p;
-        if ( ! *q ) { *p = 0; break; }
+        if ( ! *q ) { *p = 0; break; }  // the tag was last, no delim follows
         else
-            memmove(p,q,ori_len-(int)(q-p)+1);
+            memmove(p,q,ori_len-(int)(p-str)-(int)(q-p)-1);  // *q==delim
     }
     if ( len_diff==ori_len )
         str[0]='.', str[1]=0, len_diff--;
@@ -348,6 +348,25 @@ int remove_tag(char *str, char *tag, char delim)
     return len_diff;
 }
 
+
+void rm_info(kstring_t *s, const char *key)
+{
+    char *p = s->s; 
+    int n = 0;
+    while ( n<4 )
+    {
+        if ( !*p ) n++;
+        p++;
+    }
+    char *q = p+1; 
+    while ( *q && q-s->s<s->l ) q++;
+
+    int nrm = remove_tag(p, key, ';');
+    if ( nrm )
+        memmove(q-nrm, q, s->s+s->l-q+1);
+    s->l -= nrm;
+}
+
 int bcf_cpy(bcf1_t *r, const bcf1_t *b)
 {
        char *t1 = r->str;

diff --git a/bcftools/bcf.h b/bcftools/bcf.h
index 7e451a67f5dad91f95c287b6c08b7cd7d278cf40..8c52451f7a40d939a322a2737f922ac44fa9e871 100644 (file)
--- a/bcftools/bcf.h
+++ b/bcftools/bcf.h
@@ -124,7 +124,9 @@ extern "C" {
        // append more info
        int bcf_append_info(bcf1_t *b, const char *info, int l);
     // remove tag
-    int remove_tag(char *string, char *tag, char delim);
+    int remove_tag(char *string, const char *tag, char delim);
+    // remove info tag, string is the kstring holder of bcf1_t.str
+    void rm_info(kstring_t *string, const char *key);
        // copy
        int bcf_cpy(bcf1_t *r, const bcf1_t *b);
 

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 415426abcd1511f1e2cadd987501efcbf21bcb74..22ff2aca91be76f564e7bd3825abddcd4616ba42 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -48,11 +48,6 @@ void *bed_read(const char *fn);
 void bed_destroy(void *_h);
 int bed_overlap(const void *_h, const char *chr, int beg, int end);
 
-typedef struct {
-       double p[4];
-       int mq, depth, is_tested, d[4];
-} anno16_t;
-
 static double ttest(int n1, int n2, int a[4])
 {
        extern double kf_betai(double a, double b, double x);
@@ -83,7 +78,7 @@ static int test16_core(int anno[16], anno16_t *a)
        return 0;
 }
 
-static int test16(bcf1_t *b, anno16_t *a)
+int test16(bcf1_t *b, anno16_t *a)
 {
        char *p;
        int i, anno[16];
@@ -100,17 +95,6 @@ static int test16(bcf1_t *b, anno16_t *a)
        return test16_core(anno, a);
 }
 
-static void rm_info(bcf1_t *b, const char *key)
-{
-       char *p, *q;
-       if ((p = strstr(b->info, key)) == 0) return;
-       for (q = p; *q && *q != ';'; ++q);
-       if (p > b->info && *(p-1) == ';') --p;
-       memmove(p, q, b->l_str - (q - b->str));
-       b->l_str -= q - p;
-       bcf_sync(b);
-}
-
 static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10], int cons_llr, int64_t cons_gt)
 {
        kstring_t s;
@@ -119,7 +103,7 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr,
        anno16_t a;
 
        has_I16 = test16(b, &a) >= 0? 1 : 0;
-       rm_info(b, "I16="); // FIXME: probably this function has a bug. If I move it below, I16 will not be removed!
+       //rm_info(b, "I16="); // FIXME: probably this function has a bug. If I move it below, I16 will not be removed!
 
        memset(&s, 0, sizeof(kstring_t));
        kputc('\0', &s); kputs(b->ref, &s); kputc('\0', &s);
@@ -170,6 +154,8 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr,
        }
        if (has_I16 && a.is_tested) ksprintf(&s, ";PV4=%.2g,%.2g,%.2g,%.2g", a.p[0], a.p[1], a.p[2], a.p[3]);
        kputc('\0', &s);
+    rm_info(&s, "QS=");
+    rm_info(&s, "I16=");
        kputs(b->fmt, &s); kputc('\0', &s);
        free(b->str);
        b->m_str = s.m; b->l_str = s.l; b->str = s.s;
@@ -528,9 +514,9 @@ int bcfview(int argc, char *argv[])
                        int i;
                        for (i = 0; i < 9; ++i) em[i] = -1.;
                }
-        bcf_p1_set_ploidy(b, p1); // could be improved: do this per site to allow pseudo-autosomal regions
         if ( !(vc.flag&VC_KEEPALT) && vc.flag&VC_CALL && vc.min_ma_lrt>=0 )
         {
+            bcf_p1_set_ploidy(b, p1); // could be improved: do this per site to allow pseudo-autosomal regions
             int gts = call_multiallelic_gt(b,p1,vc.min_ma_lrt);
             if ( gts<=1 && vc.flag & VC_VARONLY ) continue;
         }

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index bc2aff8d1d65b9ef743d3aa9f8d60c86eb62996d..12894371c2eee9a6c5e05a5def471da413735f9b 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -201,6 +201,7 @@ void bcf_p1_destroy(bcf_p1aux_t *ma)
 }
 
 extern double kf_gammap(double s, double z);
+int test16(bcf1_t *b, anno16_t *a);
 
 int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
 {
@@ -222,7 +223,7 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
         exit(1);
     }
 
-    // set PL and PL_len
+    // find PL and DP FORMAT indexes
     uint8_t *pl = NULL;
     int npl = 0, idp=-1;
     int i;
@@ -237,16 +238,20 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
     }
     if ( !pl ) return -1;
 
+    assert(ma->q2p[0] == 1);
+
+    // Init P(D|G)
     int npdg = nals*(nals+1)/2;
-    float *pdg,*_pdg;
-    _pdg = pdg = malloc(sizeof(float)*ma->n*npdg);
+    double *pdg,*_pdg;
+    _pdg = pdg = malloc(sizeof(double)*ma->n*npdg);
     for (i=0; i<ma->n; i++)
     {
         int j; 
-        float sum = 0;
+        double sum = 0;
         for (j=0; j<npdg; j++)
         {
-            _pdg[j] = pow(10,-0.1*pl[j]);
+            //_pdg[j] = pow(10,-0.1*pl[j]); 
+            _pdg[j] = ma->q2p[pl[j]];
             sum += _pdg[j];
         }
         if ( sum )
@@ -256,89 +261,94 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
     }
 
     if ((p = strstr(b->info, "QS=")) == 0) { fprintf(stderr,"INFO/QS is required with -m, exiting\n"); exit(1); }
-    float qsum[4];
-    if ( sscanf(p+3,"%f,%f,%f,%f",&qsum[0],&qsum[1],&qsum[2],&qsum[3])!=4 ) { fprintf(stderr,"Could not parse %s\n",p); exit(1); }
+    double qsum[4];
+    if ( sscanf(p+3,"%lf,%lf,%lf,%lf",&qsum[0],&qsum[1],&qsum[2],&qsum[3])!=4 ) { fprintf(stderr,"Could not parse %s\n",p); exit(1); }
+
 
-    
+    // Calculate the most likely combination of alleles
     int ia,ib,ic, max_als=0, max_als2=0;
-    float max_lk = INT_MIN, max_lk2 = INT_MIN, lk_sum = INT_MIN;
+    double ref_lk = 0, max_lk = INT_MIN, max_lk2 = INT_MIN, lk_sum = INT_MIN;
     for (ia=0; ia<nals; ia++)
     {
-        //if ( ia && qsum[ia]==0 ) continue;
-        float lk_tot = 0;
+        double lk_tot = 0;
         int iaa = (ia+1)*(ia+2)/2-1;
         int isample;
         for (isample=0; isample<ma->n; isample++)
         {
-            float *p = pdg + isample*npdg;
-            assert( log(p[iaa]) <= 0 );
+            double *p = pdg + isample*npdg;
+            // assert( log(p[iaa]) <= 0 );
             lk_tot += log(p[iaa]);
         }
+        if ( ia==0 ) ref_lk = lk_tot;
         if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia; }
         else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia; }
         lk_sum = lk_tot>lk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum));
     }
-    for (ia=0; ia<nals; ia++)
+    if ( nals>1 )
     {
-        if ( qsum[ia]==0 ) continue;
-        //if ( ia && qsum[ia]==0 ) continue;
-        for (ib=0; ib<ia; ib++)
+        for (ia=0; ia<nals; ia++)
         {
-            if ( qsum[ib]==0 ) continue;
-            //if ( ib && qsum[ib]==0 ) continue;
-            //if ( qsum[ia]+qsum[ib]==0 ) continue;
-            float lk_tot = 0;
-            float fa  = qsum[ia]/(qsum[ia]+qsum[ib]);
-            float fb  = qsum[ib]/(qsum[ia]+qsum[ib]);
-            float fab = 2*fa*fb; fa *= fa; fb *= fb;
-            int isample, iaa = (ia+1)*(ia+2)/2-1, ibb = (ib+1)*(ib+2)/2-1, iab = iaa - ia + ib;
-            for (isample=0; isample<ma->n; isample++)
+            if ( qsum[ia]==0 ) continue;
+            int iaa = (ia+1)*(ia+2)/2-1;
+            for (ib=0; ib<ia; ib++)
             {
-                if ( b->ploidy && b->ploidy[isample]==1 ) continue;
-                float *p = pdg + isample*npdg;
-                assert( log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]) <= 0 );
-                lk_tot += log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]);
+                if ( qsum[ib]==0 ) continue;
+                double lk_tot = 0;
+                double fa  = qsum[ia]/(qsum[ia]+qsum[ib]);
+                double fb  = qsum[ib]/(qsum[ia]+qsum[ib]);
+                double fab = 2*fa*fb; fa *= fa; fb *= fb;
+                int isample, ibb = (ib+1)*(ib+2)/2-1, iab = iaa - ia + ib;
+                for (isample=0; isample<ma->n; isample++)
+                {
+                    if ( b->ploidy && b->ploidy[isample]==1 ) continue;
+                    double *p = pdg + isample*npdg;
+                    //assert( log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]) <= 0 );
+                    lk_tot += log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]);
+                }
+                if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia|1<<ib; }
+                else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia|1<<ib; }
+                lk_sum = lk_tot>lk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum));
             }
-            if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia|1<<ib; }
-            else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia|1<<ib; }
-            lk_sum = lk_tot>lk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum));
         }
     }
-    for (ia=0; ia<nals; ia++)
+    if ( nals>2 )
     {
-        if ( qsum[ia]==0 ) continue;
-        //if ( ia && qsum[ia]==0 ) continue;
-        for (ib=0; ib<ia; ib++)
+        for (ia=0; ia<nals; ia++)
         {
-            if ( qsum[ib]==0 ) continue;
-            //if ( ib && qsum[ib]==0 ) continue;
-            for (ic=0; ic<ib; ic++)
+            if ( qsum[ia]==0 ) continue;
+            int iaa = (ia+1)*(ia+2)/2-1;
+            for (ib=0; ib<ia; ib++)
             {
-                if ( qsum[ic]==0 ) continue;
-                //if ( ic && qsum[ic]==0 ) continue;
-                //if ( qsum[ia]+qsum[ib]+qsum[ic]==0 ) continue;
-                float lk_tot = 0;
-                float fa  = qsum[ia]/(qsum[ia]+qsum[ib]+qsum[ic]);
-                float fb  = qsum[ib]/(qsum[ia]+qsum[ib]+qsum[ic]);
-                float fc  = qsum[ic]/(qsum[ia]+qsum[ib]+qsum[ic]);
-                float fab = 2*fa*fb, fac = 2*fa*fc, fbc = 2*fb*fc; fa *= fa; fb *= fb; fc *= fc;
-                int isample, iaa = (ia+1)*(ia+2)/2-1, ibb = (ib+1)*(ib+2)/2-1, icc = (ic+1)*(ic+2)/2-1;
-                int iab = iaa - ia + ib, iac = iaa - ia + ic, ibc = ibb - ib + ic;
-                for (isample=0; isample<ma->n; isample++)
+                if ( qsum[ib]==0 ) continue;
+                int ibb = (ib+1)*(ib+2)/2-1; 
+                int iab = iaa - ia + ib;
+                for (ic=0; ic<ib; ic++)
                 {
-                    if ( b->ploidy && b->ploidy[isample]==1 ) continue;
-                    float *p = pdg + isample*npdg;
-                    assert( log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]) <= 0 );
-                    lk_tot += log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]);
+                    if ( qsum[ic]==0 ) continue;
+                    double lk_tot = 0;
+                    double fa  = qsum[ia]/(qsum[ia]+qsum[ib]+qsum[ic]);
+                    double fb  = qsum[ib]/(qsum[ia]+qsum[ib]+qsum[ic]);
+                    double fc  = qsum[ic]/(qsum[ia]+qsum[ib]+qsum[ic]);
+                    double fab = 2*fa*fb, fac = 2*fa*fc, fbc = 2*fb*fc; fa *= fa; fb *= fb; fc *= fc;
+                    int isample, icc = (ic+1)*(ic+2)/2-1;
+                    int iac = iaa - ia + ic, ibc = ibb - ib + ic;
+                    for (isample=0; isample<ma->n; isample++)
+                    {
+                        if ( b->ploidy && b->ploidy[isample]==1 ) continue;
+                        double *p = pdg + isample*npdg;
+                        //assert( log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]) <= 0 );
+                        lk_tot += log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]);
+                    }
+                    if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia|1<<ib|1<<ic; }
+                    else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia|1<<ib|1<<ic; }
+                    lk_sum = lk_tot>lk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum));
                 }
-                if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia|1<<ib|1<<ic; }
-                else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia|1<<ib|1<<ic; }
-                lk_sum = lk_tot>lk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum));
-        //printf("    %f\n", lk_sum);
             }
         }
     }
 
+
+    // Should we add another allele, does it increase the likelihood significantly?
     int n1=0, n2=0;
     for (i=0; i<nals; i++) if ( max_als&1<<i) n1++;
     for (i=0; i<nals; i++) if ( max_als2&1<<i) n2++;
@@ -361,14 +371,14 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
     for (isample = 0; isample < b->n_smpl; isample++) 
     {
         int ploidy = b->ploidy ? b->ploidy[isample] : 2;
-        float *p = pdg + isample*npdg;
+        double *p = pdg + isample*npdg;
         int ia, als = 0;
-        float lk = INT_MIN, lk_sum=0;
+        double lk = 0, lk_sum=0;
         for (ia=0; ia<nals; ia++)
         {
             if ( !(max_als&1<<ia) ) continue;
             int iaa = (ia+1)*(ia+2)/2-1;
-            float _lk = p[iaa]*qsum[ia]*qsum[ia];
+            double _lk = p[iaa]*qsum[ia]*qsum[ia];
             if ( _lk > lk ) { lk = _lk; als = ia<<3 | ia; }
             lk_sum += _lk;
         }
@@ -382,7 +392,7 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
                 {
                     if ( !(max_als&1<<ib) ) continue;
                     int iab = iaa - ia + ib;
-                    float _lk = 2*qsum[ia]*qsum[ib]*p[iab];
+                    double _lk = 2*qsum[ia]*qsum[ib]*p[iab];
                     if ( _lk > lk ) { lk = _lk; als = ib<<3 | ia; }
                     lk_sum += _lk;
                 }
@@ -404,42 +414,47 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
     }
     bcf_fit_alt(b,max_als);
 
-    // prepare ref, alt, filter, info, format
+
+    // Prepare BCF for output: ref, alt, filter, info, format
     memset(&s, 0, sizeof(kstring_t)); kputc('\0', &s); 
     kputs(b->ref, &s); kputc('\0', &s);
     kputs(b->alt, &s); kputc('\0', &s); kputc('\0', &s);
     {
-        kstring_t info; memset(&info, 0, sizeof(kstring_t));
         int an=0, nalts=0;
         for (i=0; i<nals; i++)
         {
             an += ac[i];
             if ( i>0 && ac[i] ) nalts++;
         }
-        ksprintf(&info, "AN=%d;", an);
+        ksprintf(&s, "AN=%d;", an);
         if ( nalts )
         {
-            kputs("AC=", &info);
+            kputs("AC=", &s);
             for (i=1; i<nals; i++)
             {
                 if ( !(gts&1<<i) ) continue;
                 nalts--;
-                ksprintf(&info,"%d", ac[i]);
-                if ( nalts>0 ) kputc(',', &info);
+                ksprintf(&s,"%d", ac[i]);
+                if ( nalts>0 ) kputc(',', &s);
             }
-            kputc(';', &info);
+            kputc(';', &s);
+        }
+        kputs(b->info, &s); 
+        anno16_t a;
+        int has_I16 = test16(b, &a) >= 0? 1 : 0;
+        if (has_I16 )
+        {
+            if ( a.is_tested) ksprintf(&s, ";PV4=%.2g,%.2g,%.2g,%.2g", a.p[0], a.p[1], a.p[2], a.p[3]);
+            ksprintf(&s, ";DP4=%d,%d,%d,%d;MQ=%d", a.d[0], a.d[1], a.d[2], a.d[3], a.mq);
         }
-        kputs(b->info, &info); 
-        info.l -= remove_tag(info.s, "I16=", ';');
-        info.l -= remove_tag(info.s, "QS=", ';');
-        kputc('\0', &info);
-        kputs(info.s, &s); kputc('\0', &s);
-        free(info.s);
+        rm_info(&s, "I16=");
+        rm_info(&s, "QS=");
+        kputc('\0', &s);
     }
     kputs(b->fmt, &s); kputc('\0', &s);
     free(b->str);
     b->m_str = s.m; b->l_str = s.l; b->str = s.s;
-    b->qual = -4.343*(log(1-exp(max_lk-lk_sum)));
+    b->qual = gts>1 ? -4.343*(ref_lk - lk_sum) : -4.343*(max_lk - lk_sum);
     if ( b->qual>999 ) b->qual = 999;
     bcf_sync(b);
 

diff --git a/bcftools/prob1.h b/bcftools/prob1.h
index 88ed5bf0861ad709364b507fdfa39a24b825eca0..eb0b145285715c54d997693c449b5262462cb904 100644 (file)
--- a/bcftools/prob1.h
+++ b/bcftools/prob1.h
@@ -14,6 +14,11 @@ typedef struct {
        double cmp[3], p_chi2, lrt; // used by contrast2()
 } bcf_p1rst_t;
 
+typedef struct {
+    double p[4];
+    int mq, depth, is_tested, d[4];
+} anno16_t;
+
 #define MC_PTYPE_FULL  1
 #define MC_PTYPE_COND2 2
 #define MC_PTYPE_FLAT  3