+ int nals = 1;
+ char *p;
+ for (p=b->alt; *p; p++)
+ {
+ if ( *p=='X' || p[0]=='.' ) break;
+ if ( p[0]==',' ) nals++;
+ }
+ if ( b->alt[0] && !*p ) nals++;
+
+ if ( nals==1 ) return 1;
+
+ if ( nals>4 )
+ {
+ if ( *b->ref=='N' ) return 0;
+ fprintf(stderr,"Not ready for this, more than 4 alleles at %d: %s, %s\n", b->pos+1, b->ref,b->alt);
+ exit(1);
+ }
+
+ // find PL and DP FORMAT indexes
+ uint8_t *pl = NULL;
+ int npl = 0, idp=-1;
+ int i;
+ for (i = 0; i < b->n_gi; ++i)
+ {
+ if (b->gi[i].fmt == bcf_str2int("PL", 2))
+ {
+ pl = (uint8_t*)b->gi[i].data;
+ npl = b->gi[i].len;
+ }
+ if (b->gi[i].fmt == bcf_str2int("DP", 2)) idp=i;
+ }
+ if ( !pl ) return -1;
+
+ assert(ma->q2p[0] == 1);
+
+ // Init P(D|G)
+ int npdg = nals*(nals+1)/2;
+ double *pdg,*_pdg;
+ _pdg = pdg = malloc(sizeof(double)*ma->n*npdg);
+ for (i=0; i<ma->n; i++)
+ {
+ int j;
+ double sum = 0;
+ for (j=0; j<npdg; j++)
+ {
+ //_pdg[j] = pow(10,-0.1*pl[j]);
+ _pdg[j] = ma->q2p[pl[j]];
+ sum += _pdg[j];
+ }
+ if ( sum )
+ for (j=0; j<npdg; j++) _pdg[j] /= sum;
+ _pdg += npdg;
+ pl += npl;
+ }
+
+ if ((p = strstr(b->info, "QS=")) == 0) { fprintf(stderr,"INFO/QS is required with -m, exiting\n"); 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;
+ double ref_lk = 0, max_lk = INT_MIN, max_lk2 = INT_MIN, lk_sum = INT_MIN;
+ for (ia=0; ia<nals; ia++)
+ {
+ double lk_tot = 0;
+ int iaa = (ia+1)*(ia+2)/2-1;
+ int isample;
+ for (isample=0; isample<ma->n; isample++)
+ {
+ 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));
+ }
+ if ( nals>1 )
+ {
+ for (ia=0; ia<nals; ia++)
+ {
+ if ( qsum[ia]==0 ) continue;
+ int iaa = (ia+1)*(ia+2)/2-1;
+ for (ib=0; ib<ia; ib++)
+ {
+ 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));
+ }
+ }