c_cpp/lib/htslib/errmod.c

   1 /*  errmod.c -- revised MAQ error model.
   2
   3     Copyright (C) 2010 Broad Institute.
   4     Copyright (C) 2012, 2013, 2016 Genome Research Ltd.
   5
   6     Author: Heng Li <lh3@sanger.ac.uk>
   7
   8 Permission is hereby granted, free of charge, to any person obtaining a copy
   9 of this software and associated documentation files (the "Software"), to deal
  10 in the Software without restriction, including without limitation the rights
  11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  12 copies of the Software, and to permit persons to whom the Software is
  13 furnished to do so, subject to the following conditions:
  14
  15 The above copyright notice and this permission notice shall be included in
  16 all copies or substantial portions of the Software.
  17
  18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  24 DEALINGS IN THE SOFTWARE.  */
  25
  26 #include <config.h>
  27
  28 #include <math.h>
  29 #include "htslib/hts.h"
  30 #include "htslib/ksort.h"
  31 KSORT_INIT_GENERIC(uint16_t)
  32
  33 struct errmod_t {
  34     double depcorr;
  35     /* table of constants generated for given depcorr and eta */
  36     double *fk, *beta, *lhet;
  37 };
  38
  39 typedef struct {
  40     double fsum[16], bsum[16];
  41     uint32_t c[16];
  42 } call_aux_t;
  43
  44 /* \Gamma(n) = (n-1)! */
  45 #define lfact(n) lgamma(n+1)
  46
  47 /* generates a success * trials table of bionomial probability densities (log transformed) */
  48 static double* logbinomial_table( const int n_size )
  49 {
  50     /* prob distribution for binom var is p(k) = {n! \over k! (n-k)! } p^k (1-p)^{n-k} */
  51     /* this calcs p(k) = {log(n!) - log(k!) - log((n-k)!) */
  52     int k, n;
  53     double *logbinom = (double*)calloc(n_size * n_size, sizeof(double));
  54     for (n = 1; n < n_size; ++n) {
  55         double lfn = lfact(n);
  56         for (k = 1; k <= n; ++k)
  57             logbinom[n<<8|k] = lfn - lfact(k) - lfact(n-k);
  58     }
  59     return logbinom;
  60 }
  61
  62 static void cal_coef(errmod_t *em, double depcorr, double eta)
  63 {
  64     int k, n, q;
  65     long double sum, sum1;
  66     double *lC;
  67
  68     // initialize ->fk
  69     em->fk = (double*)calloc(256, sizeof(double));
  70     em->fk[0] = 1.0;
  71     for (n = 1; n < 256; ++n)
  72         em->fk[n] = pow(1. - depcorr, n) * (1.0 - eta) + eta;
  73
  74     // initialize ->beta
  75     em->beta = (double*)calloc(256 * 256 * 64, sizeof(double));
  76
  77     lC = logbinomial_table( 256 );
  78
  79     for (q = 1; q < 64; ++q) {
  80         double e = pow(10.0, -q/10.0);
  81         double le = log(e);
  82         double le1 = log(1.0 - e);
  83         for (n = 1; n <= 255; ++n) {
  84             double *beta = em->beta + (q<<16|n<<8);
  85             sum1 = sum = 0.0;
  86             for (k = n; k >= 0; --k, sum1 = sum) {
  87                 sum = sum1 + expl(lC[n<<8|k] + k*le + (n-k)*le1);
  88                 beta[k] = -10. / M_LN10 * logl(sum1 / sum);
  89             }
  90         }
  91     }
  92
  93     // initialize ->lhet
  94     em->lhet = (double*)calloc(256 * 256, sizeof(double));
  95     for (n = 0; n < 256; ++n)
  96         for (k = 0; k < 256; ++k)
  97             em->lhet[n<<8|k] = lC[n<<8|k] - M_LN2 * n;
  98     free(lC);
  99 }
 100
 101 /**
 102  * Create errmod_t object with obj.depcorr set to depcorr and initialise
 103  */
 104 errmod_t *errmod_init(double depcorr)
 105 {
 106     errmod_t *em;
 107     em = (errmod_t*)calloc(1, sizeof(errmod_t));
 108     em->depcorr = depcorr;
 109     cal_coef(em, depcorr, 0.03);
 110     return em;
 111 }
 112
 113 /**
 114  * Deallocate an errmod_t object
 115  */
 116 void errmod_destroy(errmod_t *em)
 117 {
 118     if (em == 0) return;
 119     free(em->lhet); free(em->fk); free(em->beta);
 120     free(em);
 121 }
 122
 123 //
 124 // em: error model to fit to data
 125 // m: number of alleles across all samples
 126 // n: number of bases observed in sample
 127 // bases[i]: bases observed in pileup [6 bit quality|1 bit strand|4 bit base]
 128 // q[i*m+j]: (Output) phred-scaled likelihood of each genotype (i,j)
 129 int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q)
 130 {
 131     // Aux
 132     // aux.c is total count of each base observed (ignoring strand)
 133     call_aux_t aux;
 134     // Loop variables
 135     int i, j, k;
 136     // The total count of each base observed per strand
 137     int w[32];
 138
 139     memset(q, 0, m * m * sizeof(float)); // initialise q to 0
 140     if (n == 0) return 0;
 141     // This section randomly downsamples to 255 depth so as not to go beyond our precalculated matrix
 142     if (n > 255) { // if we exceed 255 bases observed then shuffle them to sample and only keep the first 255
 143         ks_shuffle(uint16_t, n, bases);
 144         n = 255;
 145     }
 146     ks_introsort(uint16_t, n, bases);
 147     /* zero out w and aux */
 148     memset(w, 0, 32 * sizeof(int));
 149     memset(&aux, 0, sizeof(call_aux_t));
 150
 151     for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
 152         uint16_t b = bases[j];
 153         /* extract quality and cap at 63 */
 154         int qual = b>>5 < 4? 4 : b>>5;
 155         if (qual > 63) qual = 63;
 156         /* extract base ORed with strand */
 157         int basestrand = b&0x1f;
 158         /* extract base */
 159         int base = b&0xf;
 160         aux.fsum[base] += em->fk[w[basestrand]];
 161         aux.bsum[base] += em->fk[w[basestrand]] * em->beta[qual<<16|n<<8|aux.c[base]];
 162         ++aux.c[base];
 163         ++w[basestrand];
 164     }
 165
 166     // generate likelihood
 167     for (j = 0; j < m; ++j) {
 168         float tmp1, tmp3;
 169         int tmp2;
 170         // homozygous
 171         for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k < m; ++k) {
 172             if (k == j) continue;
 173             tmp1 += aux.bsum[k]; tmp2 += aux.c[k]; tmp3 += aux.fsum[k];
 174         }
 175         if (tmp2) {
 176             q[j*m+j] = tmp1;
 177         }
 178         // heterozygous
 179         for (k = j + 1; k < m; ++k) {
 180             int cjk = aux.c[j] + aux.c[k];
 181             for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i < m; ++i) {
 182                 if (i == j || i == k) continue;
 183                 tmp1 += aux.bsum[i]; tmp2 += aux.c[i]; tmp3 += aux.fsum[i];
 184             }
 185             if (tmp2) {
 186                 q[j*m+k] = q[k*m+j] = -4.343 * em->lhet[cjk<<8|aux.c[k]] + tmp1;
 187             } else q[j*m+k] = q[k*m+j] = -4.343 * em->lhet[cjk<<8|aux.c[k]]; // all the bases are either j or k
 188         }
 189         /* clamp to greater than 0 */
 190         for (k = 0; k < m; ++k) if (q[j*m+k] < 0.0) q[j*m+k] = 0.0;
 191     }
 192
 193     return 0;
 194 }