Skip to content

Commit

Permalink
call: improve QUAL calculation (samtools#449)
Browse files Browse the repository at this point in the history 
- don't lose precision by summing 1+\epsilon for very small \epsilon
- never output QUAL=0
- for the test case, sites with QUAL=0 previously now have QUAL>0
  2 sites that were previously QUAL=999 now have lower QUAL. In
  constrained calling mode `-C` these sites that should be called
  non-ref with QUAL 999 are being forced into a ref call, so the
  QUAL value should be lowered in this case

Resolves
    samtools#370 (comment)
    samtools/samtools#566
  • Loading branch information
@mcshane
mcshane authored Jul 22, 2016
1 parent 74733e6 commit 66b6e76
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 14 deletions.
16 changes: 9 additions & 7 deletions mcall.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -572,9 +572,9 @@ static inline double logsumexp2(double a, double b)
}

// Macro to set the most likely alleles
#define UPDATE_MAX_LKs(als) { \
#define UPDATE_MAX_LKs(als,sum) { \
if ( max_lk<lk_tot ) { max_lk = lk_tot; max_als = (als); } \
if ( lk_tot_set ) lk_sum = logsumexp2(lk_tot,lk_sum); \
if ( sum ) lk_sum = logsumexp2(lk_tot,lk_sum); \
}

#define SWAP(type_t,x,y) {type_t tmp; tmp = x; x = y; y = tmp; }
Expand Down Expand Up @@ -603,9 +603,10 @@ static int mcall_find_best_alleles(call_t *call, int nals, int *out_als)
if ( *pdg ) { lk_tot += log(*pdg); lk_tot_set = 1; }
pdg += ngts;
}

if ( ia==0 ) ref_lk = lk_tot; // likelihood of 0/0 for all samples
else lk_tot += call->theta; // the prior
UPDATE_MAX_LKs(1<<ia);
UPDATE_MAX_LKs(1<<ia, ia>0 && lk_tot_set);
}

// Two alleles
Expand Down Expand Up @@ -637,7 +638,7 @@ static int mcall_find_best_alleles(call_t *call, int nals, int *out_als)
}
if ( ia!=0 ) lk_tot += call->theta; // the prior
if ( ib!=0 ) lk_tot += call->theta;
UPDATE_MAX_LKs(1<<ia|1<<ib);
UPDATE_MAX_LKs(1<<ia|1<<ib, lk_tot_set);
}
}
}
Expand Down Expand Up @@ -679,7 +680,7 @@ static int mcall_find_best_alleles(call_t *call, int nals, int *out_als)
if ( ia!=0 ) lk_tot += call->theta; // the prior
if ( ib!=0 ) lk_tot += call->theta; // the prior
if ( ic!=0 ) lk_tot += call->theta; // the prior
UPDATE_MAX_LKs(1<<ia|1<<ib|1<<ic);
UPDATE_MAX_LKs(1<<ia|1<<ib|1<<ic, lk_tot_set);
}
}
}
Expand Down Expand Up @@ -1528,13 +1529,14 @@ int mcall(call_t *call, bcf1_t *rec)
if ( hob != HUGE_VAL ) bcf_update_info_float(call->hdr, rec, "HOB", &hob, 1);

// Quality of a variant site. fabs() to avoid negative zeros in VCF output when CALL_KEEPALT is set
rec->qual = call->lk_sum==-HUGE_VAL || call->ref_lk==0 ? 0 : fabs(-4.343*(call->ref_lk - call->lk_sum));
rec->qual = -4.343*(call->ref_lk - logsumexp2(call->lk_sum,call->ref_lk));
}
else
{
// Set the quality of a REF site
rec->qual = call->lk_sum==-HUGE_VAL || call->ref_lk==0 ? 0 : -4.343*log(1 - exp(call->ref_lk - call->lk_sum));
rec->qual = -4.343*(call->lk_sum - logsumexp2(call->lk_sum,call->ref_lk));
}

if ( rec->qual>999 ) rec->qual = 999;
if ( rec->qual>50 ) rec->qual = rint(rec->qual);

Expand Down
14 changes: 7 additions & 7 deletions test/mpileup.cAls.out
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,14 +27,14 @@
##INFO=<ID=DP4,Number=4,Type=Integer,Description="Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases">
##INFO=<ID=MQ,Number=1,Type=Integer,Description="Average mapping quality">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102
17 1 . A G,T 0 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 2 . A T,G 0 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 3 . A C 0 . DP=11;MQ0F=0;AC=0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 1 . A G,T 52 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 2 . A T,G 52 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 3 . A C 999 . DP=11;MQ0F=0;AC=0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 4 . A G,T,C 21.815 . DP=11;MQ0F=0;AC=0,0,0;AN=2;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV 0/0:1,2,3,7,8,10,11,12,14,15:5:0 ./.:.:3:0 ./.:.:3:0
17 5 . A G,T 0 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 6 . A T,G 0 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 5 . A G,T 999 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 6 . A T,G 999 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0
17 7 . A T,G,C 21.5769 . DP=11;MQ0F=0;AC=0,0,0;AN=2;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV 0/0:1,2,3,4,5,6,2,3,5,3:5:0 ./.:.:3:0 ./.:.:3:0
17 828 . T C 409 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4
17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2
17 2220 . G C 999 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=0;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/0:139,157,255:12:6 0/0:69,75,119:4:2 0/0:131,131,131:4:4
17 2564 . A AG 999 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=0;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/0:88,98,171:6:3 0/0:57,63,117:4:2 0/0:124,124,124:4:4
17 2220 . G C 189 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=0;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/0:139,157,255:12:6 0/0:69,75,119:4:2 0/0:131,131,131:4:4
17 2564 . A AG 166 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=0;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/0:88,98,171:6:3 0/0:57,63,117:4:2 0/0:124,124,124:4:4

0 comments on commit 66b6e76

Please sign in to comment.