|
0
|
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 }
|
