← Back to branch summary

~ubuntu-branches/ubuntu/raring/samtools/raring

« back to all changes in this revision

Viewing changes to bam_maqcns.c

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-11-17 21:38:24 UTC
  • Revision ID: james.westby@ubuntu.com-20091117213824-dfouynpy3r7ismpj
Tags: 0.1.7a~dfsg-1
* New upstream release: new script sam2vcf.pl, and many other changes.
* Package converted to the format ‘3.0 (quilt)’ (debian/source/format).
* Wrote a manual page for razip (debian/razip.1).
* Better clean the example directory to make the source package
  buildable twice in a row (debian/rules).

Show diffs side-by-side

added added

removed removed

Lines of Context:
bam_maqcns.c
1
1
#include <math.h>
 
2
#include <assert.h>
2
3
#include "bam.h"
3
4
#include "bam_maqcns.h"
4
5
#include "ksort.h"
 
6
#include "kaln.h"
5
7
KSORT_INIT_GENERIC(uint32_t)
6
8
 
7
 
#define MAX_WINDOW 33
 
9
#define INDEL_WINDOW_SIZE 50
 
10
#define INDEL_EXT_DEP 0.9
8
11
 
9
12
typedef struct __bmc_aux_t {
10
13
        int max;
22
25
/*
23
26
  P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
24
27
  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]
25
 
  p_k = i/k / \sum_{i=1}^{N-1} 1/i
 
28
  p_k = 1/k / \sum_{i=1}^{N-1} 1/i
26
29
 */
27
30
static void cal_het(bam_maqcns_t *aa)
28
31
{
29
32
        int k, n1, n2;
30
33
        double sum_harmo; // harmonic sum
31
34
        double poly_rate;
32
 
        double p1 = 0.0, p3 = 0.0; // just for testing
33
35
 
34
36
        free(aa->lhet);
35
37
        aa->lhet = (double*)calloc(256 * 256, sizeof(double));
39
41
        for (n1 = 0; n1 < 256; ++n1) {
40
42
                for (n2 = 0; n2 < 256; ++n2) {
41
43
                        long double sum = 0.0;
42
 
                        double lC = lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
 
44
                        double lC = aa->is_soap? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
43
45
                        for (k = 1; k <= aa->n_hap - 1; ++k) {
44
46
                                double pk = 1.0 / k / sum_harmo;
45
47
                                double log1 = log((double)k/aa->n_hap);
47
49
                                sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
48
50
                        }
49
51
                        aa->lhet[n1<<8|n2] = lC + logl(sum);
50
 
                        if (n1 == 17 && n2 == 3) p3 = lC + logl(expl(logl(0.5) * 20));
51
 
                        if (n1 == 19 && n2 == 1) p1 = lC + logl(expl(logl(0.5) * 20));
52
52
                }
53
53
        }
54
54
        poly_rate = aa->het_rate * sum_harmo;
62
62
        long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
63
63
        double *lC;
64
64
 
65
 
        lC = (double*)calloc(256 * 256, sizeof(double));
66
65
        // aa->lhet will be allocated and initialized 
67
66
        free(aa->fk); free(aa->coef);
 
67
        aa->coef = 0;
68
68
        aa->fk = (double*)calloc(256, sizeof(double));
69
 
        aa->coef = (double*)calloc(256*256*64, sizeof(double));
70
69
        aa->fk[0] = fk2[0] = 1.0;
71
70
        for (n = 1; n != 256; ++n) {
72
71
                aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
73
72
                fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
74
73
        }
 
74
        if (aa->is_soap) return;
 
75
        aa->coef = (double*)calloc(256*256*64, sizeof(double));
 
76
        lC = (double*)calloc(256 * 256, sizeof(double));
75
77
        for (n = 1; n != 256; ++n)
76
78
                for (k = 1; k <= n; ++k)
77
79
                        lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
170
172
                        if (w[k] < 0xff) ++w[k];
171
173
                        ++b->c[k&3];
172
174
                }
173
 
                tmp = (int)(info&0x7f) < bm->cap_mapQ? (int)(info&0x7f) : bm->cap_mapQ;
 
175
                tmp = (int)(info&0xff) < bm->cap_mapQ? (int)(info&0xff) : bm->cap_mapQ;
174
176
                rms += tmp * tmp;
175
177
        }
