samtools/errmod.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <math.h>
   4 #include "errmod.h"
   5 #include "ksort.h"
   6 KSORT_INIT_GENERIC(uint16_t)
   7
   8 typedef struct __errmod_coef_t {
   9         double *fk, *beta, *lhet;
  10 } errmod_coef_t;
  11
  12 typedef struct {
  13         double fsum[16], bsum[16];
  14         uint32_t c[16];
  15 } call_aux_t;
  16
  17 static errmod_coef_t *cal_coef(double depcorr, double eta)
  18 {
  19         int k, n, q;
  20         long double sum, sum1;
  21         double *lC;
  22         errmod_coef_t *ec;
  23
  24         ec = calloc(1, sizeof(errmod_coef_t));
  25         // initialize ->fk
  26         ec->fk = (double*)calloc(256, sizeof(double));
  27         ec->fk[0] = 1.0;
  28         for (n = 1; n != 256; ++n)
  29                 ec->fk[n] = pow(1. - depcorr, n) * (1.0 - eta) + eta;
  30         // initialize ->coef
  31         ec->beta = (double*)calloc(256 * 256 * 64, sizeof(double));
  32         lC = (double*)calloc(256 * 256, sizeof(double));
  33         for (n = 1; n != 256; ++n) {
  34                 double lgn = lgamma(n+1);
  35                 for (k = 1; k <= n; ++k)
  36                         lC[n<<8|k] = lgn - lgamma(k+1) - lgamma(n-k+1);
  37         }
  38         for (q = 1; q != 64; ++q) {
  39                 double e = pow(10.0, -q/10.0);
  40                 double le = log(e);
  41                 double le1 = log(1.0 - e);
  42                 for (n = 1; n <= 255; ++n) {
  43                         double *beta = ec->beta + (q<<16|n<<8);
  44                         sum1 = sum = 0.0;
  45                         for (k = n; k >= 0; --k, sum1 = sum) {
  46                                 sum = sum1 + expl(lC[n<<8|k] + k*le + (n-k)*le1);
  47                                 beta[k] = -10. / M_LN10 * logl(sum1 / sum);
  48                         }
  49                 }
  50         }
  51         // initialize ->lhet
  52         ec->lhet = (double*)calloc(256 * 256, sizeof(double));
  53         for (n = 0; n < 256; ++n)
  54                 for (k = 0; k < 256; ++k)
  55                         ec->lhet[n<<8|k] = lC[n<<8|k] - M_LN2 * n;
  56         free(lC);
  57         return ec;
  58 }
  59
  60 errmod_t *errmod_init(float depcorr)
  61 {
  62         errmod_t *em;
  63         em = (errmod_t*)calloc(1, sizeof(errmod_t));
  64         em->depcorr = depcorr;
  65         em->coef = cal_coef(depcorr, 0.03);
  66         return em;
  67 }
  68
  69 void errmod_destroy(errmod_t *em)
  70 {
  71         if (em == 0) return;
  72         free(em->coef->lhet); free(em->coef->fk); free(em->coef->beta);
  73         free(em->coef); free(em);
  74 }
  75 // qual:6, strand:1, base:4
  76 int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q)
  77 {
  78         call_aux_t aux;
  79         int i, j, k, w[32];
  80
  81         if (m > m) return -1;
  82         memset(q, 0, m * m * sizeof(float));
  83         if (n == 0) return 0;
  84         // calculate aux.esum and aux.fsum
  85         if (n > 255) { // then sample 255 bases
  86                 ks_shuffle(uint16_t, n, bases);
  87                 n = 255;
  88         }
  89         ks_introsort(uint16_t, n, bases);
  90         memset(w, 0, 32 * sizeof(int));
  91         memset(&aux, 0, sizeof(call_aux_t));
  92         for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
  93                 uint16_t b = bases[j];
  94                 int q = b>>5 < 4? 4 : b>>5;
  95                 if (q > 63) q = 63;
  96                 k = b&0x1f;
  97                 aux.fsum[k&0xf] += em->coef->fk[w[k]];
  98                 aux.bsum[k&0xf] += em->coef->fk[w[k]] * em->coef->beta[q<<16|n<<8|aux.c[k&0xf]];
  99                 ++aux.c[k&0xf];
 100                 ++w[k];
 101         }
 102         // generate likelihood
 103         for (j = 0; j != m; ++j) {
 104                 float tmp1, tmp3;
 105                 int tmp2, bar_e;
 106                 // homozygous
 107                 for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != m; ++k) {
 108                         if (k == j) continue;
 109                         tmp1 += aux.bsum[k]; tmp2 += aux.c[k]; tmp3 += aux.fsum[k];
 110                 }
 111                 if (tmp2) {
 112                         bar_e = (int)(tmp1 / tmp3 + 0.499);
 113                         if (bar_e > 63) bar_e = 63;
 114                         q[j*m+j] = tmp1;
 115                 }
 116                 // heterozygous
 117                 for (k = j + 1; k < m; ++k) {
 118                         int cjk = aux.c[j] + aux.c[k];
 119                         for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i < m; ++i) {
 120                                 if (i == j || i == k) continue;
 121                                 tmp1 += aux.bsum[i]; tmp2 += aux.c[i]; tmp3 += aux.fsum[i];
 122                         }
 123                         if (tmp2) {
 124                                 bar_e = (int)(tmp1 / tmp3 + 0.499);
 125                                 if (bar_e > 63) bar_e = 63;
 126                                 q[j*m+k] = q[k*m+j] = -4.343 * em->coef->lhet[cjk<<8|aux.c[k]] + tmp1;
 127                         } else q[j*m+k] = q[k*m+j] = -4.343 * em->coef->lhet[cjk<<8|aux.c[k]]; // all the bases are either j or k
 128                 }
 129                 for (k = 0; k != m; ++k) if (q[j*m+k] < 0.0) q[j*m+k] = 0.0;
 130         }
 131         return 0;
 132 }