From: Heng Li <lh3@live.co.uk>
Date: Thu, 2 Sep 2010 03:24:47 +0000 (+0000)
Subject: SNP calling from the GL field
X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;ds=sidebyside;h=018b5036f1db2b928ee92931fb2ca9b81116be71;p=samtools.git

SNP calling from the GL field
---

diff --git a/bcftools/bcf.h b/bcftools/bcf.h
index e6ba225..d0f9b02 100644
--- a/bcftools/bcf.h
+++ b/bcftools/bcf.h
@@ -66,6 +66,8 @@ extern "C" {
 	int vcf_read(bcf_t *bp, bcf_hdr_t *h, bcf1_t *b);
 
 	int bcf_shrink_alt(int n_smpl, bcf1_t *b, int n);
+	int bcf_gl2pl(int n_smpl, bcf1_t *b);
+
 	void *bcf_build_refhash(bcf_hdr_t *h);
 	void bcf_str2id_destroy(void *_hash);
 	int bcf_str2id_add(void *_hash, const char *str);
diff --git a/bcftools/bcfutils.c b/bcftools/bcfutils.c
index 6908bc4..dd1e89f 100644
--- a/bcftools/bcfutils.c
+++ b/bcftools/bcfutils.c
@@ -82,3 +82,29 @@ int bcf_shrink_alt(int n_smpl, bcf1_t *b, int n)
 	bcf_sync(n_smpl, b);
 	return 0;
 }
+
+int bcf_gl2pl(int n_smpl, bcf1_t *b)
+{
+	char *p;
+	int i;
+	bcf_ginfo_t *g;
+	float *d0;
+	uint8_t *d1;
+	if (strstr(b->fmt, "PL")) return -1;
+	if ((p = strstr(b->fmt, "GL")) == 0) return -1;
+	*p = 'P';
+	for (i = 0; i < b->n_gi; ++i)
+		if (b->gi[i].fmt == bcf_str2int("GL", 2))
+			break;
+	g = b->gi + i;
+	g->fmt = bcf_str2int("PL", 2);
+	g->len /= 4; // 4 == sizeof(float)
+	d0 = (float*)g->data; d1 = (uint8_t*)g->data;
+	for (i = 0; i < n_smpl * g->len; ++i) {
+		int x = (int)(-10. * d0[i] + .499);
+		if (x > 255) x = 255;
+		if (x < 0) x = 0;
+		d1[i] = x;
+	}
+	return 0;
+}
diff --git a/bcftools/call1.c b/bcftools/call1.c
index 1060513..995d500 100644
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -287,6 +287,7 @@ int bcfview(int argc, char *argv[])
 		}
 		if (vc.flag & VC_CALL) {
 			bcf_p1rst_t pr;
+			bcf_gl2pl(h->n_smpl, b);
 			bcf_p1_cal(b, p1, &pr); // pr.g[3] is not calculated here
 			if (vc.flag&VC_HWE) bcf_p1_cal_g3(p1, pr.g);
 			if ((n_processed + 1) % 50000 == 0) bcf_p1_dump_afs(p1);