← Back to branch summary

~ubuntu-branches/ubuntu/quantal/genometools/quantal-backports

« back to all changes in this revision

Viewing changes to src/external/samtools-0.1.18/bcftools/fet.c

  • Committer: Package Import Robot
  • Author(s): Sascha Steinbiss
  • Date: 2012-07-09 14:10:23 UTC
  • Revision ID: package-import@ubuntu.com-20120709141023-juuu4spm6chqsf9o
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/external/samtools-0.1.18/bcftools/fet.c
 
1
#include <math.h>
 
2
#include <stdlib.h>
 
3
 
 
4
/* This program is implemented with ideas from this web page:
 
5
 *
 
6
 *   http://www.langsrud.com/fisher.htm
 
7
 */
 
8
 
 
9
// log\binom{n}{k}
 
10
static double lbinom(int n, int k)
 
11
{
 
12
        if (k == 0 || n == k) return 0;
 
13
        return lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
 
14
}
 
15
 
 
16
// n11  n12  | n1_
 
17
// n21  n22  | n2_
 
18
//-----------+----
 
19
// n_1  n_2  | n
 
20
 
 
21
// hypergeometric distribution
 
22
static double hypergeo(int n11, int n1_, int n_1, int n)
 
23
{
 
24
        return exp(lbinom(n1_, n11) + lbinom(n-n1_, n_1-n11) - lbinom(n, n_1));
 
25
}
 
26
 
 
27
typedef struct {
 
28
        int n11, n1_, n_1, n;
 
29
        double p;
 
30
} hgacc_t;
 
31
 
 
32
// incremental version of hypergenometric distribution
 
33
static double hypergeo_acc(int n11, int n1_, int n_1, int n, hgacc_t *aux)
 
34
{
 
35
        if (n1_ || n_1 || n) {
 
36
                aux->n11 = n11; aux->n1_ = n1_; aux->n_1 = n_1; aux->n = n;
 
37
        } else { // then only n11 changed; the rest fixed
 
38
                if (n11%11 && n11 + aux->n - aux->n1_ - aux->n_1) {
 
39
                        if (n11 == aux->n11 + 1) { // incremental
 
40
                                aux->p *= (double)(aux->n1_ - aux->n11) / n11
 
41
                                        * (aux->n_1 - aux->n11) / (n11 + aux->n - aux->n1_ - aux->n_1);
 
42
                                aux->n11 = n11;
 
43
                                return aux->p;
 
44
                        }
 
45
                        if (n11 == aux->n11 - 1) { // incremental
 
46
                                aux->p *= (double)aux->n11 / (aux->n1_ - n11)
 
47
                                        * (aux->n11 + aux->n - aux->n1_ - aux->n_1) / (aux->n_1 - n11);
 
48
                                aux->n11 = n11;
 
49
                                return aux->p;
 
50
                        }
 
51
                }
 
52
                aux->n11 = n11;
 
53
        }
 
54
        aux->p = hypergeo(aux->n11, aux->n1_, aux->n_1, aux->n);
 
55
        return aux->p;
 
56
}
 
57
 
 
58
double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two)
 
59
{
 
60
        int i, j, max, min;
 
61
        double p, q, left, right;
 
62
        hgacc_t aux;
 
63
        int n1_, n_1, n;
 
64
 
 
65
        n1_ = n11 + n12; n_1 = n11 + n21; n = n11 + n12 + n21 + n22; // calculate n1_, n_1 and n
 
66
        max = (n_1 < n1_) ? n_1 : n1_; // max n11, for right tail
 
67
        min = n1_ + n_1 - n;
 
68
        if (min < 0) min = 0; // min n11, for left tail
 
69
        *two = *_left = *_right = 1.;
 
70
        if (min == max) return 1.; // no need to do test
 
71
        q = hypergeo_acc(n11, n1_, n_1, n, &aux); // the probability of the current table
 
72
        // left tail
 
73
        p = hypergeo_acc(min, 0, 0, 0, &aux);
 
74
        for (left = 0., i = min + 1; p < 0.99999999 * q; ++i) // loop until underflow
 
75
                left += p, p = hypergeo_acc(i, 0, 0, 0, &aux);
 
76
        --i;
 
77
        if (p < 1.00000001 * q) left += p;
 
78
        else --i;
 
79
        // right tail
 
80
        p = hypergeo_acc(max, 0, 0, 0, &aux);
 
81
        for (right = 0., j = max - 1; p < 0.99999999 * q; --j) // loop until underflow
 
82
                right += p, p = hypergeo_acc(j, 0, 0, 0, &aux);
 
83
        ++j;
 
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