← Back to branch summary

~ubuntu-branches/ubuntu/quantal/genometools/quantal-backports

« back to all changes in this revision

Viewing changes to src/external/samtools-0.1.18/errmod.c

  • Committer: Package Import Robot
  • Author(s): Sascha Steinbiss
  • Date: 2012-07-09 14:10:23 UTC
  • Revision ID: package-import@ubuntu.com-20120709141023-juuu4spm6chqsf9o
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/external/samtools-0.1.18/errmod.c
 
1
#include <math.h>
 
2
#include "errmod.h"
 
3
#include "ksort.h"
 
4
KSORT_INIT_GENERIC(uint16_t)
 
5
 
 
6
typedef struct __errmod_coef_t {
 
7
        double *fk, *beta, *lhet;
 
8
} errmod_coef_t;
 
9
 
 
10
typedef struct {
 
11
        double fsum[16], bsum[16];
 
12
        uint32_t c[16];
 
13
} call_aux_t;
 
14
 
 
15
static errmod_coef_t *cal_coef(double depcorr, double eta)
 
16
{
 
17
        int k, n, q;
 
18
        long double sum, sum1;
 
19
        double *lC;
 
20
        errmod_coef_t *ec;
 
21
 
 
22
        ec = calloc(1, sizeof(errmod_coef_t));
 
23
        // initialize ->fk
 
24
        ec->fk = (double*)calloc(256, sizeof(double));
 
25
        ec->fk[0] = 1.0;
 
26
        for (n = 1; n != 256; ++n)
 
27
                ec->fk[n] = pow(1. - depcorr, n) * (1.0 - eta) + eta;
 
28
        // initialize ->coef
 
29
        ec->beta = (double*)calloc(256 * 256 * 64, sizeof(double));
 
30
        lC = (double*)calloc(256 * 256, sizeof(double));
 
31
        for (n = 1; n != 256; ++n) {
 
32
                double lgn = lgamma(n+1);
 
33
                for (k = 1; k <= n; ++k)
 
34
                        lC[n<<8|k] = lgn - lgamma(k+1) - lgamma(n-k+1);
 
35
        }
 
36
        for (q = 1; q != 64; ++q) {
 
37
                double e = pow(10.0, -q/10.0);
 
38
                double le = log(e);
 
39
                double le1 = log(1.0 - e);
 
40
                for (n = 1; n <= 255; ++n) {
 
41
                        double *beta = ec->beta + (q<<16|n<<8);
 
42
                        sum1 = sum = 0.0;
 
43
                        for (k = n; k >= 0; --k, sum1 = sum) {
 
44
                                sum = sum1 + expl(lC[n<<8|k] + k*le + (n-k)*le1);
 
45
                                beta[k] = -10. / M_LN10 * logl(sum1 / sum);
 
46
                        }
 
47
                }
 
48
        }
 
49
        // initialize ->lhet
 
50
        ec->lhet = (double*)calloc(256 * 256, sizeof(double));
 
51
        for (n = 0; n < 256; ++n)
 
52
                for (k = 0; k < 256; ++k)
 
53
                        ec->lhet[n<<8|k] = lC[n<<8|k] - M_LN2 * n;
 
54
        free(lC);
 
55
        return ec;
 
56
}
 
57
 
 
58
errmod_t *errmod_init(float depcorr)
 
59
{
 
60
        errmod_t *em;
 
61
        em = (errmod_t*)calloc(1, sizeof(errmod_t));
 
62
        em->depcorr = depcorr;
 
63
        em->coef = cal_coef(depcorr, 0.03);
 
64
        return em;
 
65
}
 
66
 
 
67
void errmod_destroy(errmod_t *em)
 
68
{
 
69
        if (em == 0) return;
 
70
        free(em->coef->lhet); free(em->coef->fk); free(em->coef->beta);
 
71
        free(em->coef); free(em);
 
72
}
 
73
// qual:6, strand:1, base:4
 
74
int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q)
 
75
{
 
76
        call_aux_t aux;
 
77
        int i, j, k, w[32];
 
78
 
 
79
        if (m > m) return -1;
 
80
        memset(q, 0, m * m * sizeof(float));
 
81
        if (n == 0) return 0;
 
82
        // calculate aux.esum and aux.fsum
 
83
        if (n > 255) { // then sample 255 bases
 
84
                ks_shuffle(uint16_t, n, bases);
 
85
                n = 255;
 
86
        }
 
87
        ks_introsort(uint16_t, n, bases);
 
88
        memset(w, 0, 32 * sizeof(int));
 
89
        memset(&aux, 0, sizeof(call_aux_t));
 
90
        for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
 
91
                uint16_t b = bases[j];
 
92
                int q = b>>5 < 4? 4 : b>>5;
 
93
                if (q > 63) q = 63;
 
94
                k = b&0x1f;
 
95
                aux.fsum[k&0xf] += em->coef->fk[w[k]];
 
96
                aux.bsum[k&0xf] += em->coef->fk[w[k]] * em->coef->beta[q<<16|n<<8|aux.c[k&0xf]];
 
97
                ++aux.c[k&0xf];
 
98
                ++w[k];
 
99
        }
 
100
        // generate likelihood
 
101
        for (j = 0; j != m; ++j) {
 
102
                float tmp1, tmp3;
 
103
                int tmp2, bar_e;
 
104
                // homozygous
 
105
                for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != m; ++k) {
 
106
                        if (k == j) continue;
 
107
                        tmp1 += aux.bsum[k]; tmp2 += aux.c[k]; tmp3 += aux.fsum[k];
 
108
                }
 
109
                if (tmp2) {
 
110
                        bar_e = (int)(tmp1 / tmp3 + 0.499);
 
111
                        if (bar_e > 63) bar_e = 63;
 
112
                        q[j*m+j] = tmp1;
 
113
                }
 
114
                // heterozygous
 
115
                for (k = j + 1; k < m; ++k) {
 
116
                        int cjk = aux.c[j] + aux.c[k];
 
117
                        for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i < m; ++i) {
 
118
                                if (i == j || i == k) continue;
 
119
                                tmp1 += aux.bsum[i]; tmp2 += aux.c[i]; tmp3 += aux.fsum[i];
 
120
                        }
 
121
                        if (tmp2) {
 
122
                                bar_e = (int)(tmp1 / tmp3 + 0.499);
 
123
                                if (bar_e > 63) bar_e = 63;
 
124
                                q[j*m+k] = q[k*m+j] = -4.343 * em->coef->lhet[cjk<<8|aux.c[k]] + tmp1;
 
125
                        } else q[j*m+k] = q[k*m+j] = -4.343 * em->coef->lhet[cjk<<8|aux.c[k]]; // all the bases are either j or k
 
126
                }
 
127
                for (k = 0; k != m; ++k) if (q[j*m+k] < 0.0) q[j*m+k] = 0.0;
 
128
        }
 
129
        return 0;
 
130
}