bcftools/fet.c

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