Mercurial > repos > siyuan > prada
comparison pyPRADA_1.2/tools/samtools-0.1.16/errmod.c @ 0:acc2ca1a3ba4
Uploaded
author | siyuan |
---|---|
date | Thu, 20 Feb 2014 00:44:58 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:acc2ca1a3ba4 |
---|---|
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 } |