3 #include "bam_maqcns.h"
6 KSORT_INIT_GENERIC(uint32_t)
8 #define INDEL_WINDOW_SIZE 50
10 typedef struct __bmc_aux_t {
16 float esum[4], fsum[4];
21 char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
24 P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
25 P(D|<b1,b2>) = \sum_{k=1}^{N-1} p_k 1/2 [(k/N)^n_2(1-k/N)^n_1 + (k/N)^n1(1-k/N)^n_2]
26 p_k = i/k / \sum_{i=1}^{N-1} 1/i
28 static void cal_het(bam_maqcns_t *aa)
31 double sum_harmo; // harmonic sum
33 double p1 = 0.0, p3 = 0.0; // just for testing
36 aa->lhet = (double*)calloc(256 * 256, sizeof(double));
38 for (k = 1; k <= aa->n_hap - 1; ++k)
40 for (n1 = 0; n1 < 256; ++n1) {
41 for (n2 = 0; n2 < 256; ++n2) {
42 long double sum = 0.0;
43 double lC = lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
44 for (k = 1; k <= aa->n_hap - 1; ++k) {
45 double pk = 1.0 / k / sum_harmo;
46 double log1 = log((double)k/aa->n_hap);
47 double log2 = log(1.0 - (double)k/aa->n_hap);
48 sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
50 aa->lhet[n1<<8|n2] = lC + logl(sum);
51 if (n1 == 17 && n2 == 3) p3 = lC + logl(expl(logl(0.5) * 20));
52 if (n1 == 19 && n2 == 1) p1 = lC + logl(expl(logl(0.5) * 20));
55 poly_rate = aa->het_rate * sum_harmo;
56 aa->q_r = -4.343 * log(2.0 * poly_rate / (1.0 - poly_rate));
59 /** initialize the helper structure */
60 static void cal_coef(bam_maqcns_t *aa)
63 long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
66 lC = (double*)calloc(256 * 256, sizeof(double));
67 // aa->lhet will be allocated and initialized
68 free(aa->fk); free(aa->coef);
69 aa->fk = (double*)calloc(256, sizeof(double));
70 aa->coef = (double*)calloc(256*256*64, sizeof(double));
71 aa->fk[0] = fk2[0] = 1.0;
72 for (n = 1; n != 256; ++n) {
73 aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
74 fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
76 for (n = 1; n != 256; ++n)
77 for (k = 1; k <= n; ++k)
78 lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
79 for (q = 1; q != 64; ++q) {
80 double e = pow(10.0, -q/10.0);
82 double le1 = log(1.0-e);
83 for (n = 1; n != 256; ++n) {
84 double *coef = aa->coef + (q<<16|n<<8);
86 for (k = n; k >= 0; --k) { // a_k = \sum_{i=k}^n C^n_k \epsilon^k (1-\epsilon)^{n-k}
87 sum_a[k] = sum_a[k+1] + expl(lC[n<<8|k] + k*le + (n-k)*le1);
88 b[k] = sum_a[k+1] / sum_a[k];
89 if (b[k] > 0.99) b[k] = 0.99;
91 for (k = 0; k != n; ++k) // log(\bar\beta_{nk}(\bar\epsilon)^{f_k})
92 q_c[k] = -4.343 * fk2[k] * logl(b[k] / e);
93 for (k = 1; k != n; ++k) q_c[k] += q_c[k-1]; // \prod_{i=0}^k c_i
94 for (k = 0; k <= n; ++k) { // powl() in 64-bit mode seems broken on my Mac OS X 10.4.9
95 tmp[k] = -4.343 * logl(1.0 - expl(fk2[k] * logl(b[k])));
96 coef[k] = (k? q_c[k-1] : 0) + tmp[k]; // this is the final c_{nk}
103 bam_maqcns_t *bam_maqcns_init()
106 bm = (bam_maqcns_t*)calloc(1, sizeof(bam_maqcns_t));
107 bm->aux = (bmc_aux_t*)calloc(1, sizeof(bmc_aux_t));
108 bm->het_rate = 0.001;
116 void bam_maqcns_prepare(bam_maqcns_t *bm)
118 cal_coef(bm); cal_het(bm);
121 void bam_maqcns_destroy(bam_maqcns_t *bm)
124 free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info);
125 free(bm->aux); free(bm);
128 glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
131 int i, j, k, w[8], c, n;
132 glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
133 float p[16], min_p = 1e30;
136 g->ref_base = ref_base;
137 if (_n == 0) return g;
139 // construct aux array
140 if (bm->aux->max < _n) {
142 kroundup32(bm->aux->max);
143 bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
145 for (i = n = 0; i < _n; ++i) {
146 const bam_pileup1_t *p = pl + i;
147 uint32_t q, x = 0, qq;
148 if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue;
149 q = (uint32_t)bam1_qual(p->b)[p->qpos];
150 x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
151 if (p->b->core.qual < q) q = p->b->core.qual;
153 qq = bam1_seqi(bam1_seq(p->b), p->qpos);
154 q = bam_nt16_nt4_table[qq? qq : ref_base];
155 if (!p->is_del && q < 4) x |= 1 << 21 | q << 16;
156 bm->aux->info[n++] = x;
158 ks_introsort(uint32_t, n, bm->aux->info);
159 // generate esum and fsum
160 b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
161 for (k = 0; k != 8; ++k) w[k] = 0;
163 for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
164 uint32_t info = bm->aux->info[j];
166 if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
169 b->esum[k&3] += bm->fk[w[k]] * (info>>24);
170 b->fsum[k&3] += bm->fk[w[k]];
171 if (w[k] < 0xff) ++w[k];
174 tmp = (int)(info&0xff) < bm->cap_mapQ? (int)(info&0xff) : bm->cap_mapQ;
177 b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
179 for (j = c = 0; j != 4; ++j) c += b->c[j];
181 for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
182 for (j = c = 0; j != 4; ++j) c += b->c[j];
184 // generate likelihood
185 for (j = 0; j != 4; ++j) {
189 for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
190 if (j == k) continue;
191 tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
194 bar_e = (int)(tmp1 / tmp3 + 0.5);
195 if (bar_e < 4) bar_e = 4; // should not happen
196 if (bar_e > 63) bar_e = 63;
197 p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
198 } else p[j<<2|j] = 0.0; // all the bases are j
200 for (k = j + 1; k < 4; ++k) {
201 for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
202 if (i == j || i == k) continue;
203 tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
206 bar_e = (int)(tmp1 / tmp3 + 0.5);
207 if (bar_e < 4) bar_e = 4;
208 if (bar_e > 63) bar_e = 63;
209 p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
210 } else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k
213 for (k = 0; k != 4; ++k)
214 if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
218 float max1, max2, min1, min2;
221 max1 = max2 = -1.0; min1 = min2 = 1e30;
222 for (k = 0; k < 4; ++k) {
223 if (b->esum[k] > max1) {
224 max2 = max1; max1 = b->esum[k]; max_k = k;
225 } else if (b->esum[k] > max2) max2 = b->esum[k];
227 for (k = 0; k < 4; ++k) {
228 if (p[k<<2|k] < min1) {
229 min2 = min1; min1 = p[k<<2|k]; min_k = k;
230 } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
232 if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
233 p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
236 // convert necessary information to glf1_t
237 g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ;
238 g->depth = n > 16777215? 16777215 : n;
239 for (j = 0; j != 4; ++j)
240 for (k = j; k < 4; ++k)
241 if (p[j<<2|k] < min_p) min_p = p[j<<2|k];
242 g->min_lk = min_p > 255.0? 255 : (int)(min_p + 0.5);
243 for (j = c = 0; j != 4; ++j)
244 for (k = j; k < 4; ++k)
245 g->lk[c++] = p[j<<2|k]-min_p > 255.0? 255 : (int)(p[j<<2|k]-min_p + 0.5);
251 uint32_t glf2cns(const glf1_t *g, int q_r)
253 int i, j, k, tmp[16], min = 10000, min2 = 10000, min3 = 10000, min_g = -1, min_g2 = -1;
255 for (i = k = 0; i < 4; ++i)
256 for (j = i; j < 4; ++j) {
258 tmp[i<<2|j] = g->lk[k++] + (i == j? 0 : q_r);
260 for (i = 0; i < 16; ++i) {
261 if (tmp[i] < 0) continue;
263 min3 = min2; min2 = min; min = tmp[i]; min_g2 = min_g; min_g = i;
264 } else if (tmp[i] < min2) {
265 min3 = min2; min2 = tmp[i]; min_g2 = i;
266 } else if (tmp[i] < min3) min3 = tmp[i];
268 x = min_g >= 0? (1U<<(min_g>>2&3) | 1U<<(min_g&3)) << 28 : 0xf << 28;
269 x |= min_g2 >= 0? (1U<<(min_g2>>2&3) | 1U<<(min_g2&3)) << 24 : 0xf << 24;
270 x |= (uint32_t)g->max_mapQ << 16;
271 x |= min2 < 10000? (min2 - min < 256? min2 - min : 255) << 8 : 0xff << 8;
272 x |= min2 < 10000 && min3 < 10000? (min3 - min2 < 256? min3 - min2 : 255) : 0xff;
276 uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm)
281 g = bam_maqcns_glfgen(n, pl, 0xf, bm);
282 x = glf2cns(g, (int)(bm->q_r + 0.5));
284 } else x = 0xfU<<28 | 0xfU<<24;
288 /************** *****************/
290 bam_maqindel_opt_t *bam_maqindel_opt_init()
292 bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
294 mi->r_indel = 0.00015;
302 void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
304 if (mir == 0) return;
305 free(mir->s[0]); free(mir->s[1]); free(mir);
308 #define MINUS_CONST 0x10000000
310 bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
311 int _n_types, int *_types)
313 int i, j, n_types, *types, left, right, max_rd_len = 0;
314 bam_maqindel_ret_t *ret = 0;
315 // if there is no proposed indel, check if there is an indel from the alignment
317 for (i = 0; i < n; ++i) {
318 const bam_pileup1_t *p = pl + i;
319 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
321 if (i == n) return 0; // no indel
323 { // calculate how many types of indels are available (set n_types and types)
326 aux = (uint32_t*)calloc(n + _n_types + 1, 4);
328 aux[m++] = MINUS_CONST; // zero indel is always a type
329 for (i = 0; i < n; ++i) {
330 const bam_pileup1_t *p = pl + i;
331 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
332 aux[m++] = MINUS_CONST + p->indel;
333 j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
334 if (j > max_rd_len) max_rd_len = j;
336 if (_n_types) // then also add this to aux[]
337 for (i = 0; i < _n_types; ++i)
338 if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
339 ks_introsort(uint32_t, m, aux);
340 // squeeze out identical types
341 for (i = 1, n_types = 1; i < m; ++i)
342 if (aux[i] != aux[i-1]) ++n_types;
343 types = (int*)calloc(n_types, sizeof(int));
345 types[j++] = aux[0] - MINUS_CONST;
346 for (i = 1; i < m; ++i) {
347 if (aux[i] != aux[i-1])
348 types[j++] = aux[i] - MINUS_CONST;
352 { // calculate left and right boundary
353 left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
354 right = pos + INDEL_WINDOW_SIZE;
357 char *ref2, *rs, *inscns = 0;
358 int k, l, *score, *pscore, max_ins = types[n_types-1];
359 ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
360 rs = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1);
361 if (max_ins > 0) { // get the consensus of inserted sequences
362 int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
364 for (i = 0; i < n_types; ++i) {
365 if (types[i] <= 0) continue; // not insertion
366 for (j = 0; j < n; ++j) {
367 const bam_pileup1_t *p = pl + j;
368 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel == types[i]) {
369 for (k = 1; k <= p->indel; ++k) {
370 int c = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos + k)];
371 if (c < 4) ++inscns_aux[i*max_ins*4 + (k-1)*4 + c];
376 // construct the consensus of inserted sequence
377 inscns = (char*)calloc(n_types * max_ins, sizeof(char));
378 for (i = 0; i < n_types; ++i) {
379 for (j = 0; j < types[i]; ++j) {
380 int max = 0, max_k = -1, *ia = inscns_aux + i*max_ins*4 + j*4;
381 for (k = 0; k < 4; ++k) {
387 inscns[i*max_ins + j] = max? 1<<max_k : 15;
393 score = (int*)calloc(n_types * n, sizeof(int));
394 pscore = (int*)calloc(n_types * n, sizeof(int));
395 for (i = 0; i < n_types; ++i) {
396 ka_param_t ap = ka_param_blast;
397 ap.band_width = 2 * types[n_types - 1] + 2;
399 for (k = 0, j = left; j <= pos; ++j)
400 ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
401 if (types[i] <= 0) j += -types[i];
402 else for (l = 0; l < types[i]; ++l)
403 ref2[k++] = inscns[i*max_ins + l];
404 for (; j < right && ref[j]; ++j)
405 ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
406 if (j < right) right = j;
407 // calculate score for each read
408 for (j = 0; j < n; ++j) {
409 const bam_pileup1_t *p = pl + j;
410 int qbeg, qend, tbeg, tend;
411 if (p->b->core.flag & BAM_FUNMAP) continue;
412 qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left, &tbeg);
413 qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, &tend);
414 for (l = qbeg; l < qend; ++l)
415 rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
417 int x, y, n_acigar, ps;
420 acigar = ka_global_core((uint8_t*)ref2 + tbeg - left, tend - tbeg + types[i], (uint8_t*)rs, qend - qbeg, &ap, &score[i*n+j], &n_acigar);
421 x = tbeg - left; y = 0;
422 for (l = 0; l < n_acigar; ++l) {
423 int op = acigar[l]&0xf;
424 int len = acigar[l]>>4;
425 if (op == BAM_CMATCH) {
427 for (k = 0; k < len; ++k)
428 if (ref2[x+k] != rs[y+k]) ps += bam1_qual(p->b)[y+k];
430 } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
431 if (op == BAM_CINS) ps += mi->q_indel * len;
433 } else if (op == BAM_CDEL) {
434 ps += mi->q_indel * len;
439 //printf("pos=%d, type=%d, j=%d, score=%d, psore=%d, %d, %d, %d, %d, ", pos+1, types[i], j, score[i*n+j], pscore[i*n+j], tbeg, tend, qbeg, qend);
440 //for (l = 0; l < n_acigar; ++l) printf("%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]); printf("\n");
445 { // get final result
446 int *sum, max1, max2, max1_i, max2_i;
447 // pick up the best two score
448 sum = (int*)calloc(n_types, sizeof(int));
449 for (i = 0; i < n_types; ++i)
450 for (j = 0; j < n; ++j)
451 sum[i] += -pscore[i*n+j];
452 max1 = max2 = -0x7fffffff; max1_i = max2_i = -1;
453 for (i = 0; i < n_types; ++i) {
455 max2 = max1; max2_i = max1_i; max1 = sum[i]; max1_i = i;
456 } else if (sum[i] > max2) {
457 max2 = sum[i]; max2_i = i;
462 ret = (bam_maqindel_ret_t*)calloc(1, sizeof(bam_maqindel_ret_t));
463 ret->indel1 = types[max1_i]; ret->indel2 = types[max2_i];
464 ret->s[0] = (char*)calloc(abs(ret->indel1) + 2, 1);
465 ret->s[1] = (char*)calloc(abs(ret->indel2) + 2, 1);
466 // write indel sequence
467 if (ret->indel1 > 0) {
469 for (k = 0; k < ret->indel1; ++k)
470 ret->s[0][k+1] = bam_nt16_rev_table[(int)inscns[max1_i*max_ins + k]];
471 } else if (ret->indel1 < 0) {
473 for (k = 0; k < -ret->indel1 && ref[pos + k + 1]; ++k)
474 ret->s[0][k+1] = ref[pos + k + 1];
475 } else ret->s[0][0] = '*';
476 if (ret->indel2 > 0) {
478 for (k = 0; k < ret->indel2; ++k)
479 ret->s[1][k+1] = bam_nt16_rev_table[(int)inscns[max2_i*max_ins + k]];
480 } else if (ret->indel2 < 0) {
482 for (k = 0; k < -ret->indel2 && ref[pos + k + 1]; ++k)
483 ret->s[1][k+1] = ref[pos + k + 1];
484 } else ret->s[1][0] = '*';
486 for (i = 0; i < n; ++i) {
487 const bam_pileup1_t *p = pl + i;
488 if (p->indel == ret->indel1) ++ret->cnt1;
489 else if (p->indel == ret->indel2) ++ret->cnt2;
490 else ++ret->cnt_anti;
493 ret->gl[0] = ret->gl[1] = 0;
494 for (j = 0; j < n; ++j) {
495 int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
496 //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2);
497 if (s1 > s2) ret->gl[0] += s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel;
498 else ret->gl[1] += s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel;
500 // write cnt_ref and cnt_ambi
501 if (max1_i != 0 && max2_i != 0) {
502 for (j = 0; j < n; ++j) {
503 int diff1 = score[j] - score[max1_i * n + j];
504 int diff2 = score[j] - score[max2_i * n + j];
505 if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
506 else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
510 free(score); free(pscore); free(ref2); free(rs); free(inscns);
513 int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5);
514 int min1, min2, min1_i;
515 q[0] = ret->gl[0] + (ret->s[0][0] != '*'? 0 : 0) * qr_indel;
516 q[1] = ret->gl[1] + (ret->s[1][0] != '*'? 0 : 0) * qr_indel;
517 q[2] = n * 3 + (ret->s[0][0] == '*' || ret->s[1][0] == '*'? 1 : 1) * qr_indel;
518 min1 = min2 = 0x7fffffff; min1_i = -1;
519 for (i = 0; i < 3; ++i) {
521 min2 = min1; min1 = q[i]; min1_i = i;
522 } else if (q[i] < min2) min2 = q[i];
525 ret->q_cns = min2 - min1;
527 if (ret->gt < 2) ret->q_ref = (ret->s[ret->gt][0] == '*')? 0 : q[1-ret->gt] - q[ret->gt] - qr_indel - 3;
528 else ret->q_ref = (ret->s[0][0] == '*')? q[0] - q[2] : q[1] - q[2];
529 if (ret->q_ref < 0) ret->q_ref = 0;