X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=blobdiff_plain;f=bcftools%2Fprob1.c;h=e3d6b5eb42fee41bb3b2e8f392b5ae38d1934bc5;hp=a024d041ea80baeb3e0c75427c148aae0db728d0;hb=da45f88bac09817fcd87916d3e104f0e4bc53874;hpb=50706bea83d8a1518485283876333279f4ae7137 diff --git a/bcftools/prob1.c b/bcftools/prob1.c index a024d04..e3d6b5e 100644 --- a/bcftools/prob1.c +++ b/bcftools/prob1.c @@ -4,7 +4,9 @@ #include #include #include +#include #include "prob1.h" +#include "kstring.h" #include "kseq.h" KSTREAM_INIT(gzFile, gzread, 16384) @@ -40,6 +42,7 @@ struct __bcf_p1aux_t { double *z, *zswap; // aux for afs double *z1, *z2, *phi1, *phi2; // only calculated when n1 is set double **hg; // hypergeometric distribution + double *lf; // log factorial double t, t1, t2; double *afs, *afs1; // afs: accumulative AFS; afs1: site posterior distribution const uint8_t *PL; // point to PL @@ -154,8 +157,10 @@ bcf_p1aux_t *bcf_p1_init(int n, uint8_t *ploidy) ma->z2 = calloc(ma->M + 1, sizeof(double)); ma->afs = calloc(ma->M + 1, sizeof(double)); ma->afs1 = calloc(ma->M + 1, sizeof(double)); + ma->lf = calloc(ma->M + 1, sizeof(double)); for (i = 0; i < 256; ++i) ma->q2p[i] = pow(10., -i / 10.); + for (i = 0; i <= ma->M; ++i) ma->lf[i] = lgamma(i + 1); bcf_p1_init_prior(ma, MC_PTYPE_FULL, 1e-3); // the simplest prior return ma; } @@ -171,10 +176,18 @@ int bcf_p1_set_n1(bcf_p1aux_t *b, int n1) return 0; } +void bcf_p1_set_ploidy(bcf1_t *b, bcf_p1aux_t *ma) +{ + // bcf_p1aux_t fields are not visible outside of prob1.c, hence this wrapper. + // Ideally, this should set ploidy per site to allow pseudo-autosomal regions + b->ploidy = ma->ploidy; +} + void bcf_p1_destroy(bcf_p1aux_t *ma) { if (ma) { int k; + free(ma->lf); if (ma->hg && ma->n1 > 0) { for (k = 0; k <= 2*ma->n1; ++k) free(ma->hg[k]); free(ma->hg); @@ -187,61 +200,292 @@ void bcf_p1_destroy(bcf_p1aux_t *ma) } } -static int cal_pdg(const bcf1_t *b, bcf_p1aux_t *ma) -{ - int i, j; - long *p, tmp; - p = alloca(b->n_alleles * sizeof(long)); - memset(p, 0, sizeof(long) * b->n_alleles); - for (j = 0; j < ma->n; ++j) { - const uint8_t *pi = ma->PL + j * ma->PL_len; - double *pdg = ma->pdg + j * 3; - pdg[0] = ma->q2p[pi[2]]; pdg[1] = ma->q2p[pi[1]]; pdg[2] = ma->q2p[pi[0]]; - for (i = 0; i < b->n_alleles; ++i) - p[i] += (int)pi[(i+1)*(i+2)/2-1]; - } - for (i = 0; i < b->n_alleles; ++i) p[i] = p[i]<<4 | i; - for (i = 1; i < b->n_alleles; ++i) // insertion sort - for (j = i; j > 0 && p[j] < p[j-1]; --j) - tmp = p[j], p[j] = p[j-1], p[j-1] = tmp; - for (i = b->n_alleles - 1; i >= 0; --i) - if ((p[i]&0xf) == 0) break; - return i; -} -// f0 is the reference allele frequency -static double mc_freq_iter(double f0, const bcf_p1aux_t *ma, int beg, int end) +extern double kf_gammap(double s, double z); +int test16(bcf1_t *b, anno16_t *a); + +int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold) { - double f, f3[3]; - int i; - f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0; - for (i = beg, f = 0.; i < end; ++i) { - double *pdg; - pdg = ma->pdg + i * 3; - f += (pdg[1] * f3[1] + 2. * pdg[2] * f3[2]) - / (pdg[0] * f3[0] + pdg[1] * f3[1] + pdg[2] * f3[2]); - } - f /= (end - beg) * 2.; - return f; + 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; in; i++) + { + int j; + double sum = 0; + for (j=0; jq2p[pl[j]]; + sum += _pdg[j]; + } + if ( sum ) + for (j=0; jinfo, "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; ian; 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_lklk_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; ian; 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_lklk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum)); + } + } + } + if ( nals>2 ) + { + for (ia=0; ian; 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 ( max_lklk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum)); + } + } + } + } + + + // Should we add another allele, does it increase the likelihood significantly? + int n1=0, n2=0; + for (i=0; in_gi; + s.m = b->m_str; s.l = b->l_str - 1; s.s = b->str; + kputs(":GT:GQ", &s); kputc('\0', &s); + b->m_str = s.m; b->l_str = s.l; b->str = s.s; + bcf_sync(b); + + // Call GTs + int isample, gts=0, ac[4] = {0,0,0,0}; + for (isample = 0; isample < b->n_smpl; isample++) + { + int ploidy = b->ploidy ? b->ploidy[isample] : 2; + double *p = pdg + isample*npdg; + int ia, als = 0; + double lk = 0, lk_sum=0; + for (ia=0; ia lk ) { lk = _lk; als = ia<<3 | ia; } + lk_sum += _lk; + } + if ( ploidy==2 ) + { + for (ia=0; ia lk ) { lk = _lk; als = ib<<3 | ia; } + lk_sum += _lk; + } + } + } + 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] = 1<<7; + ((uint8_t*)b->gi[old_n_gi+1].data)[isample] = 0; + continue; + } + ((uint8_t*)b->gi[old_n_gi].data)[isample] = als; + ((uint8_t*)b->gi[old_n_gi+1].data)[isample] = lk<100 ? (int)lk : 99; + + gts |= 1<<(als>>3&7) | 1<<(als&7); + ac[ als>>3&7 ]++; + ac[ als&7 ]++; + } + bcf_fit_alt(b,max_als); + + + // Prepare BCF for output: ref, alt, filter, info, format + memset(&s, 0, sizeof(kstring_t)); kputc('\0', &s); + kputs(b->ref, &s); kputc('\0', &s); + kputs(b->alt, &s); kputc('\0', &s); kputc('\0', &s); + { + int an=0, nalts=0; + for (i=0; i0 && ac[i] ) nalts++; + } + ksprintf(&s, "AN=%d;", an); + if ( nalts ) + { + kputs("AC=", &s); + for (i=1; i0 ) kputc(',', &s); + } + kputc(';', &s); + } + kputs(b->info, &s); + anno16_t a; + int has_I16 = test16(b, &a) >= 0? 1 : 0; + if (has_I16 ) + { + 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="); + } + kputs(b->fmt, &s); kputc('\0', &s); + free(b->str); + b->m_str = s.m; b->l_str = s.l; b->str = s.s; + b->qual = gts>1 ? -4.343*(ref_lk - lk_sum) : -4.343*(max_lk - lk_sum); + if ( b->qual>999 ) b->qual = 999; + bcf_sync(b); + + + free(pdg); + return gts; } -static double mc_gtfreq_iter(double g[3], const bcf_p1aux_t *ma, int beg, int end) +static int cal_pdg(const bcf1_t *b, bcf_p1aux_t *ma) { - double err, gg[3]; - int i; - gg[0] = gg[1] = gg[2] = 0.; - for (i = beg; i < end; ++i) { - double *pdg, sum, tmp[3]; - pdg = ma->pdg + i * 3; - tmp[0] = pdg[0] * g[0]; tmp[1] = pdg[1] * g[1]; tmp[2] = pdg[2] * g[2]; - sum = (tmp[0] + tmp[1] + tmp[2]) * (end - beg); - gg[0] += tmp[0] / sum; gg[1] += tmp[1] / sum; gg[2] += tmp[2] / sum; - } - err = fabs(gg[0] - g[0]) > fabs(gg[1] - g[1])? fabs(gg[0] - g[0]) : fabs(gg[1] - g[1]); - err = err > fabs(gg[2] - g[2])? err : fabs(gg[2] - g[2]); - g[0] = gg[0]; g[1] = gg[1]; g[2] = gg[2]; - return err; + int i, j; + long *p, tmp; + p = alloca(b->n_alleles * sizeof(long)); + memset(p, 0, sizeof(long) * b->n_alleles); + for (j = 0; j < ma->n; ++j) { + const uint8_t *pi = ma->PL + j * ma->PL_len; + double *pdg = ma->pdg + j * 3; + pdg[0] = ma->q2p[pi[2]]; pdg[1] = ma->q2p[pi[1]]; pdg[2] = ma->q2p[pi[0]]; + for (i = 0; i < b->n_alleles; ++i) + p[i] += (int)pi[(i+1)*(i+2)/2-1]; + } + for (i = 0; i < b->n_alleles; ++i) p[i] = p[i]<<4 | i; + for (i = 1; i < b->n_alleles; ++i) // insertion sort + for (j = i; j > 0 && p[j] < p[j-1]; --j) + tmp = p[j], p[j] = p[j-1], p[j-1] = tmp; + for (i = b->n_alleles - 1; i >= 0; --i) + if ((p[i]&0xf) == 0) break; + return i; } + int bcf_p1_call_gt(const bcf_p1aux_t *ma, double f0, int k) { double sum, g[3]; @@ -385,9 +629,10 @@ static inline double chi2_test(int a, int b, int c, int d) } // chi2=(a+b+c+d)(ad-bc)^2/[(a+b)(c+d)(a+c)(b+d)] -static inline double contrast2_aux(const bcf_p1aux_t *p1, double sum, int n1, int n2, int k1, int k2, double x[3]) +static inline double contrast2_aux(const bcf_p1aux_t *p1, double sum, int k1, int k2, double x[3]) { double p = p1->phi[k1+k2] * p1->z1[k1] * p1->z2[k2] / sum * p1->hg[k1][k2]; + int n1 = p1->n1, n2 = p1->n - p1->n1; if (p < CONTRAST_TINY) return -1; if (.5*k1/n1 < .5*k2/n2) x[1] += p; else if (.5*k1/n1 > .5*k2/n2) x[2] += p; @@ -415,12 +660,12 @@ static double contrast2(bcf_p1aux_t *p1, double ret[3]) p1->hg[k1][k2] = exp(lgamma(k1+k2+1) + lgamma(p1->M-k1-k2+1) - (lgamma(k1+1) + lgamma(k2+1) + lgamma(2*n1-k1+1) + lgamma(2*n2-k2+1) + tmp)); } } - { // compute sum1 and sum2 + { // compute long double suml = 0; for (k = 0; k <= p1->M; ++k) suml += p1->phi[k] * p1->z[k]; sum = suml; } - { // get the mean k1 and k2 + { // get the max k1 and k2 double max; int max_k; for (k = 0, max = 0, max_k = -1; k <= 2*n1; ++k) { @@ -436,15 +681,15 @@ static double contrast2(bcf_p1aux_t *p1, double ret[3]) } { // We can do the following with one nested loop, but that is an O(N^2) thing. The following code block is much faster for large N. double x[3], y; - long double z = 0.; - x[0] = x[1] = x[2] = 0; + long double z = 0., L[2]; + x[0] = x[1] = x[2] = 0; L[0] = L[1] = 0; for (k1 = k10; k1 >= 0; --k1) { for (k2 = k20; k2 >= 0; --k2) { - if ((y = contrast2_aux(p1, sum, n1, n2, k1, k2, x)) < 0) break; + if ((y = contrast2_aux(p1, sum, k1, k2, x)) < 0) break; else z += y; } for (k2 = k20 + 1; k2 <= 2*n2; ++k2) { - if ((y = contrast2_aux(p1, sum, n1, n2, k1, k2, x)) < 0) break; + if ((y = contrast2_aux(p1, sum, k1, k2, x)) < 0) break; else z += y; } } @@ -452,21 +697,21 @@ static double contrast2(bcf_p1aux_t *p1, double ret[3]) x[0] = x[1] = x[2] = 0; for (k1 = k10 + 1; k1 <= 2*n1; ++k1) { for (k2 = k20; k2 >= 0; --k2) { - if ((y = contrast2_aux(p1, sum, n1, n2, k1, k2, x)) < 0) break; + if ((y = contrast2_aux(p1, sum, k1, k2, x)) < 0) break; else z += y; } for (k2 = k20 + 1; k2 <= 2*n2; ++k2) { - if ((y = contrast2_aux(p1, sum, n1, n2, k1, k2, x)) < 0) break; + if ((y = contrast2_aux(p1, sum, k1, k2, x)) < 0) break; else z += y; } } ret[0] += x[0]; ret[1] += x[1]; ret[2] += x[2]; - if (ret[0] + ret[1] + ret[2] < 0.99) { // in case of bad things happened - ret[0] = ret[1] = ret[2] = 0; + if (ret[0] + ret[1] + ret[2] < 0.95) { // in case of bad things happened + ret[0] = ret[1] = ret[2] = 0; L[0] = L[1] = 0; for (k1 = 0, z = 0.; k1 <= 2*n1; ++k1) for (k2 = 0; k2 <= 2*n2; ++k2) - if ((y = contrast2_aux(p1, sum, n1, n2, k1, k2, ret)) >= 0) z += y; - if (ret[0] + ret[1] + ret[2] < 0.99) // It seems that this may be caused by floating point errors. I do not really understand why... + if ((y = contrast2_aux(p1, sum, k1, k2, ret)) >= 0) z += y; + if (ret[0] + ret[1] + ret[2] < 0.95) // It seems that this may be caused by floating point errors. I do not really understand why... z = 1.0, ret[0] = ret[1] = ret[2] = 1./3; } return (double)z; @@ -502,7 +747,7 @@ static double mc_cal_afs(bcf_p1aux_t *ma, double *p_ref_folded, double *p_var_fo return sum / ma->M; } -int bcf_p1_cal(const bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst) +int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *rst) { int i, k; long double sum = 0.; @@ -516,6 +761,7 @@ int bcf_p1_cal(const bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst) break; } } + if (i == b->n_gi) return -1; // no PL if (b->n_alleles < 2) return -1; // FIXME: find a better solution // rst->rank0 = cal_pdg(b, ma); @@ -524,6 +770,13 @@ int bcf_p1_cal(const bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst) for (k = 0, sum = 0.; k < ma->M; ++k) sum += ma->afs1[k]; rst->p_var = (double)sum; + { // compute the allele count + double max = -1; + rst->ac = -1; + for (k = 0; k <= ma->M; ++k) + if (max < ma->z[k]) max = ma->z[k], rst->ac = k; + rst->ac = ma->M - rst->ac; + } // calculate f_flat and f_em for (k = 0, sum = 0.; k <= ma->M; ++k) sum += (long double)ma->z[k]; @@ -533,46 +786,32 @@ int bcf_p1_cal(const bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst) rst->f_flat += k * p; } rst->f_flat /= ma->M; - { // calculate f_em - double flast = rst->f_flat; - for (i = 0; i < MC_MAX_EM_ITER; ++i) { - rst->f_em = mc_freq_iter(flast, ma, 0, ma->n); - if (fabs(rst->f_em - flast) < MC_EM_EPS) break; - flast = rst->f_em; - } - if (ma->n1 > 0 && ma->n1 < ma->n) { - for (k = 0; k < 2; ++k) { - flast = rst->f_em; - for (i = 0; i < MC_MAX_EM_ITER; ++i) { - rst->f_em2[k] = k? mc_freq_iter(flast, ma, ma->n1, ma->n) : mc_freq_iter(flast, ma, 0, ma->n1); - if (fabs(rst->f_em2[k] - flast) < MC_EM_EPS) break; - flast = rst->f_em2[k]; - } - } - } - } - { // compute g[3] - rst->g[0] = (1. - rst->f_em) * (1. - rst->f_em); - rst->g[1] = 2. * rst->f_em * (1. - rst->f_em); - rst->g[2] = rst->f_em * rst->f_em; - for (i = 0; i < MC_MAX_EM_ITER; ++i) - if (mc_gtfreq_iter(rst->g, ma, 0, ma->n) < MC_EM_EPS) break; - } { // estimate equal-tail credible interval (95% level) int l, h; double p; - for (i = 0, p = 0.; i < ma->M; ++i) + for (i = 0, p = 0.; i <= ma->M; ++i) if (p + ma->afs1[i] > 0.025) break; else p += ma->afs1[i]; l = i; - for (i = ma->M-1, p = 0.; i >= 0; --i) + for (i = ma->M, p = 0.; i >= 0; --i) if (p + ma->afs1[i] > 0.025) break; else p += ma->afs1[i]; h = i; rst->cil = (double)(ma->M - h) / ma->M; rst->cih = (double)(ma->M - l) / ma->M; } + if (ma->n1 > 0) { // compute LRT + double max0, max1, max2; + for (k = 0, max0 = -1; k <= ma->M; ++k) + if (max0 < ma->z[k]) max0 = ma->z[k]; + for (k = 0, max1 = -1; k <= ma->n1 * 2; ++k) + if (max1 < ma->z1[k]) max1 = ma->z1[k]; + for (k = 0, max2 = -1; k <= ma->M - ma->n1 * 2; ++k) + if (max2 < ma->z2[k]) max2 = ma->z2[k]; + rst->lrt = log(max1 * max2 / max0); + rst->lrt = rst->lrt < 0? 1 : kf_gammaq(.5, rst->lrt); + } else rst->lrt = -1.0; rst->cmp[0] = rst->cmp[1] = rst->cmp[2] = rst->p_chi2 = -1.0; - if (rst->p_var > 0.1) // skip contrast2() if the locus is a strong non-variant + if (do_contrast && rst->p_var > 0.5) // skip contrast2() if the locus is a strong non-variant rst->p_chi2 = contrast2(ma, rst->cmp); return 0; }