samtools/bcftools/fet.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <math.h>
   4 #include <stdlib.h>
   5
   6 /* This program is implemented with ideas from this web page:
   7  *
   8  *   http://www.langsrud.com/fisher.htm
   9  */
  10
  11 // log\binom{n}{k}
  12 static double lbinom(int n, int k)
  13 {
  14         if (k == 0 || n == k) return 0;
  15         return lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
  16 }
  17
  18 // n11  n12  | n1_
  19 // n21  n22  | n2_
  20 //-----------+----
  21 // n_1  n_2  | n
  22
  23 // hypergeometric distribution
  24 static double hypergeo(int n11, int n1_, int n_1, int n)
  25 {
  26         return exp(lbinom(n1_, n11) + lbinom(n-n1_, n_1-n11) - lbinom(n, n_1));
  27 }
  28
  29 typedef struct {
  30         int n11, n1_, n_1, n;
  31         double p;
  32 } hgacc_t;
  33
  34 // incremental version of hypergenometric distribution
  35 static double hypergeo_acc(int n11, int n1_, int n_1, int n, hgacc_t *aux)
  36 {
  37         if (n1_ || n_1 || n) {
  38                 aux->n11 = n11; aux->n1_ = n1_; aux->n_1 = n_1; aux->n = n;
  39         } else { // then only n11 changed; the rest fixed
  40                 if (n11%11 && n11 + aux->n - aux->n1_ - aux->n_1) {
  41                         if (n11 == aux->n11 + 1) { // incremental
  42                                 aux->p *= (double)(aux->n1_ - aux->n11) / n11
  43                                         * (aux->n_1 - aux->n11) / (n11 + aux->n - aux->n1_ - aux->n_1);
  44                                 aux->n11 = n11;
  45                                 return aux->p;
  46                         }
  47                         if (n11 == aux->n11 - 1) { // incremental
  48                                 aux->p *= (double)aux->n11 / (aux->n1_ - n11)
  49                                         * (aux->n11 + aux->n - aux->n1_ - aux->n_1) / (aux->n_1 - n11);
  50                                 aux->n11 = n11;
  51                                 return aux->p;
  52                         }
  53                 }
  54                 aux->n11 = n11;
  55         }
  56         aux->p = hypergeo(aux->n11, aux->n1_, aux->n_1, aux->n);
  57         return aux->p;
  58 }
  59
  60 double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two)
  61 {
  62         int i, j, max, min;
  63         double p, q, left, right;
  64         hgacc_t aux;
  65         int n1_, n_1, n;
  66
  67         n1_ = n11 + n12; n_1 = n11 + n21; n = n11 + n12 + n21 + n22; // calculate n1_, n_1 and n
  68         max = (n_1 < n1_) ? n_1 : n1_; // max n11, for right tail
  69         min = n1_ + n_1 - n;
  70         if (min < 0) min = 0; // min n11, for left tail
  71         *two = *_left = *_right = 1.;
  72         if (min == max) return 1.; // no need to do test
  73         q = hypergeo_acc(n11, n1_, n_1, n, &aux); // the probability of the current table
  74         // left tail
  75         p = hypergeo_acc(min, 0, 0, 0, &aux);
  76         for (left = 0., i = min + 1; p < 0.99999999 * q; ++i) // loop until underflow
  77                 left += p, p = hypergeo_acc(i, 0, 0, 0, &aux);
  78         --i;
  79         if (p < 1.00000001 * q) left += p;
  80         else --i;
  81         // right tail
  82         p = hypergeo_acc(max, 0, 0, 0, &aux);
  83         for (right = 0., j = max - 1; p < 0.99999999 * q; --j) // loop until underflow
  84                 right += p, p = hypergeo_acc(j, 0, 0, 0, &aux);
  85         ++j;
  86         if (p < 1.00000001 * q) right += p;
  87         else ++j;
  88         // two-tail
  89         *two = left + right;
  90         if (*two > 1.) *two = 1.;
  91         // adjust left and right
  92         if (abs(i - n11) < abs(j - n11)) right = 1. - left + q;
  93         else left = 1.0 - right + q;
  94         *_left = left; *_right = right;
  95         return q;
  96 }
  97
  98 #ifdef FET_MAIN
  99 #include <stdio.h>
 100
 101 int main(int argc, char *argv[])
 102 {
 103         char id[1024];
 104         int n11, n12, n21, n22;
 105         double left, right, twotail, prob;
 106
 107         while (scanf("%s%d%d%d%d", id, &n11, &n12, &n21, &n22) == 5) {
 108                 prob = kt_fisher_exact(n11, n12, n21, n22, &left, &right, &twotail);
 109                 printf("%s\t%d\t%d\t%d\t%d\t%.6g\t%.6g\t%.6g\t%.6g\n", id, n11, n12, n21, n22,
 110                                 prob, left, right, twotail);
 111         }
 112         return 0;
 113 }
 114 #endif