Make FORMAT/AD work with call -C alleles

Fixes samtools#887

Author: pd3

mcall.c

@@ -1370,6 +1370,47 @@ static int mcall_constrain_alleles(call_t *call, bcf1_t *rec, int *unseen)
         qsum[i] = call->als_map[i]<nqs ? call->qsum[call->als_map[i]] : 0;
     bcf_update_info_float(call->hdr, rec, "QS", qsum, nals);
 
+    // update any Number=R tags
+    void *tmp_ori = call->itmp, *tmp_new = call->PLs;  // reusing PLs storage which is not used at this point
+    int ntmp_ori = call->n_itmp, ntmp_new = call->mPLs;
+    for (i=0; i<rec->n_fmt; i++)
+    {
+        bcf_fmt_t *fmt = &rec->d.fmt[i];
+        int vlen = bcf_hdr_id2length(call->hdr,BCF_HL_FMT,fmt->id);
+        if ( vlen!=BCF_VL_R ) continue; // not a Number=R tag
+
+        // NB:works only for BCF_HT_INT and BCF_HT_REAL
+        int type = bcf_hdr_id2type(call->hdr,BCF_HL_FMT,fmt->id);
+        assert( type==BCF_HT_INT || type==BCF_HT_REAL );
+        assert( sizeof(float)==sizeof(int32_t) );
+
+        const char *key = bcf_hdr_int2id(call->hdr,BCF_DT_ID,fmt->id);
+        int nret = bcf_get_format_values(call->hdr, rec, key, &tmp_ori, &ntmp_ori, type);
+        if (nret<=0) continue;
+        int nsmpl = bcf_hdr_nsamples(call->hdr);
+        int size1 = sizeof(float);
+
+        for (j=0; j<nsmpl; j++)
+        {
+            void *ptr_ori = tmp_ori + j*size1*fmt->n;
+            void *ptr_new = tmp_new + j*nals*size1;
+            for (k=0; k<nals; k++)
+            {
+                void *dst = ptr_new + size1*k;
+                void *src = ptr_ori + size1*call->als_map[k]; 
+                memcpy(dst,src,size1);
+            }
+        }
+        ntmp_new = nsmpl*rec->n_allele;
+        nret = bcf_update_format(call->hdr, rec, key, tmp_new, ntmp_new, type);
+        assert( nret==0 );
+    }
+    call->PLs    = (int32_t*) tmp_new;
+    call->mPLs   = ntmp_new;
+    call->itmp   = (int32_t*) tmp_ori;
+    call->n_itmp = ntmp_ori;
+
+
     if ( *unseen ) *unseen = nals-1;
     return 0;
 }

test/mpileup.2.tab

@@ -0,0 +1,6 @@
+chr1	212740	A,G
+chr1	320055	A,G
+chr1	486173	A,T
+chr1	511277	A,G
+chr1	602567	A,G
+chr1	639707	T,A

test/mpileup.2.vcf

@@ -0,0 +1,28 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##reference=file://test/mpileup.ref.fa.gz
+##contig=<ID=chr1,length=81195210>
+##ALT=<ID=X,Description="Represents allele(s) other than observed.">
+##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
+##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">
+##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">
+##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
+##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
+##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">
+##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">
+##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">
+##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">
+##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
+##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
+##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
+##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
+##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Number of high-quality bases">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample1	sample2
+chr1	212740	.	A	G,<*>	0	.	DP=73;I16=0,0,39,4,0,0,2743,192525,0,0,2580,154800,0,0,825,18621;QS=0,2,0;VDB=0.520868;SGB=-1.38232;MQSB=1;MQ0F=0	PL:DP:AD	255,72,0,255,72,255:24:0,24,0	255,57,0,255,57,255:19:0,19,0
+chr1	320055	.	A	<*>	0	.	DP=101;I16=52,9,0,0,4116,300666,0,0,3660,219600,0,0,1281,29849,0,0;QS=2,0;MQSB=1;MQ0F=0	PL:DP:AD	0,87,255:29:29,0	0,96,255:32:32,0
+chr1	486173	.	A	T,<*>	0	.	DP=13;I16=3,1,3,0,287,21853,246,20172,240,14400,180,10800,95,2275,75,1875;QS=1.25,0.75,0;VDB=0.074936;SGB=0.620439;RPB=0.810265;MQB=1.01283;MQSB=1;BQB=0.810265;MQ0F=0	PL:DP:AD	0,9,151,9,151,151:3:3,0,0	140,0,48,143,57,194:4:1,3,0
+chr1	511277	.	A	G,<*>	0	.	DP=50;I16=0,0,25,4,0,0,1900,137374,0,0,1740,104400,0,0,672,16306;QS=0,2,0;VDB=0.0722735;SGB=-1.26186;MQSB=1;MQ0F=0	PL:DP:AD	255,30,0,255,30,255:10:0,10,0	255,57,0,255,57,255:19:0,19,0
+chr1	602567	.	A	G,<*>	0	.	DP=9;I16=3,1,1,0,328,26896,41,1681,240,14400,60,3600,99,2451,18,324;QS=1.81448,0.18552,0;SGB=-0.516033;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0	PL:DP:AD	0,3,60,3,60,60:1:1,0,0	29,0,140,38,143,175:4:3,1,0
+chr1	639707	.	T	A,<*>	0	.	DP=50;I16=0,0,23,8,0,0,1998,142356,0,0,1860,111600,0,0,612,13818;QS=0,2,0;VDB=0.563111;SGB=-1.37269;MQSB=1;MQ0F=0	PL:DP:AD	255,42,0,255,42,255:14:0,14,0	255,51,0,255,51,255:17:0,17,0

test/test.pl

@@ -203,6 +203,7 @@
 test_vcf_call($opts,in=>'mpileup.X',out=>'mpileup.X.out',args=>'-mv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.ped');
 test_vcf_call($opts,in=>'mpileup.X',out=>'mpileup.X.2.out',args=>'-mv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.2.samples');
 test_vcf_call_cAls($opts,in=>'mpileup',out=>'mpileup.cAls.out',tab=>'mpileup');
+test_vcf_call_cAls($opts,in=>'mpileup.2',out=>'mpileup.cAls.2.out',tab=>'mpileup.2');
 test_vcf_call($opts,in=>'mpileup.c',out=>'mpileup.c.1.out',args=>'-cv');
 # test_vcf_call($opts,in=>'mpileup.c',out=>'mpileup.c.2.out',args=>'-cg0');
 test_vcf_call($opts,in=>'mpileup.c.X',out=>'mpileup.c.X.out',args=>'-cv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.samples');