* vcfutils.pl: fixed a typo in help message

[samtools.git] / bcftools / prob1.c

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 85c2945779f85c9cd050fd72f00cdc78f1deeacb..e3b0f7344d096f22c37048b2131219378bcb7ba5 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -106,8 +106,10 @@ int bcf_p1_read_prior(bcf_p1aux_t *ma, const char *fn)
        for (k = 0; k <= ma->M; ++k) ma->phi[k] /= sum;
        for (k = 0; k <= ma->M; ++k) fprintf(stderr, " %d:%.3lg", k, ma->phi[ma->M - k]);
        fputc('\n', stderr);
-       for (sum = 0., k = 1; k <= ma->M; ++k) sum += k * ma->phi[ma->M - k];
-       fprintf(stderr, "[heterozygosity] %lf\n", (double)sum / ma->M);
+       for (sum = 0., k = 1; k < ma->M; ++k) sum += ma->phi[ma->M - k] * (2.* k * (ma->M - k) / ma->M / (ma->M - 1));
+       fprintf(stderr, "[%s] heterozygosity=%lf, ", __func__, (double)sum);
+       for (sum = 0., k = 1; k <= ma->M; ++k) sum += k * ma->phi[ma->M - k] / ma->M;
+       fprintf(stderr, "theta=%lf\n", (double)sum);
        return 0;
 }
 
@@ -153,8 +155,6 @@ void bcf_p1_destroy(bcf_p1aux_t *ma)
        }
 }
 
-#define char2int(s) (((int)s[0])<<8|s[1])
-
 static int cal_pdg(const bcf1_t *b, bcf_p1aux_t *ma)
 {
        int i, j, k;
@@ -281,7 +281,7 @@ static void contrast(bcf_p1aux_t *ma, double pc[4]) // mc_cal_y() must be called
 {
        int k, n1 = ma->n1, n2 = ma->n - ma->n1;
        long double sum1, sum2;
-       pc[0] = pc[1] = pc[2] = pc[3] = 0.;
+       pc[0] = pc[1] = pc[2] = pc[3] = -1.;
        if (n1 <= 0 || n2 <= 0) return;
        for (k = 0, sum1 = 0.; k <= 2*n1; ++k) sum1 += ma->phi1[k] * ma->z1[k];
        for (k = 0, sum2 = 0.; k <= 2*n2; ++k) sum2 += ma->phi2[k] * ma->z2[k];
@@ -350,7 +350,7 @@ int bcf_p1_cal(bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst)
        long double sum = 0.;
        // set PL and PL_len
        for (i = 0; i < b->n_gi; ++i) {
-               if (b->gi[i].fmt == char2int("PL")) {
+               if (b->gi[i].fmt == bcf_str2int("PL", 2)) {
                        ma->PL = (uint8_t*)b->gi[i].data;
                        ma->PL_len = b->gi[i].len;
                        break;