bam_header_t *h_tmp;
data[i] = calloc(1, sizeof(mplp_aux_t));
data[i]->fp = strcmp(fn[i], "-") == 0? bam_dopen(fileno(stdin), "r") : bam_open(fn[i], "r");
+ if ( !data[i]->fp )
+ {
+ fprintf(stderr, "[%s] failed to open %d-th input.\n", __func__, i+1);
+ exit(1);
+ }
data[i]->conf = conf;
h_tmp = bam_header_read(data[i]->fp);
data[i]->h = i? h : h_tmp; // for i==0, "h" has not been set yet
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++;