mpileup: New HWE, QBD, RPS annotations. Modified VDB (beware: this version of VDB...

[samtools.git] / bcftools / prob1.c

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index d6557223e729eddbf0b2f70bd561c2d3d326143c..3539ee3430ddbec29795354c1924a32638db2eba 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -5,6 +5,7 @@
 #include <errno.h>
 #include <assert.h>
 #include <limits.h>
+#include <zlib.h>
 #include "prob1.h"
 #include "kstring.h"
 
@@ -15,6 +16,8 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define MC_EM_EPS 1e-5
 #define MC_DEF_INDEL 0.15
 
+gzFile bcf_p1_fp_lk;
+
 unsigned char seq_nt4_table[256] = {
        4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
        4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
@@ -165,6 +168,8 @@ bcf_p1aux_t *bcf_p1_init(int n, uint8_t *ploidy)
        return ma;
 }
 
+int bcf_p1_get_M(bcf_p1aux_t *b) { return b->M; }
+
 int bcf_p1_set_n1(bcf_p1aux_t *b, int n1)
 {
        if (n1 == 0 || n1 >= b->n) return -1;
@@ -497,7 +502,7 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold, int var_o
 
     // Call GTs
     int isample, gts=0, ac[4] = {0,0,0,0};
-    int nRR = 0, nAA = 0, nRA = 0, max_dv = 0, dp_nref = 0;
+    int nRR = 0, nAA = 0, nRA = 0, max_dv = 0;
     for (isample = 0; isample < b->n_smpl; isample++) 
     {
         int ploidy = b->ploidy ? b->ploidy[isample] : 2;
@@ -546,7 +551,6 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold, int var_o
         if ( als && idv>=0 && (dv=((uint16_t*)b->gi[idv].data)[isample]) )
         {
             if ( max_dv < dv ) max_dv = dv;
-            dp_nref += dp;
         }
 
         // For HWE annotation; multiple ALT alleles treated as one
@@ -599,7 +603,6 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold, int var_o
             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);
             ksprintf(&s, ";QBD=%e", b->qual/(a.d[0] + a.d[1] + a.d[2] + a.d[3]));
-            if ( dp_nref ) ksprintf(&s, ";QBDNR=%e", b->qual/dp_nref);
             if ( max_dv ) ksprintf(&s, ";MDV=%d", max_dv);
         }
         if ( nAA+nRA )
@@ -751,6 +754,8 @@ static void mc_cal_y_core(bcf_p1aux_t *ma, int beg)
                }
        }
        if (z[0] != ma->z) memcpy(ma->z, z[0], sizeof(double) * (ma->M + 1));
+       if (bcf_p1_fp_lk)
+               gzwrite(bcf_p1_fp_lk, ma->z, sizeof(double) * (ma->M + 1));
 }
 
 static void mc_cal_y(bcf_p1aux_t *ma)