address issue #41 to also test for terms missing in 'study' - NOTE th…

…is change

increases the number of statistical tests performed, so p-values after multiple
testing correction will be numerically different from results before

goatools/go_enrichment.py

@@ -87,22 +87,20 @@ def __init__(self, pop, assoc, obo_dag, alpha=.05, study=None,
         self.results = []
 
         obo_dag.update_association(assoc)
+        self.term_study = count_terms(study, assoc, obo_dag)
         self.term_pop = count_terms(pop, assoc, obo_dag)
 
         if study:
             self.run_study(study)
 
     def run_study(self, study):
         results = self.results
-
-        term_study = count_terms(study, self.assoc, self.obo_dag)
-
         pop_n, study_n = len(self.pop), len(study)
+        allterms = set(self.term_study.keys() + self.term_pop.keys())
 
-        # Init study_count and pop_count to handle empty sets
-        study_count = pop_count = 0
-        for term, study_count in list(term_study.items()):
-            pop_count = self.term_pop[term]
+        for term in allterms:
+            study_count = self.term_study.get(term, 0)
+            pop_count = self.term_pop.get(term, 0)
             p = fisher.pvalue_population(study_count, study_n,
                                          pop_count, pop_n)
 
@@ -159,7 +157,11 @@ def update_results(self, method, corrected_pvals):
             rec.__setattr__("p_"+method, val)
 
     def print_summary(self, min_ratio=None, indent=False, pval=0.05):
-        # Header contains parameters
+        from .version import __version__ as version
+        from datetime import date
+
+        # Header contains provenance and parameters
+        print("# Generated by GOATOOLS v{0} ({1})".format(version, date.today()))
         print("# min_ratio={0} pval={1}".format(min_ratio, pval))
 
         # field names for output
@@ -170,7 +172,7 @@ def print_summary(self, min_ratio=None, indent=False, pval=0.05):
             # (over_under, is_ratio_different)
             rec.update_remaining_fields(min_ratio=min_ratio)
 
-            if pval is not None and rec.p_bonferroni > pval:
+            if pval is not None and rec.p_uncorrected >= pval:
                 continue
 
             if rec.is_ratio_different:

goatools/multiple_testing.py

@@ -21,6 +21,8 @@ def __init__(self, pvals, a=.05):
         self.a = a                  # type-1 error cutoff for each test
 
         self.set_correction()
+        # Reset all pvals > 1 to 1
+        self.corrected_pvals[self.corrected_pvals > 1] = 1
 
     def set_correction(self):
         # the purpose of multiple correction is to lower the alpha
@@ -36,9 +38,6 @@ class Bonferroni(AbstractCorrection):
     def set_correction(self):
         self.corrected_pvals *= self.n
 
-        # reset all pvals > 1 to 1
-        self.corrected_pvals[self.corrected_pvals > 1] = 1
-
 
 class Sidak(AbstractCorrection):
     """http://en.wikipedia.org/wiki/Bonferroni_correction

goatools/version.py

@@ -1 +1 @@
-__version__ = "0.5.5"
+__version__ = "0.5.6"

run.sh

@@ -28,7 +28,7 @@ do
 case $REPLY in
 
 1)
-python scripts/find_enrichment.py --alpha=0.05 --indent data/study data/population data/association
+python scripts/find_enrichment.py --pval=0.05 --indent data/study data/population data/association
 ;;
 
 2)

scripts/find_enrichment.py

@@ -84,9 +84,8 @@ def check_input_files(ns, p):
                  help='data/study data/population data/association')
     p.add_argument('--alpha', default=0.05, type=float,
                  help="Test-wise alpha for multiple testing ")
-    p.add_argument('--pval', default=None, type=float,
-                 help="Family-wise alpha (whole experiment), only print out "
-                 "Bonferroni p-value is less than this value. ")
+    p.add_argument('--pval', default=.05, type=float,
+                 help="Only print out when uncorrected p-value < this value.")
     p.add_argument('--compare', dest='compare', default=False,
                  action='store_true',
                  help="the population file as a comparison group. if this "