Changeset 1889:2ff9da5aed53 in orangebioinformatics
 Timestamp:
 10/15/13 12:03:11 (6 months ago)
 Branch:
 default
 rebase_source:
 4ad8392b085eabc1a632d544e7017e8668a82706
 Location:
 orangecontrib/bio
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

orangecontrib/bio/obiGO.py
r1873 r1889 921 921 ## print >> sys.stderr, sorted(allAnnotatedGenes) 922 922 ## return 923 923 924 if len(reference) > len(allAnnotatedGenes): 924 925 mappedReferenceGenes = reference.intersection(allAnnotatedGenes) 
orangecontrib/bio/obiProb.py
r1873 r1889 87 87 88 88 def p_value(self, k, N, m, n): 89 subtract = n  k + 1 > k 90 ## result = sum([math.exp(self._logbin(m, i) + self._logbin(N  m, n  i)) for i in (range(k) if subtract else range(k, n+1))]) 91 result = sum([self.__call__(i, N, m, n) for i in (range(k) if subtract else range(k, n+1))]) 92 ## result /= math.exp(self._logbin(N, n)) 93 return max(1.0  result if subtract else result, 0.0) 89 #From wolfram alpha: 90 #1CDF[HypergeometricDistribution[80, 50, 2000], 201] = 1.75695e16 91 #1CDF[HypergeometricDistribution[80, 50, 2000], 401] = 9.92008e52 94 92 93 if min(n,m)  k + 1 <= k: 94 #starting from k gives the shorter list of values 95 return sum(self.__call__(i, N, m, n) for i in range(k, min(n,m)+1)) 96 else: 97 value = 1.0  sum(self.__call__(i, N, m, n) for i in (range(k))) 98 #if the value is small it is probably inexact due to the limited 99 #precision of floats, as for example (1(110e20)) > 0 100 #if so, compute the result without substraction 101 if value < 1e3: #arbitary threshold 102 #print "INEXACT", value, sum(self.__call__(i, N, m, n) for i in range(k, min(n,m)+1)) 103 return sum(self.__call__(i, N, m, n) for i in range(k, min(n,m)+1)) 104 else: 105 return value 95 106 96 107 ## to speedup FDR, calculate ahead sum([1/i for i in range(1, m+1)]), for m in [1,100000]. For higher values of m use an approximation, with error less or equal to 4.99999157277e006. (sum([1/i for i in range(1, m+1)]) ~ log(m) + 0.5772..., 0.5572 is an EulerMascheroni constant)
Note: See TracChangeset
for help on using the changeset viewer.