176
178
        b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
180
182
                for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
181
183
                for (j = c = 0; j != 4; ++j) c += b->c[j];
182
184
        }
183
 
        // generate likelihood
184
 
        for (j = 0; j != 4; ++j) {
185
 
                // homozygous
186
 
                float tmp1, tmp3;
187
 
                int tmp2, bar_e;
188
 
                for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
189
 
                        if (j == k) continue;
190
 
                        tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
191
 
                }
192
 
                if (tmp2) {
193
 
                        bar_e = (int)(tmp1 / tmp3 + 0.5);
194
 
                        if (bar_e < 4) bar_e = 4; // should not happen
195
 
                        if (bar_e > 63) bar_e = 63;
196
 
                        p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
197
 
                } else p[j<<2|j] = 0.0; // all the bases are j
198
 
                // heterozygous
199
 
                for (k = j + 1; k < 4; ++k) {
200
 
                        for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
201
 
                                if (i == j || i == k) continue;
202
 
                                tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
 
185
        if (!bm->is_soap) {
 
186
                // generate likelihood
 
187
                for (j = 0; j != 4; ++j) {
 
188
                        // homozygous
 
189
                        float tmp1, tmp3;
 
190
                        int tmp2, bar_e;
 
191
                        for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
 
192
                                if (j == k) continue;
 
193
                                tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
203
194
                        }
204
195
                        if (tmp2) {
205
196
                                bar_e = (int)(tmp1 / tmp3 + 0.5);
206
 
                                if (bar_e < 4) bar_e = 4;
 
197
                                if (bar_e < 4) bar_e = 4; // should not happen
207
198
                                if (bar_e > 63) bar_e = 63;
208
 
                                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];
209
 
                        } 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
 
199
                                p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
 
200
                        } else p[j<<2|j] = 0.0; // all the bases are j
 
201
                        // heterozygous
 
202
                        for (k = j + 1; k < 4; ++k) {
 
203
                                for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
 
204
                                        if (i == j || i == k) continue;
 
205
                                        tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
 
206
                                }
 
207
                                if (tmp2) {
 
208
                                        bar_e = (int)(tmp1 / tmp3 + 0.5);
 
209
                                        if (bar_e < 4) bar_e = 4;
 
210
                                        if (bar_e > 63) bar_e = 63;
 
211
                                        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];
 
212
                                } 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
                        }
 
214
                        //
 
215
                        for (k = 0; k != 4; ++k)
 
216
                                if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
210
217
                }
211
 
                //
212
 
                for (k = 0; k != 4; ++k)
213
 
                        if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
214
 
        }
215
218
 
216
 
        { // fix p[k<<2|k]
217
 
                float max1, max2, min1, min2;
218
 
                int max_k, min_k;
219
 
                max_k = min_k = -1;
220
 
                max1 = max2 = -1.0; min1 = min2 = 1e30;
221
 
                for (k = 0; k < 4; ++k) {
222
 
                        if (b->esum[k] > max1) {
223
 
                                max2 = max1; max1 = b->esum[k]; max_k = k;
224
 
                        } else if (b->esum[k] > max2) max2 = b->esum[k];
225
 
                }
226
 
                for (k = 0; k < 4; ++k) {
227
 
                        if (p[k<<2|k] < min1) {
228
 
                                min2 = min1; min1 = p[k<<2|k]; min_k = k;
229
 
                        } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
230
 
                }
231
 
                if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
232
 
                        p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
 
219
                { // fix p[k<<2|k]
 
220
                        float max1, max2, min1, min2;
 
221
                        int max_k, min_k;
 
222
                        max_k = min_k = -1;
 
223
                        max1 = max2 = -1.0; min1 = min2 = 1e30;
 
224
                        for (k = 0; k < 4; ++k) {
 
225
                                if (b->esum[k] > max1) {
 
226
                                        max2 = max1; max1 = b->esum[k]; max_k = k;
 
227
                                } else if (b->esum[k] > max2) max2 = b->esum[k];
 
228
                        }
 
229
                        for (k = 0; k < 4; ++k) {
 
230
                                if (p[k<<2|k] < min1) {
 
231
                                        min2 = min1; min1 = p[k<<2|k]; min_k = k;
 
232
                                } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
 
233
                        }
 
234
                        if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
 
235
                                p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
 
236
                }
 
237
        } else { // apply the SOAP model
 
238
                // generate likelihood
 
239
                for (j = 0; j != 4; ++j) {
 
240
                        float tmp;
 
241
                        // homozygous
 
242
                        for (k = 0, tmp = 0.0; k != 4; ++k)
 
243
                                if (j != k) tmp += b->esum[k];
 
244
                        p[j<<2|j] = tmp;
 
245
                        // heterozygous
 
246
                        for (k = j + 1; k < 4; ++k) {
 
247
                                for (i = 0, tmp = 0.0; i != 4; ++i)
 
248
                                        if (i != j && i != k) tmp += b->esum[i];
 
249
                                p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp;
 
250
                        }
 
251
                }
233
252
        }
