gistfile1.pl

use strict;
package HyperGeometric;

use base 'Exporter';
our @EXPORT = qw( hypergeom );

####################################################################
## To calculate p-value using a hypergeometric test for compounds 
## matched to a particular pathway in a given pathway database
##
## For example: 
## There are totally 1,000 (unique) metabolites in the pathway database;
## Among them, Alanine Synthesis Pathway contains 35 (unique) metabolites;
## User uploaded 100 (unique) metabolites names in query;
## Among them, 10 metabolites matched the Alanine Synthesis Pathway; 
## To calculate this probability by random chance: hypergeom(35,1000,100,10);
####################################################################

# my $path_total = shift; # total cmpd in pathway: 20
# my $db_total = shift;   # total cmpd in the pathway DB: 800
# my $query_total = shift; # total cmpd from query: 3
# my $match_path = shift; # total cmpd that matched the pathway: 1
# 
# my $non_path_total = $db_total-$path_total; # 780
# 
# print hypergeom($path_total, $non_path_total, $query_total, $match_path);

#################sub routine#################
sub hypergeom {
   # There are m "bad" and n "good" balls in an urn.
   # Pick N of them. The probability of i or more successful selections:
   # (m!n!N!(m+n-N)!)/(i!(n-i)!(m+i-N)!(N-i)!(m+n)!)
   my ($n, $m, $N, $i) = @_;
   my $loghyp1 = logfact($m)+logfact($n)+logfact($N)+logfact($m+$n-$N);
   my $loghyp2 = logfact($i)+logfact($n-$i)+logfact($m+$i-$N)+logfact($N-$i)+logfact($m+$n);
   return exp($loghyp1 - $loghyp2);
}

sub logfact {
   return gammln(shift(@_) + 1.0);
}

sub gammln {
  my $xx = shift;
  my @cof = (76.18009172947146, -86.50532032941677,
             24.01409824083091, -1.231739572450155,
             0.12086509738661e-2, -0.5395239384953e-5);
  my $y = my $x = $xx;
  my $tmp = $x + 5.5;
  $tmp -= ($x + .5) * log($tmp);
  my $ser = 1.000000000190015;
  for my $j (0..5) {
     $ser += $cof[$j]/++$y;
  }
  -$tmp + log(2.5066282746310005*$ser/$x);
}