#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;
}
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)
+// Wigginton 2005, PMID: 15789306
+// written by Jan Wigginton
+double calc_hwe(int obs_hom1, int obs_hom2, int obs_hets)
+{
+ if (obs_hom1 + obs_hom2 + obs_hets == 0 ) return 1;
+
+ assert(obs_hom1 >= 0 && obs_hom2 >= 0 && obs_hets >= 0);
+
+ int obs_homc = obs_hom1 < obs_hom2 ? obs_hom2 : obs_hom1;
+ int obs_homr = obs_hom1 < obs_hom2 ? obs_hom1 : obs_hom2;
+
+ int rare_copies = 2 * obs_homr + obs_hets;
+ int genotypes = obs_hets + obs_homc + obs_homr;
+
+ double *het_probs = (double*) calloc(rare_copies+1, sizeof(double));
+
+ /* start at midpoint */
+ int mid = rare_copies * (2 * genotypes - rare_copies) / (2 * genotypes);
+
+ /* check to ensure that midpoint and rare alleles have same parity */
+ if ((rare_copies & 1) ^ (mid & 1)) mid++;
+
+ int curr_hets = mid;
+ int curr_homr = (rare_copies - mid) / 2;
+ int curr_homc = genotypes - curr_hets - curr_homr;
+
+ het_probs[mid] = 1.0;
+ double sum = het_probs[mid];
+ for (curr_hets = mid; curr_hets > 1; curr_hets -= 2)
+ {
+ het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0) / (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0));
+ sum += het_probs[curr_hets - 2];
+
+ /* 2 fewer heterozygotes for next iteration -> add one rare, one common homozygote */
+ curr_homr++;
+ curr_homc++;
+ }
+
+ curr_hets = mid;
+ curr_homr = (rare_copies - mid) / 2;
+ curr_homc = genotypes - curr_hets - curr_homr;
+ for (curr_hets = mid; curr_hets <= rare_copies - 2; curr_hets += 2)
+ {
+ het_probs[curr_hets + 2] = het_probs[curr_hets] * 4.0 * curr_homr * curr_homc /((curr_hets + 2.0) * (curr_hets + 1.0));
+ sum += het_probs[curr_hets + 2];
+
+ /* add 2 heterozygotes for next iteration -> subtract one rare, one common homozygote */
+ curr_homr--;
+ curr_homc--;
+ }
+ int i;
+ for (i = 0; i <= rare_copies; i++) het_probs[i] /= sum;
+
+ /* p-value calculation for p_hwe */
+ double p_hwe = 0.0;
+ for (i = 0; i <= rare_copies; i++)
+ {
+ if (het_probs[i] > het_probs[obs_hets])
+ continue;
+ p_hwe += het_probs[i];
+ }
+
+ p_hwe = p_hwe > 1.0 ? 1.0 : p_hwe;
+ free(het_probs);
+ return p_hwe;
+
+}
+
+
+static void _bcf1_set_ref(bcf1_t *b, int idp)
+{
+ kstring_t s;
+ int old_n_gi = b->n_gi;
+ s.m = b->m_str; s.l = b->l_str - 1; s.s = b->str;
+ kputs(":GT", &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, an = 0;
+ for (isample = 0; isample < b->n_smpl; isample++)
+ {
+ if ( idp>=0 && ((uint16_t*)b->gi[idp].data)[isample]==0 )
+ ((uint8_t*)b->gi[old_n_gi].data)[isample] = 1<<7;
+ else
+ {
+ ((uint8_t*)b->gi[old_n_gi].data)[isample] = 0;
+ an += b->ploidy ? b->ploidy[isample] : 2;
+ }
+ }
+ bcf_fit_alt(b,1);
+ b->qual = 999;
+
+ // 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);
+ {
+ ksprintf(&s, "AN=%d;", an);
+ 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;
+ bcf_sync(b);
+}
+
+int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold, int var_only)
{
int nals = 1;
char *p;
}
if ( b->alt[0] && !*p ) nals++;
- if ( nals==1 ) return 1;
-
- if ( nals>4 ) { fprintf(stderr,"too many alts: %d\n", nals); exit(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);
+ }
- // set PL and PL_len
+ // find PL, DV and DP FORMAT indexes
uint8_t *pl = NULL;
- int npl = 0, idp=-1;
- int i;
+ int i, npl = 0, idp = -1, idv = -1;
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;
+ else if (b->gi[i].fmt == bcf_str2int("DP", 2)) idp=i;
+ else if (b->gi[i].fmt == bcf_str2int("DV", 2)) idv=i;
+ }
+ if ( nals==1 )
+ {
+ if ( !var_only ) _bcf1_set_ref(b, idp);
+ return 1;
}
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; i<ma->n; i++)
{
int j;
- float sum = 0;
+ double sum = 0;
for (j=0; j<npdg; j++)
{
- _pdg[j] = pow(10,-0.1*pl[j]);
+ //_pdg[j] = pow(10,-0.1*pl[j]);
+ _pdg[j] = ma->q2p[pl[j]];
sum += _pdg[j];
}
if ( sum )
}
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, remembering the most and second most likely set
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, lk_sums[3];
for (ia=0; ia<nals; ia++)
{
- //if ( ia && qsum[ia]==0 ) continue;
- float lk_tot = 0;
+ double lk_tot = 0;
int iaa = (ia+1)*(ia+2)/2-1;
int isample;
for (isample=0; isample<ma->n; 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_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));
}
- for (ia=0; ia<nals; ia++)
+ lk_sums[0] = lk_sum;
+ if ( nals>1 )
{
- if ( qsum[ia]==0 ) continue;
- //if ( ia && qsum[ia]==0 ) continue;
- for (ib=0; ib<ia; ib++)
+ for (ia=0; ia<nals; ia++)
{
- if ( qsum[ib]==0 ) continue;
- //if ( ib && qsum[ib]==0 ) continue;
- //if ( qsum[ia]+qsum[ib]==0 ) continue;
- float lk_tot = 0;
- float fa = qsum[ia]/(qsum[ia]+qsum[ib]);
- float fb = qsum[ib]/(qsum[ia]+qsum[ib]);
- float fab = 2*fa*fb; fa *= fa; fb *= fb;
- int isample, iaa = (ia+1)*(ia+2)/2-1, ibb = (ib+1)*(ib+2)/2-1, iab = iaa - ia + ib;
- for (isample=0; isample<ma->n; isample++)
+ if ( qsum[ia]==0 ) continue;
+ int iaa = (ia+1)*(ia+2)/2-1;
+ for (ib=0; ib<ia; ib++)
{
- if ( b->ploidy && 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; isample<ma->n; isample++)
+ {
+ double *p = pdg + isample*npdg;
+ //assert( log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]) <= 0 );
+ 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; }
+ 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 ( 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));
}
+ lk_sums[1] = lk_sum;
}
- for (ia=0; ia<nals; ia++)
+ if ( nals>2 )
{
- if ( qsum[ia]==0 ) continue;
- //if ( ia && qsum[ia]==0 ) continue;
- for (ib=0; ib<ia; ib++)
+ for (ia=0; ia<nals; ia++)
{
- if ( qsum[ib]==0 ) continue;
- //if ( ib && qsum[ib]==0 ) continue;
- for (ic=0; ic<ib; ic++)
+ if ( qsum[ia]==0 ) continue;
+ int iaa = (ia+1)*(ia+2)/2-1;
+ for (ib=0; ib<ia; ib++)
{
- if ( qsum[ic]==0 ) continue;
- //if ( ic && qsum[ic]==0 ) continue;
- //if ( qsum[ia]+qsum[ib]+qsum[ic]==0 ) continue;
- float lk_tot = 0;
- float fa = qsum[ia]/(qsum[ia]+qsum[ib]+qsum[ic]);
- float fb = qsum[ib]/(qsum[ia]+qsum[ib]+qsum[ic]);
- float fc = qsum[ic]/(qsum[ia]+qsum[ib]+qsum[ic]);
- float fab = 2*fa*fb, fac = 2*fa*fc, fbc = 2*fb*fc; fa *= fa; fb *= fb; fc *= fc;
- int isample, iaa = (ia+1)*(ia+2)/2-1, ibb = (ib+1)*(ib+2)/2-1, icc = (ic+1)*(ic+2)/2-1;
- int iab = iaa - ia + ib, iac = iaa - ia + ic, ibc = ibb - ib + ic;
- for (isample=0; isample<ma->n; isample++)
+ if ( qsum[ib]==0 ) continue;
+ int ibb = (ib+1)*(ib+2)/2-1;
+ int iab = iaa - ia + ib;
+ for (ic=0; ic<ib; ic++)
{
- if ( b->ploidy && 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; isample<ma->n; isample++)
+ {
+ 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 );
+ 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; }
+ 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 ( 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; }
- lk_sum = lk_tot>lk_sum ? lk_tot + log(1+exp(lk_sum-lk_tot)) : lk_sum + log(1+exp(lk_tot-lk_sum));
}
}
+ lk_sums[2] = lk_sum;
}
+ // 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++;
for (i=0; i<nals; i++) if ( max_als2&1<<i) n2++;
if ( n2<n1 && kf_gammap(1,2.0*(max_lk-max_lk2))<threshold )
{
- max_lk = max_lk2;
+ // the threshold not exceeded, use the second most likely set with fewer alleles
+ max_lk = max_lk2;
max_als = max_als2;
+ n1 = n2;
}
+ lk_sum = lk_sums[n1-1];
- // Get the BCF record ready for output
+ // Get the BCF record ready for GT and GQ
kstring_t s;
- 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);
- kputs(b->info, &s); if (b->info[0]) kputc(';', &s); 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)));
- if ( b->qual>999 ) b->qual = 999;
- //bcf_sync(b);
-
- int x, old_n_gi = b->n_gi;
+ int old_n_gi = b->n_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 GT
- int isample, gts=0;
+ // Call GTs
+ int isample, gts=0, ac[4] = {0,0,0,0};
+ int nRR = 0, nAA = 0, nRA = 0, max_dv = 0, dp_nref = 0;
for (isample = 0; isample < b->n_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_s = 0;
for (ia=0; ia<nals; ia++)
{
if ( !(max_als&1<<ia) ) continue;
int iaa = (ia+1)*(ia+2)/2-1;
- float _lk = p[iaa]*qsum[ia]*qsum[ia];
+ double _lk = p[iaa]*qsum[ia]*qsum[ia];
if ( _lk > lk ) { lk = _lk; als = ia<<3 | ia; }
- lk_sum += _lk;
+ lk_s += _lk;
}
if ( ploidy==2 )
{
{
if ( !(max_als&1<<ib) ) continue;
int iab = iaa - ia + ib;
- float _lk = 2*qsum[ia]*qsum[ib]*p[iab];
+ double _lk = 2*qsum[ia]*qsum[ib]*p[iab];
if ( _lk > lk ) { lk = _lk; als = ib<<3 | ia; }
- lk_sum += _lk;
+ lk_s += _lk;
}
}
}
- lk = -log(1-lk/lk_sum)/0.2302585;
- if ( idp>=0 && ((uint16_t*)b->gi[idp].data)[isample]==0 )
+ lk = -log(1-lk/lk_s)/0.2302585;
+ int dp = 0;
+ if ( idp>=0 && (dp=((uint16_t*)b->gi[idp].data)[isample])==0 )
{
- als |= 1<<7;
- lk = 0;
+ // no coverage
+ ((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;
+ if ( lk>99 ) lk = 99;
+ ((uint8_t*)b->gi[old_n_gi].data)[isample] = als;
+ ((uint8_t*)b->gi[old_n_gi+1].data)[isample] = (int)lk;
- gts |= (als>>3&7) | (als&7);
+ // For MDV annotation
+ int dv;
+ if ( als && idv>=0 && (dv=((uint16_t*)b->gi[idv].data)[isample]) )
+ {
+ if ( max_dv < dv ) max_dv = dv;
+ dp_nref += dp;
+ }
+
+ // For HWE annotation; multiple ALT alleles treated as one
+ if ( !als ) nRR++;
+ else if ( !(als>>3&7) || !(als&7) ) nRA++;
+ else nAA++;
+
+ gts |= 1<<(als>>3&7) | 1<<(als&7);
+ ac[ als>>3&7 ]++;
+ ac[ als&7 ]++;
}
+ free(pdg);
bcf_fit_alt(b,max_als);
+
+ // The VCF spec is ambiguous about QUAL: is it the probability of anything else
+ // (that is QUAL(non-ref) = P(ref)+P(any non-ref other than ALT)) or is it
+ // QUAL(non-ref)=P(ref) and QUAL(ref)=1-P(ref)? Assuming the latter.
+ b->qual = gts>1 ? -4.343*(ref_lk - lk_sum) : -4.343*log(1-exp(ref_lk - lk_sum));
+ if ( b->qual>999 ) b->qual = 999;
+
+ // 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; i<nals; i++)
+ {
+ an += ac[i];
+ if ( i>0 && ac[i] ) nalts++;
+ }
+ ksprintf(&s, "AN=%d;", an);
+ if ( nalts )
+ {
+ kputs("AC=", &s);
+ for (i=1; i<nals; i++)
+ {
+ if ( !(gts&1<<i) ) continue;
+ nalts--;
+ ksprintf(&s,"%d", ac[i]);
+ if ( nalts>0 ) 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);
+ ksprintf(&s, ";QBD=%e", b->qual/(a.d[0] + a.d[1] + a.d[2] + a.d[3]));
+ if ( dp_nref ) ksprintf(&s, ";QBDNR=%e", b->qual/dp_nref);
+ if ( max_dv ) ksprintf(&s, ";MDV=%d", max_dv);
+ }
+ if ( nAA+nRA )
+ {
+ double hwe = calc_hwe(nAA, nRR, nRA);
+ ksprintf(&s, ";HWE=%e", hwe);
+ }
+ 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;
bcf_sync(b);
- free(pdg);
return gts;
}
}
}
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)
projects / samtools.git / blobdiff