From 213027e3cf754de8f601d3155899119f952fe617 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Mon, 28 Feb 2011 20:24:04 +0000
Subject: [PATCH]  * 0.1.12-r924:126  * fixed a bug in phase (due to recent
 changes)  * fixed a bug in vcf2fq

---
 bam.h                |  2 +-
 bcftools/vcfutils.pl |  4 ++--
 phase.c              | 12 +++++++-----
 3 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/bam.h b/bam.h
index 752d0bd..41f8d95 100644
--- a/bam.h
+++ b/bam.h
@@ -40,7 +40,7 @@
   @copyright Genome Research Ltd.
  */
 
-#define BAM_VERSION "0.1.12-r921:124"
+#define BAM_VERSION "0.1.12-r924:126"
 
 #include <stdint.h>
 #include <stdlib.h>
diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 12807c6..6ced66a 100755
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -201,7 +201,7 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 	  my @x;
 	  $x[0] = sprintf("%.4f", $c[1]-$c[2]? $c[2] / ($c[1] - $c[2]) : 100);
 	  $x[1] = sprintf("%.4f", $hsize? $c[3] / $hsize : 0);
-	  $x[2] = sprintf("%.4f", $c[1] ? $c[3] / $c[1] : 0);
+	  $x[2] = sprintf("%.4f", $c[3] / $c[1]);
 	  my $a = $c[1] - $lc[1];
 	  my $b = $c[2] - $lc[2];
 	  $x[3] = sprintf("%.4f", $a-$b? $b / ($a-$b) : 100);
@@ -512,7 +512,7 @@ Options: -d INT    minimum depth          [$opts{d}]
 	  $q = chr($q <= 126? $q : 126);
 	  $seq .= $b;
 	  $qual .= $q;
-	} else { # an INDEL
+	} elsif ($t[4] ne '.') { # an INDEL
 	  push(@gaps, [$t[1], length($t[3])]);
 	}
 	$last_pos = $t[1];
diff --git a/phase.c b/phase.c
index 5fc4ce7..aa88921 100644
--- a/phase.c
+++ b/phase.c
@@ -34,7 +34,7 @@ typedef struct {
 typedef struct {
 	int8_t seq[MAX_VARS]; // TODO: change to dynamic memory allocation!
 	int vpos, beg, end;
-	uint32_t vlen:16, single:1, flip:1, phase:1, phased:1;
+	uint32_t vlen:16, single:1, flip:1, phase:1, phased:1, ambig:1;
 	uint32_t in:16, out:16; // in-phase and out-phase
 } frag_t, *frag_p;
 
@@ -176,7 +176,8 @@ static uint64_t *fragphase(int vpos, const int8_t *path, nseq_t *hash, int flip)
 			}
 			f->phase = c[0] > c[1]? 0 : 1;
 			f->in = c[f->phase]; f->out = c[1 - f->phase];
-			if (f->in && f->out && f->in <= f->out + 1) f->phased = 0;
+			f->phased = f->in == f->out? 0 : 1;
+			f->ambig = (f->in && f->out && f->in <= f->out + 1)? 1 : 0;
 			// fix chimera
 			f->flip = 0;
 			if (flip && c[0] >= 3 && c[1] >= 3) {
@@ -321,8 +322,9 @@ static void dump_aln(phaseg_t *g, int min_pos, const nseq_t *hash)
 		if (k == kh_end(hash)) which = 3;
 		else {
 			frag_t *f = &kh_val(hash, k);
-			if (f->phased && f->flip) which = 2;
-			else if (f->phased == 0) which = 2;
+			if (f->ambig) which = 2;
+			else if (f->phased && f->flip) which = 2;
+			else if (f->phased == 0) which = 3;
 			else { // phased and not flipped
 				char c = 'Y';
 				which = f->phase;
@@ -652,7 +654,7 @@ int main_phase(int argc, char *argv[])
 				memset(f->seq, 0, MAX_VARS);
 				f->beg = p->b->core.pos;
 				f->end = bam_calend(&p->b->core, bam1_cigar(p->b));
-				f->vpos = vpos, f->vlen = 1, f->seq[0] = c, f->single = f->phased = f->flip = 0;
+				f->vpos = vpos, f->vlen = 1, f->seq[0] = c, f->single = f->phased = f->flip = f->ambig = 0;
 			}
 		}
 		if (dophase) {
-- 
2.39.2