Fix samtools#535

call.h

@@ -103,7 +103,7 @@ call_t;
 void error(const char *format, ...);
 
 /*
- *  *call() - return negative value on error or the number of non-reference
+ *  call() - return -1 value on critical error; -2 to skip the site; or the number of non-reference
  *            alleles on success.
  */
 int mcall(call_t *call, bcf1_t *rec);    // multiallic and rare-variant calling model

mcall.c

@@ -1270,7 +1270,7 @@ void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int o
 // 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)
+static int 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);
@@ -1293,7 +1293,7 @@ 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+1==*unseen ) { fprintf(stderr,"fixme? Cannot constrain to %s\n",tgt->als[i]); return -1; }
 
         if ( j>=0 )
         {
@@ -1319,7 +1319,7 @@ static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int *unseen)
         nals++;
     }
 
-    if ( !has_new && nals==rec->n_allele ) return;
+    if ( !has_new && nals==rec->n_allele ) return 0;
     bcf_update_alleles(call->hdr, rec, (const char**)call->als, nals);
 
     // create mapping from new PL to old PL
@@ -1371,6 +1371,7 @@ static void mcall_constrain_alleles(call_t *call, bcf1_t *rec, int *unseen)
     bcf_update_info_float(call->hdr, rec, "QS", qsum, nals);
 
     if ( *unseen ) *unseen = nals-1;
+    return 0;
 }
 
 
@@ -1385,7 +1386,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)!=0 ) return -2;
 
     int nsmpl = bcf_hdr_nsamples(call->hdr);
     int nals  = rec->n_allele;

vcfcall.c

@@ -852,6 +852,7 @@ int main_vcfcall(int argc, char *argv[])
         else
             ret = ccall(&args.aux, bcf_rec);
         if ( ret==-1 ) error("Something is wrong\n");
+        else if ( ret==-2 ) continue;   // skip the site
 
         // Normal output
         if ( (args.aux.flag & CALL_VARONLY) && ret==0 && !args.gvcf ) continue;     // not a variant