bam_maqcns.c

   1 #include <math.h>
   2 #include <assert.h>
   3 #include "bam.h"
   4 #include "bam_maqcns.h"
   5 #include "ksort.h"
   6 #include "errmod.h"
   7 #include "kaln.h"
   8 KSORT_INIT_GENERIC(uint32_t)
   9
  10 #define INDEL_WINDOW_SIZE 50
  11 #define INDEL_EXT_DEP 0.9
  12
  13 typedef struct __bmc_aux_t {
  14         int max;
  15         uint32_t *info;
  16         uint16_t *info16;
  17         errmod_t *em;
  18 } bmc_aux_t;
  19
  20 typedef struct {
  21         float esum[4], fsum[4];
  22         uint32_t c[4];
  23 } glf_call_aux_t;
  24
  25 char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
  26
  27 /*
  28   P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
  29   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]
  30   p_k = 1/k / \sum_{i=1}^{N-1} 1/i
  31  */
  32 static void cal_het(bam_maqcns_t *aa)
  33 {
  34         int k, n1, n2;
  35         double sum_harmo; // harmonic sum
  36         double poly_rate;
  37
  38         free(aa->lhet);
  39         aa->lhet = (double*)calloc(256 * 256, sizeof(double));
  40         sum_harmo = 0.0;
  41         for (k = 1; k <= aa->n_hap - 1; ++k)
  42                 sum_harmo += 1.0 / k;
  43         for (n1 = 0; n1 < 256; ++n1) {
  44                 for (n2 = 0; n2 < 256; ++n2) {
  45                         long double sum = 0.0;
  46                         double lC = aa->errmod == BAM_ERRMOD_SOAP? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1);
  47                         for (k = 1; k <= aa->n_hap - 1; ++k) {
  48                                 double pk = 1.0 / k / sum_harmo;
  49                                 double log1 = log((double)k/aa->n_hap);
  50                                 double log2 = log(1.0 - (double)k/aa->n_hap);
  51                                 sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
  52                         }
  53                         aa->lhet[n1<<8|n2] = lC + logl(sum);
  54                 }
  55         }
  56         poly_rate = aa->het_rate * sum_harmo;
  57         aa->q_r = -4.343 * log(2.0 * poly_rate / (1.0 - poly_rate));
  58 }
  59
  60 /** initialize the helper structure */
  61 static void cal_coef(bam_maqcns_t *aa)
  62 {
  63         int k, n, q;
  64         long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
  65         double *lC;
  66
  67         if (aa->errmod == BAM_ERRMOD_MAQ2) return; // no need to do the following
  68         // aa->lhet will be allocated and initialized
  69         free(aa->fk); free(aa->coef);
  70         aa->coef = 0;
  71         aa->fk = (double*)calloc(256, sizeof(double));
  72         aa->fk[0] = fk2[0] = 1.0;
  73         for (n = 1; n != 256; ++n) {
  74                 aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
  75                 fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
  76         }
  77         if (aa->errmod == BAM_ERRMOD_SOAP) return;
  78         aa->coef = (double*)calloc(256*256*64, sizeof(double));
  79         lC = (double*)calloc(256 * 256, sizeof(double));
  80         for (n = 1; n != 256; ++n)
  81                 for (k = 1; k <= n; ++k)
  82                         lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
  83         for (q = 1; q != 64; ++q) {
  84                 double e = pow(10.0, -q/10.0);
  85                 double le = log(e);
  86                 double le1 = log(1.0-e);
  87                 for (n = 1; n != 256; ++n) {
  88                         double *coef = aa->coef + (q<<16|n<<8);
  89                         sum_a[n+1] = 0.0;
  90                         for (k = n; k >= 0; --k) { // a_k = \sum_{i=k}^n C^n_k \epsilon^k (1-\epsilon)^{n-k}
  91                                 sum_a[k] = sum_a[k+1] + expl(lC[n<<8|k] + k*le + (n-k)*le1);
  92                                 b[k] = sum_a[k+1] / sum_a[k];
  93                                 if (b[k] > 0.99) b[k] = 0.99;
  94                         }
  95                         for (k = 0; k != n; ++k) // log(\bar\beta_{nk}(\bar\epsilon)^{f_k})
  96                                 q_c[k] = -4.343 * fk2[k] * logl(b[k] / e);
  97                         for (k = 1; k != n; ++k) q_c[k] += q_c[k-1]; // \prod_{i=0}^k c_i
  98                         for (k = 0; k <= n; ++k) { // powl() in 64-bit mode seems broken on my Mac OS X 10.4.9
  99                                 tmp[k] = -4.343 * logl(1.0 - expl(fk2[k] * logl(b[k])));
 100                                 coef[k] = (k? q_c[k-1] : 0) + tmp[k]; // this is the final c_{nk}
 101                         }
 102                 }
 103         }
 104         free(lC);
 105 }
 106
 107 bam_maqcns_t *bam_maqcns_init()
 108 {
 109         bam_maqcns_t *bm;
 110         bm = (bam_maqcns_t*)calloc(1, sizeof(bam_maqcns_t));
 111         bm->aux = (bmc_aux_t*)calloc(1, sizeof(bmc_aux_t));
 112         bm->het_rate = 0.001;
 113         bm->theta = 0.83f;
 114         bm->n_hap = 2;
 115         bm->eta = 0.03;
 116         bm->cap_mapQ = 60;
 117         return bm;
 118 }
 119
 120 void bam_maqcns_prepare(bam_maqcns_t *bm)
 121 {
 122         if (bm->errmod == BAM_ERRMOD_MAQ2) bm->aux->em = errmod_init(1. - bm->theta);
 123         cal_coef(bm); cal_het(bm);
 124 }
 125
 126 void bam_maqcns_destroy(bam_maqcns_t *bm)
 127 {
 128         if (bm == 0) return;
 129         free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info); free(bm->aux->info16);
 130         if (bm->aux->em) errmod_destroy(bm->aux->em);
 131         free(bm->aux); free(bm);
 132 }
 133
 134 glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
 135 {
 136         glf_call_aux_t *b = 0;
 137         int i, j, k, w[8], c, n;
 138         glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
 139         float p[16], min_p = 1e30;
 140         uint64_t rms;
 141
 142         g->ref_base = ref_base;
 143         if (_n == 0) return g;
 144
 145         // construct aux array
 146         if (bm->aux->max < _n) {
 147                 bm->aux->max = _n;
 148                 kroundup32(bm->aux->max);
 149                 bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
 150                 bm->aux->info16 = (uint16_t*)realloc(bm->aux->info16, 2 * bm->aux->max);
 151         }
 152         for (i = n = 0, rms = 0; i < _n; ++i) {
 153                 const bam_pileup1_t *p = pl + i;
 154                 uint32_t q, x = 0, qq;
 155                 uint16_t y = 0;
 156                 if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue;
 157                 q = (uint32_t)bam1_qual(p->b)[p->qpos];
 158                 x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
 159                 y |= bam1_strand(p->b)<<4;
 160                 if (p->b->core.qual < q) q = p->b->core.qual;
 161                 c = p->b->core.qual < bm->cap_mapQ? p->b->core.qual : bm->cap_mapQ;
 162                 rms += c * c;
 163                 x |= q << 24;
 164                 y |= q << 5;
 165                 qq = bam1_seqi(bam1_seq(p->b), p->qpos);
 166                 q = bam_nt16_nt4_table[qq? qq : ref_base];
 167                 if (!p->is_del && q < 4) x |= 1 << 21 | q << 16, y |= q;
 168                 bm->aux->info16[n] = y;
 169                 bm->aux->info[n++] = x;
 170         }
 171         rms = (uint8_t)(sqrt((double)rms / n) + .499);
 172         if (bm->errmod == BAM_ERRMOD_MAQ2) {
 173                 errmod_cal(bm->aux->em, n, 4, bm->aux->info16, p);
 174                 goto goto_glf;
 175         }
 176         ks_introsort(uint32_t, n, bm->aux->info);
 177         // generate esum and fsum
 178         b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
 179         for (k = 0; k != 8; ++k) w[k] = 0;
 180         for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
 181                 uint32_t info = bm->aux->info[j];
 182                 if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
 183                 k = info>>16&7;
 184                 if (info>>24 > 0) {
 185                         b->esum[k&3] += bm->fk[w[k]] * (info>>24);
 186                         b->fsum[k&3] += bm->fk[w[k]];
 187                         if (w[k] < 0xff) ++w[k];
 188                         ++b->c[k&3];
 189                 }
 190         }
 191         // rescale ->c[]
 192         for (j = c = 0; j != 4; ++j) c += b->c[j];
 193         if (c > 255) {
 194                 for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
 195                 for (j = c = 0; j != 4; ++j) c += b->c[j];
 196         }
 197         if (bm->errmod == BAM_ERRMOD_MAQ) {
 198                 // generate likelihood
 199                 for (j = 0; j != 4; ++j) {
 200                         // homozygous
 201                         float tmp1, tmp3;
 202                         int tmp2, bar_e;
 203                         for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
 204                                 if (j == k) continue;
 205                                 tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
 206                         }
 207                         if (tmp2) {
 208                                 bar_e = (int)(tmp1 / tmp3 + 0.5);
 209                                 if (bar_e < 4) bar_e = 4; // should not happen
 210                                 if (bar_e > 63) bar_e = 63;
 211                                 p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
 212                         } else p[j<<2|j] = 0.0; // all the bases are j
 213                         // heterozygous
 214                         for (k = j + 1; k < 4; ++k) {
 215                                 for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
 216                                         if (i == j || i == k) continue;
 217                                         tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
 218                                 }
 219                                 if (tmp2) {
 220                                         bar_e = (int)(tmp1 / tmp3 + 0.5);
 221                                         if (bar_e < 4) bar_e = 4;
 222                                         if (bar_e > 63) bar_e = 63;
 223                                         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];
 224                                 } 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
 225                         }
 226                         //
 227                         for (k = 0; k != 4; ++k)
 228                                 if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
 229                 }
 230
 231                 { // fix p[k<<2|k]
 232                         float max1, max2, min1, min2;
 233                         int max_k, min_k;
 234                         max_k = min_k = -1;
 235                         max1 = max2 = -1.0; min1 = min2 = 1e30;
 236                         for (k = 0; k < 4; ++k) {
 237                                 if (b->esum[k] > max1) {
 238                                         max2 = max1; max1 = b->esum[k]; max_k = k;
 239                                 } else if (b->esum[k] > max2) max2 = b->esum[k];
 240                         }
 241                         for (k = 0; k < 4; ++k) {
 242                                 if (p[k<<2|k] < min1) {
 243                                         min2 = min1; min1 = p[k<<2|k]; min_k = k;
 244                                 } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
 245                         }
 246                         if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
 247                                 p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
 248                 }
 249         } else if (bm->errmod == BAM_ERRMOD_SOAP) { // apply the SOAP model
 250                 // generate likelihood
 251                 for (j = 0; j != 4; ++j) {
 252                         float tmp;
 253                         // homozygous
 254                         for (k = 0, tmp = 0.0; k != 4; ++k)
 255                                 if (j != k) tmp += b->esum[k];
 256                         p[j<<2|j] = tmp;
 257                         // heterozygous
 258                         for (k = j + 1; k < 4; ++k) {
 259                                 for (i = 0, tmp = 0.0; i != 4; ++i)
 260                                         if (i != j && i != k) tmp += b->esum[i];
 261                                 p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp;
 262                         }
 263                 }
 264         }
 265
 266 goto_glf:
 267         // convert necessary information to glf1_t
 268         g->ref_base = ref_base; g->max_mapQ = rms;
 269         g->depth = n > 16777215? 16777215 : n;
 270         for (j = 0; j != 4; ++j)
 271                 for (k = j; k < 4; ++k)
 272                         if (p[j<<2|k] < min_p) min_p = p[j<<2|k];
 273         g->min_lk = min_p > 255.0? 255 : (int)(min_p + 0.5);
 274         for (j = c = 0; j != 4; ++j)
 275                 for (k = j; k < 4; ++k)
 276                         g->lk[c++] = p[j<<2|k]-min_p > 255.0? 255 : (int)(p[j<<2|k]-min_p + 0.5);
 277
 278         free(b);
 279         return g;
 280 }
 281
 282 uint32_t glf2cns(const glf1_t *g, int q_r)
 283 {
 284         int i, j, k, tmp[16], min = 10000, min2 = 10000, min3 = 10000, min_g = -1, min_g2 = -1;
 285         uint32_t x = 0;
 286         for (i = k = 0; i < 4; ++i)
 287                 for (j = i; j < 4; ++j) {
 288                         tmp[j<<2|i] = -1;
 289                         tmp[i<<2|j] = g->lk[k++] + (i == j? 0 : q_r);
 290                 }
 291         for (i = 0; i < 16; ++i) {
 292                 if (tmp[i] < 0) continue;
 293                 if (tmp[i] < min) {
 294                         min3 = min2; min2 = min; min = tmp[i]; min_g2 = min_g; min_g = i;
 295                 } else if (tmp[i] < min2) {
 296                         min3 = min2; min2 = tmp[i]; min_g2 = i;
 297                 } else if (tmp[i] < min3) min3 = tmp[i];
 298         }
 299         x = min_g >= 0? (1U<<(min_g>>2&3) | 1U<<(min_g&3)) << 28 : 0xf << 28;
 300         x |= min_g2 >= 0? (1U<<(min_g2>>2&3) | 1U<<(min_g2&3)) << 24 : 0xf << 24;
 301         x |= (uint32_t)g->max_mapQ << 16;
 302         x |= min2 < 10000? (min2 - min < 256? min2 - min : 255) << 8 : 0xff << 8;
 303         x |= min2 < 10000 && min3 < 10000? (min3 - min2 < 256? min3 - min2 : 255) : 0xff;
 304         return x;
 305 }
 306
 307 uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm)
 308 {
 309         glf1_t *g;
 310         uint32_t x;
 311         if (n) {
 312                 g = bam_maqcns_glfgen(n, pl, 0xf, bm);
 313                 x = glf2cns(g, (int)(bm->q_r + 0.5));
 314                 free(g);
 315         } else x = 0xfU<<28 | 0xfU<<24;
 316         return x;
 317 }
 318
 319 /************** *****************/
 320
 321 bam_maqindel_opt_t *bam_maqindel_opt_init()
 322 {
 323         bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
 324         mi->q_indel = 40;
 325         mi->r_indel = 0.00015;
 326         mi->r_snp = 0.001;
 327         //
 328         mi->mm_penalty = 3;
 329         mi->indel_err = 4;
 330         mi->ambi_thres = 10;
 331         return mi;
 332 }
 333
 334 void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
 335 {
 336         if (mir == 0) return;
 337         free(mir->s[0]); free(mir->s[1]); free(mir);
 338 }
 339
 340 int bam_tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
 341 {
 342         int k, x = c->pos, y = 0, last_y = 0;
 343         *_tpos = c->pos;
 344         for (k = 0; k < c->n_cigar; ++k) {
 345                 int op = cigar[k] & BAM_CIGAR_MASK;
 346                 int l = cigar[k] >> BAM_CIGAR_SHIFT;
 347                 if (op == BAM_CMATCH) {
 348                         if (c->pos > tpos) return y;
 349                         if (x + l > tpos) {
 350                                 *_tpos = tpos;
 351                                 return y + (tpos - x);
 352                         }
 353                         x += l; y += l;
 354                         last_y = y;
 355                 } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
 356                 else if (op == BAM_CDEL || op == BAM_CREF_SKIP) {
 357                         if (x + l > tpos) {
 358                                 *_tpos = is_left? x : x + l;
 359                                 return y;
 360                         }
 361                         x += l;
 362                 }
 363         }
 364         *_tpos = x;
 365         return last_y;
 366 }
 367
 368 #define MINUS_CONST 0x10000000
 369
 370 bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
 371                                                                  int _n_types, int *_types)
 372 {
 373         int i, j, n_types, *types, left, right, max_rd_len = 0;
 374         bam_maqindel_ret_t *ret = 0;
 375         // if there is no proposed indel, check if there is an indel from the alignment
 376         if (_n_types == 0) {
 377                 for (i = 0; i < n; ++i) {
 378                         const bam_pileup1_t *p = pl + i;
 379                         if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
 380                 }
 381                 if (i == n) return 0; // no indel
 382         }
 383         { // calculate how many types of indels are available (set n_types and types)
 384                 int m;
 385                 uint32_t *aux;
 386                 aux = (uint32_t*)calloc(n + _n_types + 1, 4);
 387                 m = 0;
 388                 aux[m++] = MINUS_CONST; // zero indel is always a type
 389                 for (i = 0; i < n; ++i) {
 390                         const bam_pileup1_t *p = pl + i;
 391                         if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
 392                                 aux[m++] = MINUS_CONST + p->indel;
 393                         j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
 394                         if (j > max_rd_len) max_rd_len = j;
 395                 }
 396                 if (_n_types) // then also add this to aux[]
 397                         for (i = 0; i < _n_types; ++i)
 398                                 if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
 399                 ks_introsort(uint32_t, m, aux);
 400                 // squeeze out identical types
 401                 for (i = 1, n_types = 1; i < m; ++i)
 402                         if (aux[i] != aux[i-1]) ++n_types;
 403                 types = (int*)calloc(n_types, sizeof(int));
 404                 j = 0;
 405                 types[j++] = aux[0] - MINUS_CONST;
 406                 for (i = 1; i < m; ++i) {
 407                         if (aux[i] != aux[i-1])
 408                                 types[j++] = aux[i] - MINUS_CONST;
 409                 }
 410                 free(aux);
 411         }
 412         { // calculate left and right boundary
 413                 left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
 414                 right = pos + INDEL_WINDOW_SIZE;
 415                 if (types[0] < 0) right -= types[0];
 416                 // in case the alignments stand out the reference
 417                 for (i = pos; i < right; ++i)
 418                         if (ref[i] == 0) break;
 419                 right = i;
 420         }
 421         { // the core part
 422                 char *ref2, *rs, *inscns = 0;
 423                 int qr_snp, k, l, *score, *pscore, max_ins = types[n_types-1];
 424                 qr_snp = (int)(-4.343 * log(mi->r_snp) + .499);
 425                 if (max_ins > 0) { // get the consensus of inserted sequences
 426                         int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
 427                         // count occurrences
 428                         for (i = 0; i < n_types; ++i) {
 429                                 if (types[i] <= 0) continue; // not insertion
 430                                 for (j = 0; j < n; ++j) {
 431                                         const bam_pileup1_t *p = pl + j;
 432                                         if (!(p->b->core.flag&BAM_FUNMAP) && p->indel == types[i]) {
 433                                                 for (k = 1; k <= p->indel; ++k) {
 434                                                         int c = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos + k)];
 435                                                         if (c < 4) ++inscns_aux[i*max_ins*4 + (k-1)*4 + c];
 436                                                 }
 437                                         }
 438                                 }
 439                         }
 440                         // construct the consensus of inserted sequence
 441                         inscns = (char*)calloc(n_types * max_ins, sizeof(char));
 442                         for (i = 0; i < n_types; ++i) {
 443                                 for (j = 0; j < types[i]; ++j) {
 444                                         int max = 0, max_k = -1, *ia = inscns_aux + i*max_ins*4 + j*4;
 445                                         for (k = 0; k < 4; ++k) {
 446                                                 if (ia[k] > max) {
 447                                                         max = ia[k];
 448                                                         max_k = k;
 449                                                 }
 450                                         }
 451                                         inscns[i*max_ins + j] = max? 1<<max_k : 15;
 452                                 }
 453                         }
 454                         free(inscns_aux);
 455                 }
 456                 // calculate score
 457                 ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
 458                 rs   = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1);
 459                 score = (int*)calloc(n_types * n, sizeof(int));
 460                 pscore = (int*)calloc(n_types * n, sizeof(int));
 461                 for (i = 0; i < n_types; ++i) {
 462                         ka_param_t ap = ka_param_blast;
 463                         ap.band_width = 2 * types[n_types - 1] + 2;
 464                         ap.gap_end_ext = 0;
 465                         // write ref2
 466                         for (k = 0, j = left; j <= pos; ++j)
 467                                 ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
 468                         if (types[i] <= 0) j += -types[i];
 469                         else for (l = 0; l < types[i]; ++l)
 470                                          ref2[k++] = bam_nt16_nt4_table[(int)inscns[i*max_ins + l]];
 471                         if (types[0] < 0) { // mask deleted sequences
 472                                 int jj, tmp = types[i] >= 0? -types[0] : -types[0] + types[i];
 473                                 for (jj = 0; jj < tmp && j < right && ref[j]; ++jj, ++j)
 474                                         ref2[k++] = 4;
 475                         }
 476                         for (; j < right && ref[j]; ++j)
 477                                 ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
 478                         if (j < right) right = j;
 479                         // calculate score for each read
 480                         for (j = 0; j < n; ++j) {
 481                                 const bam_pileup1_t *p = pl + j;
 482                                 int qbeg, qend, tbeg, tend;
 483                                 if (p->b->core.flag & BAM_FUNMAP) continue;
 484                                 qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left,  0, &tbeg);
 485                                 qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend);
 486                                 assert(tbeg >= left);
 487                                 for (l = qbeg; l < qend; ++l)
 488                                         rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
 489                                 {
 490                                         int x, y, n_acigar, ps;
 491                                         uint32_t *acigar;
 492                                         ps = 0;
 493                                         if (tend - tbeg + types[i] <= 0) {
 494                                                 score[i*n+j] = -(1<<20);
 495                                                 pscore[i*n+j] = 1<<20;
 496                                                 continue;
 497                                         }
 498                                         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);
 499                                         x = tbeg - left; y = 0;
 500                                         for (l = 0; l < n_acigar; ++l) {
 501                                                 int op = acigar[l]&0xf;
 502                                                 int len = acigar[l]>>4;
 503                                                 if (op == BAM_CMATCH) {
 504                                                         int k;
 505                                                         for (k = 0; k < len; ++k)
 506                                                                 if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4)
 507                                                                         ps += bam1_qual(p->b)[y+k] < qr_snp? bam1_qual(p->b)[y+k] : qr_snp;
 508                                                         x += len; y += len;
 509                                                 } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
 510                                                         if (op == BAM_CINS && l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
 511                                                         y += len;
 512                                                 } else if (op == BAM_CDEL) {
 513                                                         if (l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
 514                                                         x += len;
 515                                                 }
 516                                         }
 517                                         pscore[i*n+j] = ps;
 518                                         /*if (1) { // for debugging only
 519                                                 fprintf(stderr, "id=%d, pos=%d, type=%d, j=%d, score=%d, psore=%d, %d, %d, %d, %d, %d, ",
 520                                                                 j, pos+1, types[i], j, score[i*n+j], pscore[i*n+j], tbeg, tend, qbeg, qend, mi->q_indel);
 521                                                 for (l = 0; l < n_acigar; ++l) fprintf(stderr, "%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]);
 522                                                 fprintf(stderr, "\n");
 523                                                 for (l = 0; l < tend - tbeg + types[i]; ++l) fputc("ACGTN"[ref2[l+tbeg-left]], stderr);
 524                                                 fputc('\n', stderr);
 525                                                 for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[rs[l]], stderr);
 526                                                 fputc('\n', stderr);
 527                                                 }*/
 528                                         free(acigar);
 529                                 }
 530                         }
 531                 }
 532                 { // get final result
 533                         int *sum, max1, max2, max1_i, max2_i;
 534                         // pick up the best two score
 535                         sum = (int*)calloc(n_types, sizeof(int));
 536                         for (i = 0; i < n_types; ++i)
 537                                 for (j = 0; j < n; ++j)
 538                                         sum[i] += -pscore[i*n+j];
 539                         max1 = max2 = -0x7fffffff; max1_i = max2_i = -1;
 540                         for (i = 0; i < n_types; ++i) {
 541                                 if (sum[i] > max1) {
 542                                         max2 = max1; max2_i = max1_i; max1 = sum[i]; max1_i = i;
 543                                 } else if (sum[i] > max2) {
 544                                         max2 = sum[i]; max2_i = i;
 545                                 }
 546                         }
 547                         free(sum);
 548                         // write ret
 549                         ret = (bam_maqindel_ret_t*)calloc(1, sizeof(bam_maqindel_ret_t));
 550                         ret->indel1 = types[max1_i]; ret->indel2 = types[max2_i];
 551                         ret->s[0] = (char*)calloc(abs(ret->indel1) + 2, 1);
 552                         ret->s[1] = (char*)calloc(abs(ret->indel2) + 2, 1);
 553                         // write indel sequence
 554                         if (ret->indel1 > 0) {
 555                                 ret->s[0][0] = '+';
 556                                 for (k = 0; k < ret->indel1; ++k)
 557                                         ret->s[0][k+1] = bam_nt16_rev_table[(int)inscns[max1_i*max_ins + k]];
 558                         } else if (ret->indel1 < 0) {
 559                                 ret->s[0][0] = '-';
 560                                 for (k = 0; k < -ret->indel1 && ref[pos + k + 1]; ++k)
 561                                         ret->s[0][k+1] = ref[pos + k + 1];
 562                         } else ret->s[0][0] = '*';
 563                         if (ret->indel2 > 0) {
 564                                 ret->s[1][0] = '+';
 565                                 for (k = 0; k < ret->indel2; ++k)
 566                                         ret->s[1][k+1] = bam_nt16_rev_table[(int)inscns[max2_i*max_ins + k]];
 567                         } else if (ret->indel2 < 0) {
 568                                 ret->s[1][0] = '-';
 569                                 for (k = 0; k < -ret->indel2 && ref[pos + k + 1]; ++k)
 570                                         ret->s[1][k+1] = ref[pos + k + 1];
 571                         } else ret->s[1][0] = '*';
 572                         // write count
 573                         for (i = 0; i < n; ++i) {
 574                                 const bam_pileup1_t *p = pl + i;
 575                                 if (p->indel == ret->indel1) ++ret->cnt1;
 576                                 else if (p->indel == ret->indel2) ++ret->cnt2;
 577                                 else ++ret->cnt_anti;
 578                         }
 579                         { // write gl[]
 580                                 int tmp, seq_err = 0;
 581                                 double x = 1.0;
 582                                 tmp = max1_i - max2_i;
 583                                 if (tmp < 0) tmp = -tmp;
 584                                 for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP;
 585                                 seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP);
 586                                 ret->gl[0] = ret->gl[1] = 0;
 587                                 for (j = 0; j < n; ++j) {
 588                                         int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
 589                                         //fprintf(stderr, "id=%d, %d, %d, %d, %d, %d\n", j, pl[j].b->core.pos+1, types[max1_i], types[max2_i], s1, s2);
 590                                         if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err;
 591                                         else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err;
 592                                 }
 593                         }
 594                         // write cnt_ref and cnt_ambi
 595                         if (max1_i != 0 && max2_i != 0) {
 596                                 for (j = 0; j < n; ++j) {
 597                                         int diff1 = score[j] - score[max1_i * n + j];
 598                                         int diff2 = score[j] - score[max2_i * n + j];
 599                                         if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
 600                                         else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
 601                                 }
 602                         }
 603                 }
 604                 free(score); free(pscore); free(ref2); free(rs); free(inscns);
 605         }
 606         { // call genotype
 607                 int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5);
 608                 int min1, min2, min1_i;
 609                 q[0] = ret->gl[0] + (ret->s[0][0] != '*'? 0 : 0) * qr_indel;
 610                 q[1] = ret->gl[1] + (ret->s[1][0] != '*'? 0 : 0) * qr_indel;
 611                 q[2] = n * 3 + (ret->s[0][0] == '*' || ret->s[1][0] == '*'? 1 : 1) * qr_indel;
 612                 min1 = min2 = 0x7fffffff; min1_i = -1;
 613                 for (i = 0; i < 3; ++i) {
 614                         if (q[i] < min1) {
 615                                 min2 = min1; min1 = q[i]; min1_i = i;
 616                         } else if (q[i] < min2) min2 = q[i];
 617                 }
 618                 ret->gt = min1_i;
 619                 ret->q_cns = min2 - min1;
 620                 // set q_ref
 621                 if (ret->gt < 2) ret->q_ref = (ret->s[ret->gt][0] == '*')? 0 : q[1-ret->gt] - q[ret->gt] - qr_indel - 3;
 622                 else ret->q_ref = (ret->s[0][0] == '*')? q[0] - q[2] : q[1] - q[2];
 623                 if (ret->q_ref < 0) ret->q_ref = 0;
 624         }
 625         free(types);
 626         return ret;
 627 }