bam_maqcns.c

   1 #include <math.h>
   2 #include "bam.h"
   3 #include "bam_maqcns.h"
   4 #include "ksort.h"
   5 KSORT_INIT_GENERIC(uint32_t)
   6
   7 typedef struct __bmc_aux_t {
   8         int max;
   9         uint32_t *info;
  10 } bmc_aux_t;
  11
  12 typedef struct {
  13         float esum[4], fsum[4];
  14         uint32_t c[4];
  15         uint32_t mapQ_max;
  16 } glf_call_aux_t;
  17
  18 char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
  19
  20 /*
  21   P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
  22   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]
  23   p_k = i/k / \sum_{i=1}^{N-1} 1/i
  24  */
  25 static void cal_het(bam_maqcns_t *aa)
  26 {
  27         int k, n1, n2;
  28         double sum_harmo; // harmonic sum
  29         double poly_rate;
  30         double p1 = 0.0, p3 = 0.0; // just for testing
  31
  32         free(aa->lhet);
  33         aa->lhet = (double*)calloc(256 * 256, sizeof(double));
  34         sum_harmo = 0.0;
  35         for (k = 1; k <= aa->n_hap - 1; ++k)
  36                 sum_harmo += 1.0 / k;
  37         for (n1 = 0; n1 < 256; ++n1) {
  38                 for (n2 = 0; n2 < 256; ++n2) {
  39                         long double sum = 0.0;
  40                         double lC = lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
  41                         for (k = 1; k <= aa->n_hap - 1; ++k) {
  42                                 double pk = 1.0 / k / sum_harmo;
  43                                 double log1 = log((double)k/aa->n_hap);
  44                                 double log2 = log(1.0 - (double)k/aa->n_hap);
  45                                 sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
  46                         }
  47                         aa->lhet[n1<<8|n2] = lC + logl(sum);
  48                         if (n1 == 17 && n2 == 3) p3 = lC + logl(expl(logl(0.5) * 20));
  49                         if (n1 == 19 && n2 == 1) p1 = lC + logl(expl(logl(0.5) * 20));
  50                 }
  51         }
  52         poly_rate = aa->het_rate * sum_harmo;
  53         aa->q_r = -4.343 * log(2.0 * poly_rate / (1.0 - poly_rate));
  54 }
  55
  56 /** initialize the helper structure */
  57 static void cal_coef(bam_maqcns_t *aa)
  58 {
  59         int k, n, q;
  60         long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
  61         double *lC;
  62
  63         lC = (double*)calloc(256 * 256, sizeof(double));
  64         // aa->lhet will be allocated and initialized
  65         free(aa->fk); free(aa->coef);
  66         aa->fk = (double*)calloc(256, sizeof(double));
  67         aa->coef = (double*)calloc(256*256*64, sizeof(double));
  68         aa->fk[0] = fk2[0] = 1.0;
  69         for (n = 1; n != 256; ++n) {
  70                 aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
  71                 fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
  72         }
  73         for (n = 1; n != 256; ++n)
  74                 for (k = 1; k <= n; ++k)
  75                         lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
  76         for (q = 1; q != 64; ++q) {
  77                 double e = pow(10.0, -q/10.0);
  78                 double le = log(e);
  79                 double le1 = log(1.0-e);
  80                 for (n = 1; n != 256; ++n) {
  81                         double *coef = aa->coef + (q<<16|n<<8);
  82                         sum_a[n+1] = 0.0;
  83                         for (k = n; k >= 0; --k) { // a_k = \sum_{i=k}^n C^n_k \epsilon^k (1-\epsilon)^{n-k}
  84                                 sum_a[k] = sum_a[k+1] + expl(lC[n<<8|k] + k*le + (n-k)*le1);
  85                                 b[k] = sum_a[k+1] / sum_a[k];
  86                                 if (b[k] > 0.99) b[k] = 0.99;
  87                         }
  88                         for (k = 0; k != n; ++k) // log(\bar\beta_{nk}(\bar\epsilon)^{f_k})
  89                                 q_c[k] = -4.343 * fk2[k] * logl(b[k] / e);
  90                         for (k = 1; k != n; ++k) q_c[k] += q_c[k-1]; // \prod_{i=0}^k c_i
  91                         for (k = 0; k <= n; ++k) { // powl() in 64-bit mode seems broken on my Mac OS X 10.4.9
  92                                 tmp[k] = -4.343 * logl(1.0 - expl(fk2[k] * logl(b[k])));
  93                                 coef[k] = (k? q_c[k-1] : 0) + tmp[k]; // this is the final c_{nk}
  94                         }
  95                 }
  96         }
  97         free(lC);
  98 }
  99
 100 bam_maqcns_t *bam_maqcns_init()
 101 {
 102         bam_maqcns_t *bm;
 103         bm = (bam_maqcns_t*)calloc(1, sizeof(bam_maqcns_t));
 104         bm->aux = (bmc_aux_t*)calloc(1, sizeof(bmc_aux_t));
 105         bm->het_rate = 0.001;
 106         bm->theta = 0.85;
 107         bm->n_hap = 2;
 108         bm->eta = 0.03;
 109         return bm;
 110 }
 111
 112 void bam_maqcns_prepare(bam_maqcns_t *bm)
 113 {
 114         cal_coef(bm); cal_het(bm);
 115 }
 116
 117 void bam_maqcns_destroy(bam_maqcns_t *bm)
 118 {
 119         if (bm == 0) return;
 120         free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info);
 121         free(bm->aux); free(bm);
 122 }
 123
 124 glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
 125 {
 126         glf_call_aux_t *b;
 127         int i, j, k, w[8], c, n;
 128         glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
 129         float p[16], min_p = 1e30;
 130
 131         g->ref_base = ref_base;
 132         if (_n == 0) return g;
 133
 134         // construct aux array
 135         if (bm->aux->max < _n) {
 136                 bm->aux->max = _n;
 137                 kroundup32(bm->aux->max);
 138                 bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
 139         }
 140         for (i = n = 0; i < _n; ++i) {
 141                 const bam_pileup1_t *p = pl + i;
 142                 uint32_t q, x = 0;
 143                 if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue;
 144                 q = (uint32_t)bam1_qual(p->b)[p->qpos];
 145                 x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
 146                 if (p->b->core.qual < q) q = p->b->core.qual;
 147                 x |= q << 24;
 148                 q = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
 149                 if (!p->is_del && q < 4) x |= 1 << 21 | q << 16;
 150                 bm->aux->info[n++] = x;
 151         }
 152         ks_introsort(uint32_t, n, bm->aux->info);
 153         // generate esum and fsum
 154         b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
 155         for (k = 0; k != 8; ++k) w[k] = 0;
 156         b->mapQ_max = 0;
 157         for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
 158                 uint32_t info = bm->aux->info[j];
 159                 if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
 160                 k = info>>16&7;
 161                 if (info>>24 > 0) {
 162                         b->esum[k&3] += bm->fk[w[k]] * (info>>24);
 163                         b->fsum[k&3] += bm->fk[w[k]];
 164                         if (w[k] < 0xff) ++w[k];
 165                         ++b->c[k&3];
 166                 }
 167                 if (b->mapQ_max < (info&0x7f)) b->mapQ_max = info&0x7f;
 168         }
 169         // rescale ->c[]
 170         for (j = c = 0; j != 4; ++j) c += b->c[j];
 171         if (c > 255) {
 172                 for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
 173                 for (j = c = 0; j != 4; ++j) c += b->c[j];
 174         }
 175         // generate likelihood
 176         for (j = 0; j != 4; ++j) {
 177                 // homozygous
 178                 float tmp1, tmp3;
 179                 int tmp2, bar_e;
 180                 for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
 181                         if (j == k) continue;
 182                         tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
 183                 }
 184                 if (tmp2) {
 185                         bar_e = (int)(tmp1 / tmp3 + 0.5);
 186                         if (bar_e < 4) bar_e = 4; // should not happen
 187                         if (bar_e > 63) bar_e = 63;
 188                         p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
 189                 } else p[j<<2|j] = 0.0; // all the bases are j
 190                 // heterozygous
 191                 for (k = j + 1; k < 4; ++k) {
 192                         for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
 193                                 if (i == j || i == k) continue;
 194                                 tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
 195                         }
 196                         if (tmp2) {
 197                                 bar_e = (int)(tmp1 / tmp3 + 0.5);
 198                                 if (bar_e < 4) bar_e = 4;
 199                                 if (bar_e > 63) bar_e = 63;
 200                                 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];
 201                         } 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
 202                 }
 203                 //
 204                 for (k = 0; k != 4; ++k)
 205                         if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
 206         }
 207
 208         // convert necessary information to glf1_t
 209         g->ref_base = ref_base; g->max_mapQ = b->mapQ_max;
 210         g->depth = n > 16777215? 16777215 : n;
 211         for (j = 0; j != 4; ++j)
 212                 for (k = j; k < 4; ++k)
 213                         if (p[j<<2|k] < min_p) min_p = p[j<<2|k];
 214         g->min_lk = min_p > 255.0? 255 : (int)(min_p + 0.5);
 215         for (j = c = 0; j != 4; ++j)
 216                 for (k = j; k < 4; ++k)
 217                         g->lk[c++] = p[j<<2|k]-min_p > 255.0? 255 : (int)(p[j<<2|k]-min_p + 0.5);
 218
 219         free(b);
 220         return g;
 221 }
 222
 223 uint32_t glf2cns(const glf1_t *g, int q_r)
 224 {
 225         int i, j, k, tmp[16], min = 10000, min2 = 10000, min3 = 10000, min_g = -1, min_g2 = -1;
 226         uint32_t x = 0;
 227         for (i = k = 0; i < 4; ++i)
 228                 for (j = i; j < 4; ++j) {
 229                         tmp[j<<2|i] = -1;
 230                         tmp[i<<2|j] = g->lk[k++] + (i == j? 0 : q_r);
 231                 }
 232         for (i = 0; i < 16; ++i) {
 233                 if (tmp[i] < 0) continue;
 234                 if (tmp[i] < min) {
 235                         min3 = min2; min2 = min; min = tmp[i]; min_g2 = min_g; min_g = i;
 236                 } else if (tmp[i] < min2) {
 237                         min3 = min2; min2 = tmp[i]; min_g2 = i;
 238                 } else if (tmp[i] < min3) min3 = tmp[i];
 239         }
 240         x = min_g >= 0? (1U<<(min_g>>2&3) | 1U<<(min_g&3)) << 28 : 0xf << 28;
 241         x |= min_g2 >= 0? (1U<<(min_g2>>2&3) | 1U<<(min_g2&3)) << 24 : 0xf << 24;
 242         x |= (uint32_t)g->max_mapQ << 16;
 243         x |= min2 < 10000? (min2 - min < 256? min2 - min : 255) << 8 : 0xff << 8;
 244         x |= min2 < 10000 && min3 < 10000? (min3 - min2 < 256? min3 - min2 : 255) : 0xff;
 245         return x;
 246 }
 247
 248 uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm)
 249 {
 250         glf1_t *g;
 251         uint32_t x;
 252         if (n) {
 253                 g = bam_maqcns_glfgen(n, pl, 0xf, bm);
 254                 x = glf2cns(g, (int)(bm->q_r + 0.5));
 255                 free(g);
 256         } else x = 0xfU<<28 | 0xfU<<24;
 257         return x;
 258 }
 259
 260 /************** *****************/
 261
 262 bam_maqindel_opt_t *bam_maqindel_opt_init()
 263 {
 264         bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
 265         mi->mm_penalty = 3;
 266         mi->indel_err = 4;
 267         mi->ambi_thres = 10;
 268         return mi;
 269 }
 270
 271 void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
 272 {
 273         if (mir == 0) return;
 274         free(mir->s1); free(mir->s2); free(mir);
 275 }
 276
 277 #define MINUS_CONST 0x10000000
 278
 279 bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref)
 280 {
 281         int i, j, n_types, *types, left, right;
 282         bam_maqindel_ret_t *ret = 0;
 283         for (i = 0; i < n; ++i) {
 284                 const bam_pileup1_t *p = pl + i;
 285                 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
 286         }
 287         if (i == n) return 0; // no indel
 288         { // calculate how many types of indels are available (set n_types and types)
 289                 int m;
 290                 uint32_t *aux;
 291                 aux = (uint32_t*)calloc(n+1, 4);
 292                 m = 0;
 293                 aux[m++] = MINUS_CONST; // zero indel is always a type
 294                 for (i = 0; i < n; ++i) {
 295                         const bam_pileup1_t *p = pl + i;
 296                         if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
 297                                 aux[m++] = MINUS_CONST + p->indel;
 298                 }
 299                 ks_introsort(uint32_t, m, aux);
 300                 n_types = 1;
 301                 for (i = 1; i < m; ++i)
 302                         if (aux[i] != aux[i-1]) ++n_types;
 303                 types = (int*)calloc(n_types, sizeof(int));
 304                 j = 0;
 305                 types[j++] = aux[0] - MINUS_CONST;
 306                 for (i = 1; i < m; ++i) {
 307                         if (aux[i] != aux[i-1])
 308                                 types[j++] = aux[i] - MINUS_CONST;
 309                 }
 310                 free(aux);
 311         }
 312         { // calculate left and right boundary
 313                 bam_segreg_t seg;
 314                 left = 0x7fffffff; right = 0;
 315                 for (i = 0; i < n; ++i) {
 316                         const bam_pileup1_t *p = pl + i;
 317                         if (!(p->b->core.flag&BAM_FUNMAP)) {
 318                                 bam_segreg(pos, &p->b->core, bam1_cigar(p->b), &seg);
 319                                 if (seg.tbeg < left) left = seg.tbeg;
 320                                 if (seg.tend > right) right = seg.tend;
 321                         }
 322                 }
 323         }
 324         { // the core part
 325                 char *ref2, *inscns = 0;
 326                 int k, l, *score, max_ins = types[n_types-1];
 327                 ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
 328                 if (max_ins > 0) { // get the consensus of inserted sequences
 329                         int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
 330                         // count occurrences
 331                         for (i = 0; i < n_types; ++i) {
 332                                 if (types[i] <= 0) continue; // not insertion
 333                                 for (j = 0; j < n; ++j) {
 334                                         const bam_pileup1_t *p = pl + j;
 335                                         if (!(p->b->core.flag&BAM_FUNMAP) && p->indel == types[i]) {
 336                                                 for (k = 1; k <= p->indel; ++k) {
 337                                                         int c = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos + k)];
 338                                                         if (c < 4) ++inscns_aux[i*max_ins*4 + (k-1)*4 + c];
 339                                                 }
 340                                         }
 341                                 }
 342                         }
 343                         // construct the consensus
 344                         inscns = (char*)calloc(n_types * max_ins, sizeof(char));
 345                         for (i = 0; i < n_types; ++i) {
 346                                 for (j = 0; j < types[i]; ++j) {
 347                                         int max = 0, max_k = -1, *ia = inscns_aux + i*max_ins*4 + j*4;
 348                                         for (k = 0; k < 4; ++k) {
 349                                                 if (ia[k] > max) {
 350                                                         max = ia[k];
 351                                                         max_k = k;
 352                                                 }
 353                                         }
 354                                         inscns[i*max_ins + j] = max? 1<<max_k : 15;
 355                                 }
 356                         }
 357                         free(inscns_aux);
 358                 }
 359                 // calculate score
 360                 score = (int*)calloc(n_types * n, sizeof(int));
 361                 for (i = 0; i < n_types; ++i) {
 362                         // write ref2
 363                         for (k = 0, j = left; j <= pos; ++j)
 364                                 ref2[k++] = bam_nt16_table[(int)ref[j]];
 365                         if (types[i] <= 0) j += -types[i];
 366                         else for (l = 0; l < types[i]; ++l)
 367                                          ref2[k++] = inscns[i*max_ins + l];
 368                         for (; j < right && ref[j]; ++j)
 369                                 ref2[k++] = bam_nt16_table[(int)ref[j]];
 370                         // calculate score for each read
 371                         for (j = 0; j < n; ++j) {
 372                                 const bam_pileup1_t *p = pl + j;
 373                                 uint32_t *cigar;
 374                                 bam1_core_t *c = &p->b->core;
 375                                 int s;
 376                                 bam_segreg_t seg;
 377                                 if (c->flag&BAM_FUNMAP) continue;
 378                                 cigar = bam1_cigar(p->b);
 379                                 bam_segreg(pos, c, cigar, &seg);
 380                                 for (s = 0, l = seg.qbeg; c->pos + l < right && l < seg.qend; ++l) {
 381                                         int cq = bam1_seqi(bam1_seq(p->b), l), ct;
 382                                         ct = c->pos + l >= left? ref2[c->pos + l - left] : 15; // "<" should not happen if there is no bug
 383                                         if (cq < 15 && ct < 15)
 384                                                 s += cq == ct? 1 : -mi->mm_penalty;
 385                                 }
 386                                 score[i*n + j] = s;
 387                                 if (types[i] != 0) { // then try the other way to calculate the score
 388                                         for (s = 0, l = seg.qbeg; c->pos + l + types[i] < right && l < seg.qend; ++l) {
 389                                                 int cq = bam1_seqi(bam1_seq(p->b), l), ct;
 390                                                 ct = c->pos + l + types[i] >= left? ref2[c->pos + l + types[i] - left] : 15;
 391                                                 if (cq < 15 && ct < 15)
 392                                                         s += cq == ct? 1 : -mi->mm_penalty;
 393                                         }
 394                                 }
 395                                 if (score[i*n+j] < s) score[i*n+j] = s; // choose the higher of the two scores
 396                                 if (types[i] != 0) score[i*n+j] -= mi->indel_err;
 397                                 //printf("%d, %d, %d, %d\n", i, types[i], j, score[i*n+j]);
 398                         }
 399                 }
 400                 { // get final result
 401                         int *sum, max1, max2, max1_i, max2_i;
 402                         // pick up the best two score
 403                         sum = (int*)calloc(n_types, sizeof(int));
 404                         for (i = 0; i < n_types; ++i)
 405                                 for (j = 0; j < n; ++j)
 406                                         sum[i] += score[i*n+j];
 407                         max1 = max2 = -0x7fffffff; max1_i = max2_i = -1;
 408                         for (i = 0; i < n_types; ++i) {
 409                                 if (sum[i] > max1) {
 410                                         max2 = max1; max2_i = max1_i; max1 = sum[i]; max1_i = i;
 411                                 } else if (sum[i] > max2) {
 412                                         max2 = sum[i]; max2_i = i;
 413                                 }
 414                         }
 415                         free(sum);
 416                         // write ret
 417                         ret = (bam_maqindel_ret_t*)calloc(1, sizeof(bam_maqindel_ret_t));
 418                         ret->indel1 = types[max1_i]; ret->indel2 = types[max2_i];
 419                         ret->s1 = (char*)calloc(abs(ret->indel1) + 2, 1);
 420                         ret->s2 = (char*)calloc(abs(ret->indel2) + 2, 1);
 421                         if (ret->indel1 > 0) {
 422                                 ret->s1[0] = '+';
 423                                 for (k = 0; k < ret->indel1; ++k)
 424                                         ret->s1[k+1] = bam_nt16_rev_table[(int)inscns[max1_i*max_ins + k]];
 425                         } else if (ret->indel1 < 0) {
 426                                 ret->s1[0] = '-';
 427                                 for (k = 0; k < -ret->indel1 && ref[pos + k + 1]; ++k)
 428                                         ret->s1[k+1] = ref[pos + k + 1];
 429                         } else ret->s1[0] = '*';
 430                         if (ret->indel2 > 0) {
 431                                 ret->s2[0] = '+';
 432                                 for (k = 0; k < ret->indel2; ++k)
 433                                         ret->s2[k+1] = bam_nt16_rev_table[(int)inscns[max2_i*max_ins + k]];
 434                         } else if (ret->indel2 < 0) {
 435                                 ret->s2[0] = '-';
 436                                 for (k = 0; k < -ret->indel2 && ref[pos + k + 1]; ++k)
 437                                         ret->s2[k+1] = ref[pos + k + 1];
 438                         } else ret->s2[0] = '*';
 439                         for (j = 0; j < n; ++j) {
 440                                 if (score[max1_i*n+j] < 0 && score[max2_i*n+j] < 0) ++ret->cnt_anti;
 441                                 else {
 442                                         int diff = score[max1_i*n+j] - score[max2_i*n+j];
 443                                         if (diff > mi->ambi_thres) ++ret->cnt1;
 444                                         else if (diff < -mi->ambi_thres) ++ret->cnt2;
 445                                         else ++ret->cnt_ambi;
 446                                 }
 447                         }
 448                 }
 449                 free(score); free(ref2); free(inscns);
 450         }
 451         free(types);
 452         return ret;
 453 }