X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fld.c;h=0207819bc731762b30908cef687b35134cee0f65;hb=674ffee7adcfc928f5039777180206fdebc5539b;hp=01145b84ef81538e9310a2cc61feb295dd9481d8;hpb=bd35dfe8e4c21fddd73199cef3e689e1c4d98494;p=samtools.git

diff --git a/bcftools/ld.c b/bcftools/ld.c
index 01145b8..0207819 100644
--- a/bcftools/ld.c
+++ b/bcftools/ld.c
@@ -70,7 +70,7 @@ double bcf_ld_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4])
 		for (i = 0; i < n_smpl; ++i) {
 			const uint8_t *pi = PL[j] + i * PL_len[j];
 			double *p = pdg[j] + i * 3;
-			p[0] = g_q2p[pi[b[j]->n_alleles]]; p[1] = g_q2p[pi[1]]; p[2] = g_q2p[pi[0]];
+			p[0] = g_q2p[pi[2]]; p[1] = g_q2p[pi[1]]; p[2] = g_q2p[pi[0]];
 		}
 	}
 	// iteration
@@ -80,7 +80,7 @@ double bcf_ld_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4])
 		memcpy(flast, f, 4 * sizeof(double));
 		freq_iter(n_smpl, pdg, f);
 		for (i = 0; i < 4; ++i) {
-			double x = fabs(f[0] - flast[0]);
+			double x = fabs(f[i] - flast[i]);
 			if (x > eps) eps = x;
 		}
 		if (eps < LD_ITER_EPS) break;
@@ -90,7 +90,7 @@ double bcf_ld_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4])
 	{ // calculate r^2
 		double p[2], q[2], D;
 		p[0] = f[0] + f[1]; q[0] = 1 - p[0];
-		p[1] = f[2] + f[3]; q[1] = 1 - p[1];
+		p[1] = f[0] + f[2]; q[1] = 1 - p[1];
 		D = f[0] * f[3] - f[1] * f[2];
 		r = sqrt(D * D / (p[0] * p[1] * q[0] * q[1]));
 		// fprintf(stderr, "R(%lf,%lf,%lf,%lf)=%lf\n", f[0], f[1], f[2], f[3], r2);
@@ -98,3 +98,46 @@ double bcf_ld_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4])
 	}
 	return r;
 }
+
+int bcf_main_pwld(int argc, char *argv[])
+{
+	bcf_t *fp;
+	bcf_hdr_t *h;
+	bcf1_t **b, *b0;
+	int i, j, m, n;
+	double f[4];
+	if (argc == 1) {
+		fprintf(stderr, "Usage: bcftools pwld <in.bcf>\n");
+		return 1;
+	}
+	fp = bcf_open(argv[1], "rb");
+	h = bcf_hdr_read(fp);
+	// read the entire BCF
+	m = n = 0; b = 0;
+	b0 = calloc(1, sizeof(bcf1_t));
+	while (bcf_read(fp, h, b0) >= 0) {
+		if (m == n) {
+			m = m? m<<1 : 16;
+			b = realloc(b, sizeof(void*) * m);
+		}
+		b[n] = calloc(1, sizeof(bcf1_t));
+		bcf_cpy(b[n++], b0);
+	}
+	bcf_destroy(b0);
+	// compute pair-wise r^2
+	printf("%d\n", n); // the number of loci
+	for (i = 0; i < n; ++i) {
+		printf("%s:%d", h->ns[b[i]->tid], b[i]->pos + 1);
+		for (j = 0; j < i; ++j) {
+			double r = bcf_ld_freq(b[i], b[j], f);
+			printf("\t%.3f", r*r);
+		}
+		printf("\t1.000\n");
+	}
+	// free
+	for (i = 0; i < n; ++i) bcf_destroy(b[i]);
+	free(b);
+	bcf_hdr_destroy(h);
+	bcf_close(fp);
+	return 0;
+}