X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fprob1.c;h=ffa608e29d6408c03a858f3db3daba162a1e9ed2;hb=5193fbf4944f67fd2d2bb0be5a54b5bbd27c5e6b;hp=e3d6b5eb42fee41bb3b2e8f392b5ae38d1934bc5;hpb=da45f88bac09817fcd87916d3e104f0e4bc53874;p=samtools.git

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index e3d6b5e..ffa608e 100644
--- 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;
@@ -300,10 +305,12 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
                 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 ( b->ploidy && b->ploidy[isample]==1 )
+                        lk_tot +=  log(fa*p[iaa] + fb*p[ibb]);
+                    else 
+                        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; }
@@ -334,10 +341,12 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
                     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 ( b->ploidy && b->ploidy[isample]==1 ) 
+                            lk_tot += log(fa*p[iaa] + fb*p[ibb] + fc*p[icc]);
+                        else
+                            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; }
@@ -347,7 +356,6 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
         }
     }
 
-
     // 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++;
@@ -595,6 +603,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)