X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fcall1.c;h=cb33e06ef05438fe302ae398ad8de9dbd9ec3130;hb=50706bea83d8a1518485283876333279f4ae7137;hp=7aca1595cf7de74147a98492c08e34bd560c162f;hpb=e6ca3f825fea0a18b6068905e3173ff19104e756;p=samtools.git

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 7aca159..cb33e06 100644
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -38,6 +38,23 @@ void *bed_read(const char *fn);
 void bed_destroy(void *_h);
 int bed_overlap(const void *_h, const char *chr, int beg, int end);
 
+static double test_hwe(const double g[3])
+{
+	extern double kf_gammaq(double p, double x);
+	double fexp, chi2, f[3], n;
+	int i;
+	n = g[0] + g[1] + g[2];
+	fexp = (2. * g[2] + g[1]) / (2. * n);
+	if (fexp > 1. - 1e-10) fexp = 1. - 1e-10;
+	if (fexp < 1e-10) fexp = 1e-10;
+	f[0] = n * (1. - fexp) * (1. - fexp);
+	f[1] = n * 2. * fexp * (1. - fexp);
+	f[2] = n * fexp * fexp;
+	for (i = 0, chi2 = 0.; i < 3; ++i)
+		chi2 += (g[i] - f[i]) * (g[i] - f[i]) / f[i];
+	return kf_gammaq(.5, chi2 / 2.);
+}
+
 typedef struct {
 	double p[4];
 	int mq, depth, is_tested, d[4];
@@ -117,7 +134,11 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
 	kputs(b->info, &s);
 	if (b->info[0]) kputc(';', &s);
 //	ksprintf(&s, "AF1=%.4lg;AFE=%.4lg;CI95=%.4lg,%.4lg", 1.-pr->f_em, 1.-pr->f_exp, pr->cil, pr->cih);
-	ksprintf(&s, "AF1=%.4g;CI95=%.4g,%.4g", 1.-pr->f_em, pr->cil, pr->cih);
+	ksprintf(&s, "AF1=%.4g;CI95=%.4g,%.4g;G3=%.4g,%.4g,%.4g", 1.-pr->f_em, pr->cil, pr->cih, pr->g[2], pr->g[1], pr->g[0]);
+	if (n_smpl > 5) {
+		double hwe = test_hwe(pr->g);
+		if (hwe < 0.1) ksprintf(&s, ";HWE=%.4g", hwe);
+	}
 	if (has_I16) ksprintf(&s, ";DP4=%d,%d,%d,%d;MQ=%d", a.d[0], a.d[1], a.d[2], a.d[3], a.mq);
 	fq = pr->p_ref_folded < 0.5? -4.343 * log(pr->p_ref_folded) : 4.343 * log(pr->p_var_folded);
 	if (fq < -999) fq = -999;
@@ -216,6 +237,10 @@ static void write_header(bcf_hdr_t *h)
 		kputs("##INFO=<ID=FQ,Number=1,Type=Float,Description=\"Phred probability of all samples being the same\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=AF1,"))
 		kputs("##INFO=<ID=AF1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the site allele frequency of the first ALT allele\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=G3,"))
+		kputs("##INFO=<ID=G3,Number=3,Type=Float,Description=\"ML estimate of genotype frequencies\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=HWE,"))
+		kputs("##INFO=<ID=HWE,Number=1,Type=Float,Description=\"Chi^2 based HWE test P-value based on G3\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=CI95,"))
 		kputs("##INFO=<ID=CI95,Number=2,Type=Float,Description=\"Equal-tail Bayesian credible interval of the site allele frequency at the 95% level\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=PV4,"))
@@ -454,7 +479,7 @@ int bcfview(int argc, char *argv[])
 				kstring_t s;
 				s.m = s.l = 0; s.s = 0;
 				if (*b->info) kputc(';', &s);
-				ksprintf(&s, "NEIR=%.3f;NEIF=%.3f,%.3f", r2, f[0]+f[2], f[0]+f[1]);
+				ksprintf(&s, "NEIR=%.3f;NEIF4=%.3f,%.3f,%.3f,%.3f", r2, f[0], f[1], f[2], f[3]);
 				bcf_append_info(b, s.s, s.l);
 				free(s.s);
 			}