X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fmain.c;h=07d9e93de8fa102a5453adfd5279195a813f72d2;hb=0f5ce5adb9029850cbea17c2997019459f84324c;hp=75125f872d0117d287186a72251b0e2eaeb5853f;hpb=19bf588b6b3425b1e925a2f837f0b8f351d057d7;p=samtools.git

diff --git a/bcftools/main.c b/bcftools/main.c
index 75125f8..07d9e93 100644
--- a/bcftools/main.c
+++ b/bcftools/main.c
@@ -5,7 +5,6 @@
 
 int bcfview(int argc, char *argv[]);
 int bcf_main_index(int argc, char *argv[]);
-int bcf_main_pwld(int argc, char *argv[]);
 
 #define BUF_SIZE 0x10000
 
@@ -43,6 +42,50 @@ int bcf_cat(int n, char * const *fn)
 	return 0;
 }
 
+int bcf_main_pwld(int argc, char *argv[])
+{
+	extern double bcf_pair_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4]);
+	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_pair_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;
+}
+
 int main(int argc, char *argv[])
 {
 	if (argc == 1) {