234
253
 
235
254
        // convert necessary information to glf1_t
304
323
        free(mir->s[0]); free(mir->s[1]); free(mir);
305
324
}
306
325
 
 
326
int bam_tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
 
327
{
 
328
        int k, x = c->pos, y = 0, last_y = 0;
 
329
        *_tpos = c->pos;
 
330
        for (k = 0; k < c->n_cigar; ++k) {
 
331
                int op = cigar[k] & BAM_CIGAR_MASK;
 
332
                int l = cigar[k] >> BAM_CIGAR_SHIFT;
 
333
                if (op == BAM_CMATCH) {
 
334
                        if (c->pos > tpos) return y;
 
335
                        if (x + l > tpos) {
 
336
                                *_tpos = tpos;
 
337
                                return y + (tpos - x);
 
338
                        }
 
339
                        x += l; y += l;
 
340
                        last_y = y;
 
341
                } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
 
342
                else if (op == BAM_CDEL || op == BAM_CREF_SKIP) {
 
343
                        if (x + l > tpos) {
 
344
                                *_tpos = is_left? x : x + l;
 
345
                                return y;
 
346
                        }
 
347
                        x += l;
 
348
                }
 
349
        }
 
350
        *_tpos = x;
 
351
        return last_y;
 
352
}
 
353
 
307
354
#define MINUS_CONST 0x10000000
308
355
 
309
356
bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
310
357
                                                                 int _n_types, int *_types)
311
358
{
312
 
        int i, j, n_types, *types, left, right;
 
359
        int i, j, n_types, *types, left, right, max_rd_len = 0;
313
360
        bam_maqindel_ret_t *ret = 0;
314
361
        // if there is no proposed indel, check if there is an indel from the alignment
315
362
        if (_n_types == 0) {
329
376
                        const bam_pileup1_t *p = pl + i;
330
377
                        if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
331
378
                                aux[m++] = MINUS_CONST + p->indel;
 
379
                        j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
 
380
                        if (j > max_rd_len) max_rd_len = j;
332
381
                }
333
382
                if (_n_types) // then also add this to aux[]
334
383
                        for (i = 0; i < _n_types; ++i)
347
396
                free(aux);
348
397
        }
349
398
        { // calculate left and right boundary
350
 
                bam_segreg_t seg;
351
 
                left = 0x7fffffff; right = 0;
352
 
                for (i = 0; i < n; ++i) {
353
 
                        const bam_pileup1_t *p = pl + i;
354
 
                        if (!(p->b->core.flag&BAM_FUNMAP)) {
355
 
                                bam_segreg(pos, &p->b->core, bam1_cigar(p->b), &seg);
356
 
                                if (seg.tbeg < left) left = seg.tbeg;
357
 
                                if (seg.tend > right) right = seg.tend;
358
 
                        }
359
 
                }
360
 
                if (pos - left > MAX_WINDOW) left = pos - MAX_WINDOW;
361
 
                if (right - pos> MAX_WINDOW) right = pos + MAX_WINDOW;
 
399
                left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
 
400
                right = pos + INDEL_WINDOW_SIZE;
 
401
                if (types[0] < 0) right -= types[0];
 
402
                // in case the alignments stand out the reference
 
403
                for (i = pos; i < right; ++i)
 
404
                        if (ref[i] == 0) break;
 
405
                right = i;
362
406
        }
