From 70c740facc966321754c6bfcc6d61ea056480638 Mon Sep 17 00:00:00 2001
From: Petr Danecek <pd3@sanger.ac.uk>
Date: Wed, 10 Oct 2012 12:47:41 +0100
Subject: [PATCH] Complain when BAM cannot be open. Severe bug fixed in -m
 haploid calling.

---
 bam_plcmd.c      |  5 +++++
 bcftools/prob1.c | 13 ++++++++-----
 2 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/bam_plcmd.c b/bam_plcmd.c
index ed00d2e..f77be04 100644
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -208,6 +208,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
 		bam_header_t *h_tmp;
 		data[i] = calloc(1, sizeof(mplp_aux_t));
 		data[i]->fp = strcmp(fn[i], "-") == 0? bam_dopen(fileno(stdin), "r") : bam_open(fn[i], "r");
+        if ( !data[i]->fp )
+        {
+            fprintf(stderr, "[%s] failed to open %d-th input.\n", __func__, i+1);
+            exit(1);
+        }
 		data[i]->conf = conf;
 		h_tmp = bam_header_read(data[i]->fp);
 		data[i]->h = i? h : h_tmp; // for i==0, "h" has not been set yet
diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index e3d6b5e..5051cc3 100644
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -300,10 +300,12 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
                 int isample, ibb = (ib+1)*(ib+2)/2-1, iab = iaa - ia + ib;
                 for (isample=0; isample<ma->n; isample++)
                 {
-                    if ( b->ploidy && b->ploidy[isample]==1 ) continue;
                     double *p = pdg + isample*npdg;
                     //assert( log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]) <= 0 );
-                    lk_tot += log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]);
+                    if ( b->ploidy && b->ploidy[isample]==1 )
+                        lk_tot +=  log(fa*p[iaa] + fb*p[ibb]);
+                    else 
+                        lk_tot +=  log(fa*p[iaa] + fb*p[ibb] + fab*p[iab]);
                 }
                 if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia|1<<ib; }
                 else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia|1<<ib; }
@@ -334,10 +336,12 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
                     int iac = iaa - ia + ic, ibc = ibb - ib + ic;
                     for (isample=0; isample<ma->n; isample++)
                     {
-                        if ( b->ploidy && b->ploidy[isample]==1 ) continue;
                         double *p = pdg + isample*npdg;
                         //assert( log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]) <= 0 );
-                        lk_tot += log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]);
+                        if ( b->ploidy && b->ploidy[isample]==1 ) 
+                            lk_tot += log(fa*p[iaa] + fb*p[ibb] + fc*p[icc]);
+                        else
+                            lk_tot += log(fa*p[iaa] + fb*p[ibb] + fc*p[icc] + fab*p[iab] + fac*p[iac] + fbc*p[ibc]);
                     }
                     if ( max_lk<lk_tot ) { max_lk2 = max_lk; max_als2 = max_als; max_lk = lk_tot; max_als = 1<<ia|1<<ib|1<<ic; }
                     else if ( max_lk2<lk_tot ) { max_lk2 = lk_tot; max_als2 = 1<<ia|1<<ib|1<<ic; }
@@ -347,7 +351,6 @@ int call_multiallelic_gt(bcf1_t *b, bcf_p1aux_t *ma, double threshold)
         }
     }
 
-
     // Should we add another allele, does it increase the likelihood significantly?
     int n1=0, n2=0;
     for (i=0; i<nals; i++) if ( max_als&1<<i) n1++;
-- 
2.39.2