* change the order of PL/GL according to the latest VCF spec

[samtools.git] / bcftools / ld.c

diff --git a/bcftools/ld.c b/bcftools/ld.c
index 01145b84ef81538e9310a2cc61feb295dd9481d8..11474ed422aa40e5b839831a5c0b5ab80f5c73ad 100644 (file)
--- 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);