#include <errno.h>
#include <assert.h>
#include <limits.h>
+#include <zlib.h>
#include "prob1.h"
#include "kstring.h"
#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,
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;
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; }
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; }
}
}
-
// 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++;
lk = -log(1-lk/lk_sum)/0.2302585;
if ( idp>=0 && ((uint16_t*)b->gi[idp].data)[isample]==0 )
{
- ((uint8_t*)b->gi[old_n_gi].data)[isample] = als | 1<<7;
+ ((uint8_t*)b->gi[old_n_gi].data)[isample] = 1<<7;
((uint8_t*)b->gi[old_n_gi+1].data)[isample] = 0;
continue;
}
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);
}
+ kputc('\0', &s);
rm_info(&s, "I16=");
rm_info(&s, "QS=");
- kputc('\0', &s);
}
kputs(b->fmt, &s); kputc('\0', &s);
free(b->str);
}
}
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)