363
407
        { // the core part
364
 
                char *ref2, *inscns = 0;
 
408
                char *ref2, *rs, *inscns = 0;
365
409
                int k, l, *score, *pscore, max_ins = types[n_types-1];
366
 
                ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
367
410
                if (max_ins > 0) { // get the consensus of inserted sequences
368
411
                        int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
369
412
                        // count occurrences
396
439
                        free(inscns_aux);
397
440
                }
398
441
                // calculate score
 
442
                ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
 
443
                rs   = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1);
399
444
                score = (int*)calloc(n_types * n, sizeof(int));
400
445
                pscore = (int*)calloc(n_types * n, sizeof(int));
401
446
                for (i = 0; i < n_types; ++i) {
 
447
                        ka_param_t ap = ka_param_blast;
 
448
                        ap.band_width = 2 * types[n_types - 1] + 2;
402
449
                        // write ref2
403
450
                        for (k = 0, j = left; j <= pos; ++j)
404
 
                                ref2[k++] = bam_nt16_table[(int)ref[j]];
 
451
                                ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
405
452
                        if (types[i] <= 0) j += -types[i];
406
453
                        else for (l = 0; l < types[i]; ++l)
407
 
                                         ref2[k++] = inscns[i*max_ins + l];
 
454
                                         ref2[k++] = bam_nt16_nt4_table[(int)inscns[i*max_ins + l]];
408
455
                        for (; j < right && ref[j]; ++j)
409
 
                                ref2[k++] = bam_nt16_table[(int)ref[j]];
 
456
                                ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
 
457
                        if (j < right) right = j;
410
458
                        // calculate score for each read
411
459
                        for (j = 0; j < n; ++j) {
412
460
                                const bam_pileup1_t *p = pl + j;
413
 
                                uint32_t *cigar;
414
 
                                bam1_core_t *c = &p->b->core;
415
 
                                int s, ps;
416
 
                                bam_segreg_t seg;
417
 
                                if (c->flag&BAM_FUNMAP) continue;
418
 
                                cigar = bam1_cigar(p->b);
419
 
                                bam_segreg(pos, c, cigar, &seg);
420
 
                                for (ps = s = 0, l = seg.qbeg; c->pos + l < right && l < seg.qend; ++l) {
421
 
                                        int cq = bam1_seqi(bam1_seq(p->b), l), ct;
422
 
                                        // in the following line, "<" will happen if reads are too long
423
 
                                        ct = c->pos + l - seg.qbeg >= left? ref2[c->pos + l - seg.qbeg - left] : 15;
424
 
                                        if (cq < 15 && ct < 15) {
425
 
                                                s += cq == ct? 1 : -mi->mm_penalty;
426
 
                                                if (cq != ct) ps += bam1_qual(p->b)[l];
 
461
                                int qbeg, qend, tbeg, tend;
 
462
                                if (p->b->core.flag & BAM_FUNMAP) continue;
 
463
                                qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left,  0, &tbeg);
 
464
                                qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend);
 
465
                                assert(tbeg >= left);
 
466
                                for (l = qbeg; l < qend; ++l)
 
467
                                        rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
 
468
                                {
 
469
                                        int x, y, n_acigar, ps;
 
470
                                        uint32_t *acigar;
 
471
                                        ps = 0;
 
472
                                        if (tend - tbeg + types[i] <= 0) {
 
473
                                                score[i*n+j] = -(1<<20);
 
474
                                                pscore[i*n+j] = 1<<20;
 
475
                                                continue;
427
476
                                        }
428
 
                                }
429
 
                                score[i*n + j] = s; pscore[i*n + j] = ps;
430
 
                                if (types[i] != 0) { // then try the other way to calculate the score
431
 
                                        for (ps = s = 0, l = seg.qbeg; c->pos + l + types[i] < right && l < seg.qend; ++l) {
432
 
                                                int cq = bam1_seqi(bam1_seq(p->b), l), ct;
433
 
                                                ct = c->pos + l - seg.qbeg + types[i] >= left? ref2[c->pos + l - seg.qbeg + types[i] - left] : 15;
434
 
                                                if (cq < 15 && ct < 15) {
435
 
                                                        s += cq == ct? 1 : -mi->mm_penalty;
436
 
                                                        if (cq != ct) ps += bam1_qual(p->b)[l];
 
477
                                        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);
 
478
                                        x = tbeg - left; y = 0;
 
479
                                        for (l = 0; l < n_acigar; ++l) {
 
480
                                                int op = acigar[l]&0xf;
 
481
                                                int len = acigar[l]>>4;
 
482
                                                if (op == BAM_CMATCH) {
 
483
                                                        int k;
 
484
                                                        for (k = 0; k < len; ++k)
 
485
                                                                if (ref2[x+k] != rs[y+k]) ps += bam1_qual(p->b)[y+k];
 
486
                                                        x += len; y += len;
 
487
                                                } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
 
488
                                                        if (op == BAM_CINS) ps += mi->q_indel * len;
 
489
                                                        y += len;
 
490
                                                } else if (op == BAM_CDEL) {
 
491
                                                        ps += mi->q_indel * len;
 
492
                                                        x += len;
437
493
                                                }
438
494
                                        }
 
495
                                        pscore[i*n+j] = ps;
 
496
                                        /*if (pos == 2618517) { // for debugging only
 
497
                                                fprintf(stderr, "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);
 
498
                                                for (l = 0; l < n_acigar; ++l) fprintf(stderr, "%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]); fprintf(stderr, "\n");
 
499
                                                for (l = 0; l < tend - tbeg + types[i]; ++l) fputc("ACGTN"[ref2[l]], stderr); fputc('\n', stderr);
 
500
                                                for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[rs[l]], stderr); fputc('\n', stderr);
 
501
                                                }*/
 
502
                                        free(acigar);
439
503
                                }
440
 
                                if (score[i*n+j] < s) score[i*n+j] = s; // choose the higher of the two scores
441
 
                                if (pscore[i*n+j] > ps) pscore[i*n+j] = ps;
442
 
                                //if (types[i] != 0) score[i*n+j] -= mi->indel_err;
443
 
                                //printf("%d, %d, %d, %d, %d, %d, %d\n", p->b->core.pos + 1, seg.qbeg, i, types[i], j,
444
 
                                //         score[i*n+j], pscore[i*n+j]);
445
504
                        }
446
505
                }
