From faf76ae18e5ddf2d03253e2ed77f444bcdf4914c Mon Sep 17 00:00:00 2001 From: Petr Danecek Date: Thu, 5 Jul 2012 09:58:00 +0100 Subject: [PATCH] VCF annotations with and without -m and more robust rm_info --- bcftools/bcf.c | 27 ++++++-- bcftools/bcf.h | 4 +- bcftools/call1.c | 24 ++----- bcftools/prob1.c | 169 ++++++++++++++++++++++++++--------------------- bcftools/prob1.h | 5 ++ 5 files changed, 128 insertions(+), 101 deletions(-) diff --git a/bcftools/bcf.c b/bcftools/bcf.c index 87ff6c7..6695c74 100644 --- a/bcftools/bcf.c +++ b/bcftools/bcf.c @@ -323,7 +323,7 @@ int bcf_append_info(bcf1_t *b, const char *info, int l) return 0; } -int remove_tag(char *str, char *tag, char delim) +int remove_tag(char *str, const char *tag, char delim) { char *tmp = str, *p; int len_diff = 0, ori_len = strlen(str); @@ -336,11 +336,11 @@ int remove_tag(char *str, char *tag, char delim) } char *q=p+1; while ( *q && *q!=delim ) q++; - if ( p==str && *q ) q++; + if ( p==str && *q ) q++; // the tag is first, don't move the delim char len_diff += q-p; - if ( ! *q ) { *p = 0; break; } + if ( ! *q ) { *p = 0; break; } // the tag was last, no delim follows else - memmove(p,q,ori_len-(int)(q-p)+1); + memmove(p,q,ori_len-(int)(p-str)-(int)(q-p)-1); // *q==delim } if ( len_diff==ori_len ) str[0]='.', str[1]=0, len_diff--; @@ -348,6 +348,25 @@ int remove_tag(char *str, char *tag, char delim) return len_diff; } + +void rm_info(kstring_t *s, const char *key) +{ + char *p = s->s; + int n = 0; + while ( n<4 ) + { + if ( !*p ) n++; + p++; + } + char *q = p+1; + while ( *q && q-s->sl ) q++; + + int nrm = remove_tag(p, key, ';'); + if ( nrm ) + memmove(q-nrm, q, s->s+s->l-q+1); + s->l -= nrm; +} + int bcf_cpy(bcf1_t *r, const bcf1_t *b) { char *t1 = r->str; diff --git a/bcftools/bcf.h b/bcftools/bcf.h index 7e451a6..8c52451 100644 --- a/bcftools/bcf.h +++ b/bcftools/bcf.h @@ -124,7 +124,9 @@ extern "C" { // append more info int bcf_append_info(bcf1_t *b, const char *info, int l); // remove tag - int remove_tag(char *string, char *tag, char delim); + int remove_tag(char *string, const char *tag, char delim); + // remove info tag, string is the kstring holder of bcf1_t.str + void rm_info(kstring_t *string, const char *key); // copy int bcf_cpy(bcf1_t *r, const bcf1_t *b); diff --git a/bcftools/call1.c b/bcftools/call1.c index 415426a..22ff2ac 100644 --- a/bcftools/call1.c +++ b/bcftools/call1.c @@ -48,11 +48,6 @@ void *bed_read(const char *fn); void bed_destroy(void *_h); int bed_overlap(const void *_h, const char *chr, int beg, int end); -typedef struct { - double p[4]; - int mq, depth, is_tested, d[4]; -} anno16_t; - static double ttest(int n1, int n2, int a[4]) { extern double kf_betai(double a, double b, double x); @@ -83,7 +78,7 @@ static int test16_core(int anno[16], anno16_t *a) return 0; } -static int test16(bcf1_t *b, anno16_t *a) +int test16(bcf1_t *b, anno16_t *a) { char *p; int i, anno[16]; @@ -100,17 +95,6 @@ static int test16(bcf1_t *b, anno16_t *a) return test16_core(anno, a); } -static void rm_info(bcf1_t *b, const char *key) -{ - char *p, *q; - if ((p = strstr(b->info, key)) == 0) return; - for (q = p; *q && *q != ';'; ++q); - if (p > b->info && *(p-1) == ';') --p; - memmove(p, q, b->l_str - (q - b->str)); - b->l_str -= q - p; - bcf_sync(b); -} - static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10], int cons_llr, int64_t cons_gt) { kstring_t s; @@ -119,7 +103,7 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, anno16_t a; has_I16 = test16(b, &a) >= 0? 1 : 0; - rm_info(b, "I16="); // FIXME: probably this function has a bug. If I move it below, I16 will not be removed! + //rm_info(b, "I16="); // FIXME: probably this function has a bug. If I move it below, I16 will not be removed! memset(&s, 0, sizeof(kstring_t)); kputc('\0', &s); kputs(b->ref, &s); kputc('\0', &s); @@ -170,6 +154,8 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, } if (has_I16 && a.is_tested) ksprintf(&s, ";PV4=%.2g,%.2g,%.2g,%.2g", a.p[0], a.p[1], a.p[2], a.p[3]); kputc('\0', &s); + rm_info(&s, "QS="); + rm_info(&s, "I16="); kputs(b->fmt, &s); kputc('\0', &s); free(b->str); b->m_str = s.m; b->l_str = s.l; b->str = s.s; @@ -528,9 +514,9 @@ int bcfview(int argc, char *argv[]) int i; for (i = 0; i < 9; ++i) em[i] = -1.; } - bcf_p1_set_ploidy(b, p1); // could be improved: do this per site to allow pseudo-autosomal regions if ( !(vc.flag&VC_KEEPALT) && vc.flag&VC_CALL && vc.min_ma_lrt>=0 ) { + bcf_p1_set_ploidy(b, p1); // could be improved: do this per site to allow pseudo-autosomal regions int gts = call_multiallelic_gt(b,p1,vc.min_ma_lrt); if ( gts<=1 && vc.flag & VC_VARONLY ) continue; } diff --git a/bcftools/prob1.c b/bcftools/prob1.c index bc2aff8..1289437 100644 --- a/bcftools/prob1.c +++ b/bcftools/prob1.c @@ -201,6 +201,7 @@ void bcf_p1_destroy(bcf_p1aux_t *ma) } 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) { @@ -222,7 +223,7 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold) exit(1); } - // set PL and PL_len + // find PL and DP FORMAT indexes uint8_t *pl = NULL; int npl = 0, idp=-1; int i; @@ -237,16 +238,20 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold) } if ( !pl ) return -1; + assert(ma->q2p[0] == 1); + + // Init P(D|G) int npdg = nals*(nals+1)/2; - float *pdg,*_pdg; - _pdg = pdg = malloc(sizeof(float)*ma->n*npdg); + double *pdg,*_pdg; + _pdg = pdg = malloc(sizeof(double)*ma->n*npdg); for (i=0; in; i++) { int j; - float sum = 0; + double sum = 0; for (j=0; jq2p[pl[j]]; sum += _pdg[j]; } if ( sum ) @@ -256,89 +261,94 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold) } if ((p = strstr(b->info, "QS=")) == 0) { fprintf(stderr,"INFO/QS is required with -m, exiting\n"); exit(1); } - float qsum[4]; - if ( sscanf(p+3,"%f,%f,%f,%f",&qsum[0],&qsum[1],&qsum[2],&qsum[3])!=4 ) { fprintf(stderr,"Could not parse %s\n",p); 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; - float max_lk = INT_MIN, max_lk2 = INT_MIN, lk_sum = INT_MIN; + double ref_lk = 0, max_lk = INT_MIN, max_lk2 = INT_MIN, lk_sum = INT_MIN; for (ia=0; ian; isample++) { - float *p = pdg + isample*npdg; - assert( log(p[iaa]) <= 0 ); + 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)); } - for (ia=0; ia1 ) { - if ( qsum[ia]==0 ) continue; - //if ( ia && qsum[ia]==0 ) continue; - for (ib=0; ibn; isample++) + if ( qsum[ia]==0 ) continue; + int iaa = (ia+1)*(ia+2)/2-1; + for (ib=0; ibploidy && b->ploidy[isample]==1 ) continue; - float *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 ( 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; isamplen; 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 ( max_lklk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum)); } } - for (ia=0; ia2 ) { - if ( qsum[ia]==0 ) continue; - //if ( ia && qsum[ia]==0 ) continue; - for (ib=0; ibn; isample++) + if ( qsum[ib]==0 ) continue; + int ibb = (ib+1)*(ib+2)/2-1; + int iab = iaa - ia + ib; + for (ic=0; icploidy && b->ploidy[isample]==1 ) continue; - float *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 ( qsum[ic]==0 ) continue; + double lk_tot = 0; + double fa = qsum[ia]/(qsum[ia]+qsum[ib]+qsum[ic]); + double fb = qsum[ib]/(qsum[ia]+qsum[ib]+qsum[ic]); + double fc = qsum[ic]/(qsum[ia]+qsum[ib]+qsum[ic]); + double fab = 2*fa*fb, fac = 2*fa*fc, fbc = 2*fb*fc; fa *= fa; fb *= fb; fc *= fc; + int isample, icc = (ic+1)*(ic+2)/2-1; + int iac = iaa - ia + ic, ibc = ibb - ib + ic; + for (isample=0; isamplen; 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)); } - if ( max_lklk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum)); - //printf(" %f\n", lk_sum); } } } + + // Should we add another allele, does it increase the likelihood significantly? int n1=0, n2=0; for (i=0; in_smpl; isample++) { int ploidy = b->ploidy ? b->ploidy[isample] : 2; - float *p = pdg + isample*npdg; + double *p = pdg + isample*npdg; int ia, als = 0; - float lk = INT_MIN, lk_sum=0; + double lk = 0, lk_sum=0; for (ia=0; ia lk ) { lk = _lk; als = ia<<3 | ia; } lk_sum += _lk; } @@ -382,7 +392,7 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold) { if ( !(max_als&1< lk ) { lk = _lk; als = ib<<3 | ia; } lk_sum += _lk; } @@ -404,42 +414,47 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold) } bcf_fit_alt(b,max_als); - // prepare ref, alt, filter, info, format + + // 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); { - kstring_t info; memset(&info, 0, sizeof(kstring_t)); int an=0, nalts=0; for (i=0; i0 && ac[i] ) nalts++; } - ksprintf(&info, "AN=%d;", an); + ksprintf(&s, "AN=%d;", an); if ( nalts ) { - kputs("AC=", &info); + kputs("AC=", &s); for (i=1; i0 ) kputc(',', &info); + ksprintf(&s,"%d", ac[i]); + if ( nalts>0 ) kputc(',', &s); } - kputc(';', &info); + 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); } - kputs(b->info, &info); - info.l -= remove_tag(info.s, "I16=", ';'); - info.l -= remove_tag(info.s, "QS=", ';'); - kputc('\0', &info); - kputs(info.s, &s); kputc('\0', &s); - free(info.s); + rm_info(&s, "I16="); + rm_info(&s, "QS="); + kputc('\0', &s); } 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 = -4.343*(log(1-exp(max_lk-lk_sum))); + 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); diff --git a/bcftools/prob1.h b/bcftools/prob1.h index 88ed5bf..eb0b145 100644 --- a/bcftools/prob1.h +++ b/bcftools/prob1.h @@ -14,6 +14,11 @@ typedef struct { double cmp[3], p_chi2, lrt; // used by contrast2() } bcf_p1rst_t; +typedef struct { + double p[4]; + int mq, depth, is_tested, d[4]; +} anno16_t; + #define MC_PTYPE_FULL 1 #define MC_PTYPE_COND2 2 #define MC_PTYPE_FLAT 3 -- 2.39.2