Skip to content

Commit

Permalink
call -mC alleles fix: Do not take X as last in ALT list for granted
Browse files Browse the repository at this point in the history
  • Loading branch information
@pd3
pd3 committed Nov 30, 2015
1 parent f0aff23 commit 47ab239
Show file tree
Hide file tree
Showing 4 changed files with 43 additions and 14 deletions.
29 changes: 22 additions & 7 deletions mcall.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1235,7 +1235,11 @@ void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int o
}
}

static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int unseen)

// NB: in this function we temporarily use calls->als_map for a different
// purpose to store mapping from new (target) alleles to original alleles.
//
static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int *unseen)
{
bcf_sr_regions_t *tgt = call->srs->targets;
if ( tgt->nals>5 ) error("Maximum accepted number of alleles is 5, got %d\n", tgt->nals);
Expand All @@ -1258,6 +1262,8 @@ static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int unseen)
call->als[nals] = tgt->als[i];
j = vcmp_find_allele(call->vcmp, rec->d.allele+1, rec->n_allele - 1, tgt->als[i]);

if ( j+1==*unseen ) error("Cannot constrain to %s\n",tgt->als[i]);

if ( j>=0 )
{
// existing allele
Expand All @@ -1270,11 +1276,18 @@ static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int unseen)
// present at multiallelic indels sites. In that case we use the
// last allele anyway, because the least likely allele comes last
// in mpileup's ALT output.
call->als_map[nals] = unseen>=0 ? unseen : rec->n_allele - 1;
call->als_map[nals] = (*unseen)>=0 ? *unseen : rec->n_allele - 1;
has_new = 1;
}
nals++;
}
if ( *unseen )
{
call->als_map[nals] = *unseen;
call->als[nals] = rec->d.allele[*unseen];
nals++;
}

if ( !has_new && nals==rec->n_allele ) return;
bcf_update_alleles(call->hdr, rec, (const char**)call->als, nals);

Expand All @@ -1301,15 +1314,15 @@ static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int unseen)
for (k=0; k<npls_new; k++)
{
new_pl[k] = ori_pl[call->pl_map[k]];
if ( new_pl[k]==bcf_int32_missing && unseen>=0 )
if ( new_pl[k]==bcf_int32_missing && *unseen>=0 )
{
// missing value, and there is an unseen allele: identify the
// alleles and use the lk of either AX or XX
int k_ori = call->pl_map[k], ia, ib;
bcf_gt2alleles(k_ori, &ia, &ib);
k_ori = bcf_alleles2gt(ia,unseen);
if ( ori_pl[k_ori]==bcf_int32_missing ) k_ori = bcf_alleles2gt(ib,unseen);
if ( ori_pl[k_ori]==bcf_int32_missing ) k_ori = bcf_alleles2gt(unseen,unseen);
k_ori = bcf_alleles2gt(ia,*unseen);
if ( ori_pl[k_ori]==bcf_int32_missing ) k_ori = bcf_alleles2gt(ib,*unseen);
if ( ori_pl[k_ori]==bcf_int32_missing ) k_ori = bcf_alleles2gt(*unseen,*unseen);
new_pl[k] = ori_pl[k_ori];
}
if ( !k && new_pl[k]==bcf_int32_vector_end ) new_pl[k]=bcf_int32_missing;
Expand All @@ -1325,6 +1338,8 @@ static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int unseen)
for (i=0; i<nals; i++)
qsum[i] = call->als_map[i]<nqs ? call->qsum[call->als_map[i]] : 0;
bcf_update_info_float(call->hdr, rec, "QS", qsum, nals);

if ( *unseen ) *unseen = nals-1;
}


Expand All @@ -1339,7 +1354,7 @@ int mcall(call_t *call, bcf1_t *rec)
int i, unseen = call->unseen;

// Force alleles when calling genotypes given alleles was requested
if ( call->flag & CALL_CONSTR_ALLELES ) mcall_constrain_alleles(call, rec, unseen);
if ( call->flag & CALL_CONSTR_ALLELES ) mcall_constrain_alleles(call, rec, &unseen);

int nsmpl = bcf_hdr_nsamples(call->hdr);
int nals = rec->n_allele;
Expand Down
7 changes: 7 additions & 0 deletions test/mpileup.cAls.out
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,13 @@
##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 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 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
Expand Down
7 changes: 7 additions & 0 deletions test/mpileup.tab
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,10 @@
17 1 A,G,T
17 2 A,T,G
17 3 A,C
17 4 A,C,T,G
17 5 A,G,T
17 6 A,T,G
17 7 A,T,G,C
17 828 T,C
17 1665 T,C
17 2220 G,C
Expand Down
14 changes: 7 additions & 7 deletions test/mpileup.vcf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,13 @@
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Number of high-quality bases">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="Number of high-quality non-reference bases">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102
17 1 . A <X> 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 2 . A <X> 0 . DP=11;I16=11,0,0,0,439,17587,0,0,319,9251,0,0,226,5030,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 3 . G <X> 0 . DP=11;I16=11,0,0,0,431,16971,0,0,319,9251,0,0,229,5111,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 4 . C <X> 0 . DP=11;I16=11,0,0,0,423,16417,0,0,319,9251,0,0,232,5202,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,71:3:0
17 5 . T <X> 0 . DP=11;I16=11,0,0,0,450,18520,0,0,319,9251,0,0,234,5252,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 6 . T <X> 0 . DP=11;I16=11,0,0,0,403,14847,0,0,319,9251,0,0,236,5310,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 7 . C <X> 0 . DP=11;I16=11,0,0,0,446,18114,0,0,319,9251,0,0,237,5327,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 1 . A G,X,T,C 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,1,0,1,1;MQ0F=0 PL:DP:DV 0,0,0,0,0,0,.,.,.,.,.,.,.,.,.:5:0 .:3:0 .:3:0
17 2 . A G,X,T,C 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,1,0,1,1;MQ0F=0 PL:DP:DV 0,0,0,0,0,0,.,.,.,.,.,.,.,.,.:5:0 .:3:0 .:3:0
17 3 . A G,X,T,C 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,1,0,1,1;MQ0F=0 PL:DP:DV 0,0,0,0,0,0,.,.,.,.,.,.,.,.,.:5:0 .:3:0 .:3:0
17 4 . A G,X,T,C 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,1,0,1,1;MQ0F=0 PL:DP:DV 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15:5:0 .:3:0 .:3:0
17 5 . A X,G 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,0,1;MQ0F=0 PL:DP:DV 0,0,0,0,0,0:5:0 .:3:0 .:3:0
17 6 . A X,G 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,0,1;MQ0F=0 PL:DP:DV 0,0,0,0,0,0:5:0 .:3:0 .:3:0
17 7 . A X,G 0 . DP=11;I16=11,0,0,0,452,18594,0,0,319,9251,0,0,223,4959,0,0;QS=1,0,1;MQ0F=0 PL:DP:DV 1,2,3,4,5,6:5:0 .:3:0 .:3:0
17 8 . T <X> 0 . DP=11;I16=11,0,0,0,465,19677,0,0,319,9251,0,0,238,5354,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 9 . C <X> 0 . DP=11;I16=11,0,0,0,447,18205,0,0,319,9251,0,0,239,5391,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,72:3:0
17 10 . A <X> 0 . DP=11;I16=11,0,0,0,426,16756,0,0,319,9251,0,0,240,5438,0,0;QS=3,0;MQ0F=0 PL:DP:DV 0,15,100:5:0 0,9,72:3:0 0,9,69:3:0
Expand Down

0 comments on commit 47ab239

Please sign in to comment.