1 |
#!/usr/bin/perl |
2 |
use strict; |
3 |
use Getopt::Std; |
4 |
use FindBin;use lib $FindBin::Bin; |
5 |
use Fisher qw(fishers_exact log_factorial); |
6 |
|
7 |
|
8 |
my $usage = q/Usage: |
9 |
fisher_exact_test.pl [-T] a1 b1 a2 b2 |
10 |
Computes Fisher's exact test p-values for 2x2 matrix: |
11 |
|
12 |
a1 b1 |
13 |
|
14 |
a2 b2 |
15 |
|
16 |
-T option will make the input values to be taken in terms of totals: |
17 |
|
18 |
t1 t2 n1 n2 |
19 |
|
20 |
..where n1 is a subset of t1 and n2 a subset of t2, i.e. |
21 |
t1=a1+b1 and t2=a2+b2 |
22 |
|
23 |
/; |
24 |
umask 0002; |
25 |
getopts('To:') || die($usage."\n"); |
26 |
die($usage." Error: exactly 4 values are needed (totals first)!\n") unless @ARGV==4; |
27 |
my $totals=$Getopt::Std::opt_T; |
28 |
my @fncache=(0,0); |
29 |
|
30 |
# -- needed by Audic-Calverie: calc_audic() |
31 |
my $MAXIT = 500; |
32 |
my $EPS = 3.0E-30; |
33 |
my $FPMIN = 1.0E-30; |
34 |
my @COF = ( 76.18009172947146, -86.50532032941677, |
35 |
24.01409824083091, -1.231739572450155, |
36 |
0.1208650973866179E-2,-0.5395239384953E-5 ); |
37 |
|
38 |
|
39 |
|
40 |
{ |
41 |
my ($a, $b , $c, $d)=@ARGV; |
42 |
if ($totals) { |
43 |
my ($t1, $t2, $n1,$n2)=@ARGV; |
44 |
die("Error: totals should be larger than parts!\n") |
45 |
unless $n1<$t1 && $n2<$t2; |
46 |
($a, $b , $c, $d)=($n1, $t1-$n1, $n2, $t2-$n2); |
47 |
print STDERR " a b c d = $a $b $c $d\n"; |
48 |
#my $ap=audic($t1, $t2, $n1, $n2); |
49 |
my ($ap, $apdir)=calc_audic($n1, $n2, $t1, $t2, 1); |
50 |
print ">>>> Audic & Claverie : $ap ($apdir)\n"; |
51 |
} |
52 |
my $FisherProb = ProbCTable($a, $b , $c, $d); |
53 |
print "-----> Fisher Probability: $FisherProb\n"; |
54 |
my $fet1=fishers_exact($a,$b,$c,$d); |
55 |
my $fet2=fishers_exact($a,$b,$c,$d, 1); |
56 |
print "Fisher.pm : one-tailed: $fet1\n"; |
57 |
print "Fisher.pm : two-tailed: $fet2\n"; |
58 |
exit; |
59 |
#--- rest is optional: |
60 |
my $a_b = $a + $b; |
61 |
my $c_d = $c + $d; |
62 |
my $a_c = $a + $c; |
63 |
my $b_d = $b + $d; |
64 |
my $total = $a + $b + $c + $d; |
65 |
my $p1 = $a/$a_b; |
66 |
my $p2 = $c/$c_d; |
67 |
my $deltap = abs($p1 - $p2); |
68 |
my $p = $a_c / $total; |
69 |
my $sd = sqrt($p * (1-$p) * (1/$a_b + 1/$c_d)); |
70 |
my $z = abs($p1 - $p2) / $sd; |
71 |
my $zp = abs(ltqnorm($FisherProb)); |
72 |
my $z95 = abs(ltqnorm(0.95)); |
73 |
|
74 |
my $z_ratio = $z / $zp; |
75 |
my $deltap95 = $deltap * ($z95 / $zp); |
76 |
my $ratio95 = $z95 / $zp; |
77 |
print "Deltap:$deltap\tProjectedDeltap95:$deltap95\tRatio:$ratio95\n\n"; |
78 |
|
79 |
print "$a\t$b\t$a_b\n$c\t$d\t$c_d\n". |
80 |
"$a_c\t$b_d\t$total\nFisherProb: $FisherProb <----\n" . |
81 |
"z:$z\nzp:$zp\nratio:$z_ratio\n$p1\t$p2\t$p\n\n"; |
82 |
|
83 |
my $p_this_table = ProbOneTable($a , $b , $c , $d); |
84 |
print "This one table p:$p_this_table\n"; |
85 |
|
86 |
my $mid_p = $FisherProb - $p_this_table/2; |
87 |
print "Mid-p:$mid_p\n\n"; |
88 |
} |
89 |
|
90 |
#================ SUBROUTINES ============ |
91 |
|
92 |
sub LnFactorial { |
93 |
my $n = shift; |
94 |
return $fncache[$n] if defined $fncache[$n]; |
95 |
return undef if $n < 0; |
96 |
for (my $i = scalar(@fncache); $i <= $n; $i++) { |
97 |
$fncache[$i] = $fncache[$i-1] + log($i); |
98 |
} |
99 |
|
100 |
return $fncache[$n]; |
101 |
} |
102 |
sub factorial { |
103 |
my $r=1; |
104 |
$r *= $_ foreach 2..$_[0]; |
105 |
return $r; |
106 |
} |
107 |
|
108 |
#compute Audic & Claverie probability |
109 |
sub audic { |
110 |
my ($n1, $n2, $x, $y)=@_; |
111 |
my $v=$y*log($n2/$n1)+log_factorial($x+$y)-log_factorial($x)-log_factorial($y)- |
112 |
($x+$y+1)*log(1+$n2/$n1); |
113 |
return exp($v); |
114 |
#my $nratio=$n2/$n1; |
115 |
#my $v=($nratio ** $y)*factorial($x+$y)/(factorial($x)*factorial($y)*((1+$nratio)**($x+$y+1))); |
116 |
|
117 |
} |
118 |
|
119 |
sub calc_audic { |
120 |
=pod |
121 |
|
122 |
=head2 calc_audic $x, $y, $Nx, $Ny, <$signedValue> |
123 |
|
124 |
Determines the statistical significance of the difference |
125 |
in tag count (expression) between two libraries. This |
126 |
function uses the method described by Audic and |
127 |
Claverie (1997). This method can be called on an |
128 |
instantiated object, as well as statically. |
129 |
|
130 |
B<Arguments> |
131 |
|
132 |
I<$x,$y> |
133 |
|
134 |
The number of tags in the x- and y-axis |
135 |
libraries, respectively. |
136 |
|
137 |
I<$Nx,$Ny> |
138 |
|
139 |
The total number of tags in the x- and y-axis |
140 |
libraries, respectively. |
141 |
|
142 |
I<$signedValue> (optional) |
143 |
|
144 |
A boolean value (>=1 is FALSE). If this flag is |
145 |
set to TRUE, downregulated comparisons will return |
146 |
both a p-value and either +1, -1, or 0 to indicate |
147 |
up/down/same-regulation (i.e. -1 if the expression |
148 |
ratio of tags in the x-axis library(s) is greater |
149 |
than that of the y-axis library(s)). This flag |
150 |
is FALSE by default. |
151 |
|
152 |
B<Returns> |
153 |
|
154 |
The p-value for the observation. A lower number is |
155 |
more significant. Typically, observations with |
156 |
p <= 0.05 are considered statistically significant. |
157 |
|
158 |
If $signedValue is set to TRUE, the function also |
159 |
returns a 0, -1 or +1 to indicate same/down/up-regulation. |
160 |
|
161 |
B<Usage> |
162 |
|
163 |
# the function is static, so it can be accessed directly |
164 |
my $p = calc_audic( 3, 10, 50000, 60000 ); |
165 |
|
166 |
# or: |
167 |
my ( $p, $sign ) = calc_audic( 3, 10, 50000, 60000, 1 ); |
168 |
if( $p <= 0.05 ) { |
169 |
if( $sign == +1 ) { print "Significantly upregulated.\n"; } |
170 |
if( $sign == -1 ) { print "Significantly downregulated.\n"; } |
171 |
if( $sign == 0 ) { die( "Same expression should never be significant!" ); } |
172 |
} |
173 |
|
174 |
=cut |
175 |
my $x = shift; |
176 |
if( !defined( $x ) ) { die( " calc_audic no arguments provided." ); } |
177 |
my $y = shift; |
178 |
my $M = shift; # cf n1 |
179 |
my $N = shift; # cf n2 |
180 |
my $bSign = shift || 0; |
181 |
|
182 |
my $p = $M / ( $M+$N ); |
183 |
|
184 |
my $thisproba = __betai( $x+1, $y+1, $p ); |
185 |
my $thisproba2 = __betai( $y+1, $x+1, 1.0-$p ); |
186 |
|
187 |
if( $bSign >= 1 ) { |
188 |
my $ratio1 = $x / $M; |
189 |
my $ratio2 = $y / $N; |
190 |
my $sign = 0; |
191 |
if( $ratio1 > $ratio2 ) { $sign = -1; } |
192 |
if( $ratio1 < $ratio2 ) { $sign = +1; } |
193 |
#return ( ( $thisproba < $thisproba2 ? ( 2*$thisproba ) : ( 2*$thisproba2 ) ), $sign ); |
194 |
return ( ($thisproba < $thisproba2) ? $thisproba : $thisproba2, $sign ); |
195 |
} |
196 |
|
197 |
# return ( $thisproba < $thisproba2 ? ( 2*$thisproba ) : ( 2*$thisproba2 ) ); |
198 |
return ( $thisproba < $thisproba2) ? $thisproba : $thisproba2 ; |
199 |
|
200 |
} |
201 |
|
202 |
################################### |
203 |
# Audic and Claverie C->Perl Port # |
204 |
################################### |
205 |
|
206 |
sub __gammln { |
207 |
|
208 |
my $xx = shift; |
209 |
|
210 |
my $x = $xx; |
211 |
my $y = $xx; |
212 |
|
213 |
my $tmp = $x + 5.5; |
214 |
$tmp -= ( $x + 0.5 ) * log( $tmp ); |
215 |
my $ser = 1.000000000190015; |
216 |
for( my $j = 0; $j <= 5; $j++ ) { $ser += $COF[$j] / ++$y; } |
217 |
|
218 |
return -$tmp + log( 2.5066282746310005 * $ser / $x ); |
219 |
|
220 |
} |
221 |
|
222 |
sub __betai { |
223 |
my $a = shift; |
224 |
my $b = shift; |
225 |
my $x = shift; |
226 |
|
227 |
if( $x < 0.0 || $x > 1.0 ) { |
228 |
die( "Bad x in routine betai." ); |
229 |
} |
230 |
|
231 |
my $bt; |
232 |
|
233 |
if( $x == 0.0 || $x == 1.0 ) { |
234 |
$bt = 0.0; |
235 |
} else { |
236 |
$bt = exp( __gammln( $a+$b ) - __gammln( $a ) - __gammln( $b ) + $a*log( $x ) + $b*log( 1.0-$x ) ); |
237 |
} |
238 |
|
239 |
if( $x < ( $a+1.0 )/( $a+$b+2.0 ) ) { |
240 |
return $bt * __betacf( $a, $b, $x ) / $a; |
241 |
} |
242 |
|
243 |
return 1.0 - $bt * __betacf( $b, $a, 1.0-$x ) / $b; |
244 |
|
245 |
} |
246 |
|
247 |
sub __fabs { |
248 |
my $x = shift; |
249 |
return ( $x < 0 ? -$x : $x ); |
250 |
} |
251 |
|
252 |
sub __betacf { |
253 |
|
254 |
my $a = shift; |
255 |
my $b = shift; |
256 |
my $x = shift; |
257 |
|
258 |
my $qab = $a + $b; |
259 |
my $qap = $a + 1.0; |
260 |
my $qam = $a - 1.0; |
261 |
my $c = 1.0; |
262 |
my $d = 1.0 - $qab * $x / $qap; |
263 |
|
264 |
if( __fabs( $d ) < $FPMIN ) { $d = $FPMIN; } |
265 |
|
266 |
$d = 1.0 / $d; # inverse d |
267 |
my $h = $d; |
268 |
my $m; |
269 |
for( $m = 1; $m <= $MAXIT; $m++ ) { |
270 |
my $m2 = 2 * $m; |
271 |
my $aa = $m * ( $b-$m ) * $x / ( ( $qam + $m2 ) * ( $a + $m2 ) ); |
272 |
|
273 |
$d = 1.0 + $aa*$d; |
274 |
if( __fabs( $d ) < $FPMIN ) { $d = $FPMIN; } |
275 |
|
276 |
$c = 1.0 + $aa/$c; |
277 |
if( __fabs( $c ) < $FPMIN ) { $c = $FPMIN; } |
278 |
|
279 |
$d = 1.0 / $d; # inverse d |
280 |
|
281 |
$h *= $d*$c; |
282 |
|
283 |
$aa = -($a+$m)*($qab+$m)*$x/(($a+$m2)*($qap+$m2)); |
284 |
|
285 |
$d = 1.0 + $aa * $d; |
286 |
if( __fabs( $d ) < $FPMIN ) { $d = $FPMIN; } |
287 |
|
288 |
$c = 1.0 + $aa / $c; |
289 |
if( __fabs( $c ) < $FPMIN ) { $c = $FPMIN; } |
290 |
|
291 |
$d = 1.0 / $d; # inverse d; |
292 |
|
293 |
my $del = $d*$c; |
294 |
$h *= $del; |
295 |
if( __fabs( $del-1.0 ) < $EPS ) { last; } |
296 |
} |
297 |
|
298 |
if( $m > $MAXIT ) { |
299 |
# die( "a or b too big, or MAXIT too small in __betacf" ); |
300 |
print STDERR "a($a) or b($b) too big, or MAXIT($MAXIT) too small in __betacf!\n"; |
301 |
} |
302 |
return $h; |
303 |
} |
304 |
|
305 |
|
306 |
|
307 |
# Compute the probability of getting this exact table |
308 |
# using the hypergeometric distribution |
309 |
sub ProbOneTable { |
310 |
my ($a , $b , $c, $d) = @_; |
311 |
my $n = $a + $b + $c + $d; |
312 |
my $LnNumerator = LnFactorial($a+$b)+ |
313 |
LnFactorial($c+$d)+ |
314 |
LnFactorial($a+$c)+ |
315 |
LnFactorial($b+$d); |
316 |
|
317 |
my $LnDenominator = LnFactorial($a) + |
318 |
LnFactorial($b) + |
319 |
LnFactorial($c) + |
320 |
LnFactorial($d) + |
321 |
LnFactorial($n); |
322 |
|
323 |
my $LnP = $LnNumerator - $LnDenominator; |
324 |
return exp($LnP); |
325 |
} |
326 |
|
327 |
# Compute the cumulative probability by adding up individual |
328 |
# probabilities |
329 |
sub ProbCTable { |
330 |
my ($a, $b, $c, $d) = @_; |
331 |
|
332 |
my $min; |
333 |
|
334 |
my $n = $a + $b + $c + $d; |
335 |
|
336 |
my $p = 0; |
337 |
$p += ProbOneTable($a, $b, $c, $d); |
338 |
if( ($a * $d) >= ($b * $c) ) { |
339 |
$min = ($c < $b) ? $c : $b; |
340 |
for(my $i = 0; $i < $min; $i++) { |
341 |
$p += ProbOneTable(++$a, --$b, --$c, ++$d); |
342 |
} |
343 |
} |
344 |
|
345 |
if ( ($a * $d) < ($b * $c) ) { |
346 |
$min = ($a < $d) ? $a : $d; |
347 |
for(my $i = 0; $i < $min; $i++) { |
348 |
$p += ProbOneTable(--$a, ++$b, ++$c, --$d); |
349 |
} |
350 |
} |
351 |
return $p; |
352 |
} |
353 |
|
354 |
# Lower tail quantile for standard normal distribution function. |
355 |
# |
356 |
# This function returns an approximation of the inverse cumulative |
357 |
# standard normal distribution function. I.e., given P, it returns |
358 |
# an approximation to the X satisfying P = Pr{Z <= X} where Z is a |
359 |
# random variable from the standard normal distribution. |
360 |
# |
361 |
# The algorithm uses a minimax approximation by rational functions |
362 |
# and the result has a relative error whose absolute value is less |
363 |
# than 1.15e-9. |
364 |
# |
365 |
# Author: Peter J. Acklam |
366 |
# Time-stamp: 2000-07-19 18:26:14 |
367 |
# E-mail: pjacklam@online.no |
368 |
# WWW URL: http://home.online.no/~pjacklam |
369 |
sub ltqnorm { |
370 |
my $p = shift; |
371 |
die "input argument must be in (0,1)\n" unless 0 < $p && $p < 1; |
372 |
|
373 |
# Coefficients in rational approximations. |
374 |
my @a = (-3.969683028665376e+01, 2.209460984245205e+02, |
375 |
-2.759285104469687e+02, 1.383577518672690e+02, |
376 |
-3.066479806614716e+01, 2.506628277459239e+00); |
377 |
my @b = (-5.447609879822406e+01, 1.615858368580409e+02, |
378 |
-1.556989798598866e+02, 6.680131188771972e+01, |
379 |
-1.328068155288572e+01 ); |
380 |
my @c = (-7.784894002430293e-03, -3.223964580411365e-01, |
381 |
-2.400758277161838e+00, -2.549732539343734e+00, |
382 |
4.374664141464968e+00, 2.938163982698783e+00); |
383 |
my @d = ( 7.784695709041462e-03, 3.224671290700398e-01, |
384 |
2.445134137142996e+00, 3.754408661907416e+00); |
385 |
|
386 |
# Define break-points. |
387 |
my $plow = 0.02425; |
388 |
my $phigh = 1 - $plow; |
389 |
|
390 |
# Rational approximation for lower region: |
391 |
if ( $p < $plow ) { |
392 |
my $q = sqrt(-2*log($p)); |
393 |
return ((((($c[0]*$q+$c[1])*$q+$c[2])*$q+$c[3])*$q+$c[4])*$q+$c[5]) / |
394 |
(((($d[0]*$q+$d[1])*$q+$d[2])*$q+$d[3])*$q+1); |
395 |
} |
396 |
|
397 |
# Rational approximation for upper region: |
398 |
if ( $phigh < $p ) { |
399 |
my $q = sqrt(-2*log(1-$p)); |
400 |
return -((((($c[0]*$q+$c[1])*$q+$c[2])*$q+$c[3])*$q+$c[4])*$q+$c[5]) / |
401 |
(((($d[0]*$q+$d[1])*$q+$d[2])*$q+$d[3])*$q+1); |
402 |
} |
403 |
|
404 |
# Rational approximation for central region: |
405 |
my $q = $p - 0.5; |
406 |
my $r = $q*$q; |
407 |
return ((((($a[0]*$r+$a[1])*$r+$a[2])*$r+$a[3])*$r+$a[4])*$r+$a[5])*$q / |
408 |
((((($b[0]*$r+$b[1])*$r+$b[2])*$r+$b[3])*$r+$b[4])*$r+1); |
409 |
} |