DPR trimmed with call -c

overerd · Oct 8, 2014 · f53ec1a · f53ec1a
1 parent 15c718e
commit f53ec1a
Show file tree

Hide file tree

Showing 5 changed files with 21 additions and 3 deletions.
diff --git a/Makefile b/Makefile
@@ -105,7 +105,7 @@ PLUGINP = $(PLUGINC:.c=.P)
 	echo -e "\t\\" >> $@; \
 	echo '$(CC) $(CFLAGS) $(INCLUDES) -fPIC -shared -o $$@ $$< version.c \' >> $@; \
 	echo '`if [ -e $*.dep ]; then cat $*.dep; fi` \' >> $@; \
-	echo '-L$(HTSDIR) -lhts' >> $@;
+	echo '-L$(HTSDIR) -lhts' >> $@; cat $@
 
 -include $(PLUGINP)
 

diff --git a/call.h b/call.h
@@ -132,4 +132,7 @@ uint32_t *call_trio_prep(int is_x, int is_son);
 /** gVCF */
 void gvcf_write(htsFile *fh, gvcf_t *gvcf, bcf_hdr_t *hdr, bcf1_t *rec, int is_ref);
 
+void init_allele_trimming_maps(call_t *call, int als, int nals);
+void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als);
+
 #endif
diff --git a/ccall.c b/ccall.c
@@ -43,6 +43,8 @@ void ccall_init(call_t *call)
     call_init_pl2p(call);
     call->cdat->p1 = bcf_p1_init(bcf_hdr_nsamples(call->hdr), call->ploidy);
     call->gts = (int*) calloc(bcf_hdr_nsamples(call->hdr)*2,sizeof(int));   // assuming at most diploid everywhere
+    call->nals_map = 5;
+    call->als_map  = (int*) malloc(sizeof(int)*call->nals_map);
 
     bcf_hdr_append(call->hdr,"##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">");
     bcf_hdr_append(call->hdr,"##INFO=<ID=AF1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the first ALT allele frequency (assuming HWE)\">");
@@ -62,6 +64,8 @@ void ccall_init(call_t *call)
 }
 void ccall_destroy(call_t *call)
 {
+    free(call->itmp);
+    free(call->als_map);
     free(call->gts);
     free(call->anno16);
     free(call->PLs);
@@ -219,7 +223,7 @@ static int update_bcf1(call_t *call, bcf1_t *rec, const bcf_p1rst_t *pr, double
     if (rec->qual > 999) rec->qual = 999;
 
     // Remove unused alleles
-    int nals = !is_var ? 1 : pr->rank0 < 2? 2 : pr->rank0+1;
+    int nals_ori = rec->n_allele, nals = !is_var ? 1 : pr->rank0 < 2? 2 : pr->rank0+1;
     if ( nals<rec->n_allele )
     {
         bcf_update_alleles(call->hdr, rec, (const char**)rec->d.allele, nals);
@@ -283,6 +287,13 @@ static int update_bcf1(call_t *call, bcf1_t *rec, const bcf_p1rst_t *pr, double
         // GQ: todo
     }
     bcf_update_genotypes(call->hdr, rec, call->gts, rec->n_sample*2);
+
+    // trim Number=R tags
+    int out_als = 0;
+    for (i=0; i<nals; i++) out_als |= 1<<i;
+    init_allele_trimming_maps(call, out_als, nals_ori);
+    mcall_trim_numberR(call, rec, nals_ori, nals, out_als);
+
     return is_var;
 }
 

diff --git a/mcall.c b/mcall.c
@@ -468,6 +468,8 @@ void init_allele_trimming_maps(call_t *call, int als, int nals)
         else call->als_map[i] = -1;
     }
 
+    if ( !call->pl_map ) return;
+
     // pl_map: new(k) -> old(l)
     int k = 0, l = 0;
     for (i=0; i<nals; i++)
@@ -1178,7 +1180,7 @@ static void mcall_trim_PLs(call_t *call, bcf1_t *rec, int nals, int nout_als, in
     bcf_update_format_int32(call->hdr, rec, "PL", call->PLs, npls_dst*nsmpl);
 }
 
-static void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als)
+void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als)
 {
     int i, ret;
 

diff --git a/plugins/vcf2sex.c b/plugins/vcf2sex.c
@@ -113,6 +113,8 @@ int process_region_precise(args_t *args, char *seq, regitr_t *itr)
                 if ( gts[igt]==bcf_gt_missing || gts[igt]==bcf_int32_missing || gts[igt]==bcf_int32_vector_end ) break;
                 else ploidy++;
             args->counts[ismpl*(args->max_ploidy+1) + ploidy]++;
+            if ( args->verbose )
+                fprintf(stderr,"%s:%d\t%s\tploidy=%d\n", seq,rec->pos+1,args->hdr->samples[ismpl],ploidy);
         }
     }