447
506
                { // get final result
491
550
                                else if (p->indel == ret->indel2) ++ret->cnt2;
492
551
                                else ++ret->cnt_anti;
493
552
                        }
494
 
                        // write gl[]
495
 
                        ret->gl[0] = ret->gl[1] = 0;
496
 
                        for (j = 0; j < n; ++j) {
497
 
                                int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
498
 
                                //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2);
499
 
                                if (s1 > s2) ret->gl[0] += s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel;
500
 
                                else ret->gl[1] += s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel;
 
553
                        { // write gl[]
 
554
                                int tmp, seq_err = 0;
 
555
                                double x = 1.0;
 
556
                                tmp = max1_i - max2_i;
 
557
                                if (tmp < 0) tmp = -tmp;
 
558
                                for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP;
 
559
                                seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP);
 
560
                                ret->gl[0] = ret->gl[1] = 0;
 
561
                                for (j = 0; j < n; ++j) {
 
562
                                        int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
 
563
                                        //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2);
 
564
                                        if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err;
 
565
                                        else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err;
 
566
                                }
501
567
                        }
502
568
                        // write cnt_ref and cnt_ambi
503
569
                        if (max1_i != 0 && max2_i != 0) {
509
575
                                }
510
576
                        }
511
577
                }
512
 
                free(score); free(pscore); free(ref2); free(inscns);
 
578
                free(score); free(pscore); free(ref2); free(rs); free(inscns);
513
579
        }
514
580
        { // call genotype
515
581
                int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5);