Repository: monarch-initiative/gpsea (language: python)

Issue:
Streamline predicate/analysis workflow
Currently, all `PolyPredicate` implementations return these values as available categories:
```
HETEROZYGOUS, HOMOZYGOUS, NO_VARIANT
```

However, the categories do not make intuitive sense in context of several predicates. 

For instance `VariantEffectPredicate` should test if the patient has a variant with the `query: VariantEffect` on the given transcript or not. So, the predicate should in fact be a simple predicate (not poly) and return `True`/`False` for a patient that has or has not such variant.

`ExonPredicate` is another example. The predicate should return `True`/`False` based on if the variant affects the exon in question or not.

Having more accurate predicate categories, the `CohortAnalyzer` should be able to treat the predicates the same, which simplifies the code. It should be easy to  apply an additional logic. For instance, filter out the patients that do not meet the expected mode of inheritance, count the groups, do statistical testing, apply multiple-testing correction, etc..



Fix change length bug
Fix #93 

The PR fixes a bug in change length calculation, removes trimming of the end coordinate before sending out to VEP for annotation, and prettifies the `GenomicInterpretation` JSON snippets.

Moreover, a documentation section regarding tests is added. From now on, we have 2 kinds of test, online and offline. The _online_ tests must be marked by `@pytest.mark.online` to not run them in standard CI due to spurious failures. In the long run, we should remove the online tests altogether and use mocks instead. The _offline_ tests are just the regular tests.


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 8591fc91..11799af5 100644
--- a/README.md
+++ b/README.md
@@ -3,35 +3,41 @@
 ![PyPi downloads](https://img.shields.io/pypi/dm/genophenocorr.svg?label=Pypi%20downloads)
 ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/genophenocorr)
 
-Genophenocorr is a Python library for genotype-phenotype association analysis. 
+Genophenocorr is a Python library for genotype-phenotype association analysis.
 
 An example of simple genotype-phenotype association analysis
+
 ```python
 # Load HPO
 import hpotk
+
 hpo = hpotk.load_minimal_ontology('http://purl.obolibrary.org/obo/hp.json')
 
-# Load a cohort of phenopackets 
+# Load a cohort of phenopackets
 from genophenocorr.data import get_toy_cohort
+
 cohort = get_toy_cohort()
 
-# Analyze genotype-phenotype associations 
-from genophenocorr.analysis import CohortAnalysis
-from genophenocorr.constants import VariantEffect
+# Analyze genotype-phenotype associations
+from genophenocorr.analysis import configure_cohort_analysis
+from genophenocorr.analysis.predicate import BooleanPredicate
+from genophenocorr.model import VariantEffect
+
+cohort_analysis = configure_cohort_analysis(cohort, hpo)
+frameshift = cohort_analysis.compare_by_variant_effect(VariantEffect.FRAMESHIFT_VARIANT, tx_id='NM_1234.5')
 
-cohort_analysis = CohortAnalysis(cohort, 'NM_1234.5', hpo)
-frameshift = cohort_analysis.compare_by_variant_type(VariantEffect.FRAMESHIFT_VARIANT)
-print(frameshift)
+frameshift.summarize(hpo, phenotype_category=BooleanPredicate.YES)
 ```
 
-prints a table with genotype-phenotype correlations:
+provides a pandas data frame with genotype-phenotype correlations:
 
 ```text
-                            With frameshift_variant         Without frameshift_variant
-                                              Count Percent                      Count Percent  p-value
-HP:0001166 (Arachnodactyly)                       4  30.77%                         10  76.92%  0.04718
-HP:0001250 (Seizure)                             11  84.62%                          9  69.23%  0.64472
-HP:0001257 (Spasticity)                           8  61.54%                          9  69.23%  1.00000
+FRAMESHIFT_VARIANT on NM_1234.5        No             Yes
+                                    Count   Percent Count Percent   p value Corrected p value
+    Arachnodactyly [HP:0001166]         1      3.84    13    50.0   0.00078          0.020299
+    Seizure [HP:0001250]               11     84.62     9    69.2   0.64472          0.913432
+    Spasticity [HP:0001257]             8     61.54     9    69.2   1.00000          1.000000
+    ...                               ...       ...    ...    ...       ...               ...
 ```
 
 ## Documentation
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index 2c32821d..d455c514 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -37,10 +37,10 @@ with the package:
 
   See :ref:`input-data` section to learn about preparing your data for the analysis.
 
-We can then view the data using the list commands. 
+We can then view the data using the list commands.
 
 .. doctest:: tutorial
-  
+
   >>> sorted(cohort.list_all_patients())
   ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z']
   >>> sorted(cohort.list_all_phenotypes())
@@ -54,34 +54,30 @@ We can then view the data using the list commands.
   [('FRAMESHIFT_VARIANT', 1), ('MISSENSE_VARIANT', 1)]
 
 Using the counts, we can choose and run what analyses we want.
-For instance, we can partition the patients into two groups based on presence/absence of a *frameshift* variant:
+For instance, we can partition the patients into two groups based on presence/absence of a *missense* variant:
 
 .. doctest:: tutorial
 
-  >>> from genophenocorr.analysis import CohortAnalysis
+  >>> from genophenocorr.analysis import configure_cohort_analysis
+  >>> from genophenocorr.analysis.predicate import BooleanPredicate  # TODO - explain the predicate or update the API
   >>> from genophenocorr.model import VariantEffect
-  >>> cohort_analysis = CohortAnalysis(cohort, 'NM_1234.5', hpo, include_unmeasured=False)
-  >>> frameshift = cohort_analysis.compare_by_variant_type(VariantEffect.FRAMESHIFT_VARIANT)
-  >>> frameshift # doctest: +NORMALIZE_WHITESPACE
-                              With frameshift_variant         Without frameshift_variant
-                                                Count Percent                      Count Percent  p-value Corrected p-values
-  HP:0001166 (Arachnodactyly)                       4  30.77%                         10  76.92%  0.04718            0.14154
-  HP:0001250 (Seizure)                             11  84.62%                          9  69.23%  0.64472            1.00000
-  HP:0001257 (Spasticity)                           8  61.54%                          9  69.23%  1.00000            1.00000
 
+  >>> cohort_analysis = configure_cohort_analysis(cohort, hpo)
+  >>> missense = cohort_analysis.compare_by_variant_effect(VariantEffect.MISSENSE_VARIANT, tx_id='NM_1234.5')
+  >>> summary_df = missense.summarize(hpo, BooleanPredicate.YES)
+  >>> summary_df.head(1) # doctest: +NORMALIZE_WHITESPACE
+    MISSENSE_VARIANT on NM_1234.5    No             Yes
+                                    Count   Percent Count Percent   p value Corrected p value
+    Arachnodactyly [HP:0001166]         1  3.846154    13    50.0  0.000781          0.020299
 
-Or perform similar partitioning based on presence/absence of a *missense* variant:
+..
 
-.. doctest:: tutorial
+  We're showing just 1 row above. This is due to 2-.. rows all having corrected p value of `1.000` resulting
+  in unstable sort order. We can show more rows with a better cohort, as soon as we have it!
 
-  >>> missense = cohort_analysis.compare_by_variant_type(VariantEffect.MISSENSE_VARIANT)
-  >>> missense # doctest: +NORMALIZE_WHITESPACE
-                              With missense_variant         Without missense_variant
-                                              Count Percent                    Count Percent   p-value Corrected p-values
-  HP:0001166 (Arachnodactyly)                    13  81.25%                        1  10.00%  0.000781           0.002342
-  HP:0001257 (Spasticity)                        11  68.75%                        6  60.00%  0.692449           1.000000
-  HP:0001250 (Seizure)                           12  75.00%                        8  80.00%  1.000000           1.000000
+..
 
+  We can show analysis for `VariantEffect.FRAMESHIFT_VARIANT` as well..
 
-The tables present the HPO terms that annotate the cohort members and report their counts and p values
-for each genotype group. The rows are sorted by the p value in ascending order.
+The table presents the HPO terms that annotate the cohort members and report their counts and p values
+for each genotype group. The rows are sorted by the corrected p value in ascending order.
diff --git a/src/genophenocorr/analysis/__init__.py b/src/genophenocorr/analysis/__init__.py
index e29bbd50..7fa96f4d 100644
--- a/src/genophenocorr/analysis/__init__.py
+++ b/src/genophenocorr/analysis/__init__.py
@@ -1,7 +1,10 @@
 from . import predicate
 
-from ._analyzers import CohortAnalysis
+from ._api import CohortAnalysis, GenotypePhenotypeAnalysisResult
+from ._config import CohortAnalysisConfiguration, CohortAnalysisConfigurationBuilder, configure_cohort_analysis
 
 __all__ = [
-    'CohortAnalysis'
+    'configure_cohort_analysis',
+    'CohortAnalysis', 'GenotypePhenotypeAnalysisResult',
+    'CohortAnalysisConfiguration', 'CohortAnalysisConfigurationBuilder'
 ]
diff --git a/src/genophenocorr/analysis/_analyzers.py b/src/genophenocorr/analysis/_analyzers.py
index 83fc0310..82e6ebcd 100644
--- a/src/genophenocorr/analysis/_analyzers.py
+++ b/src/genophenocorr/analysis/_analyzers.py
@@ -1,33 +1,48 @@
-import abc
 import logging
-import math
 import typing
-from decimal import Decimal
+import warnings
 
-import numpy as np
+import hpotk
 
-from scipy import stats
 from statsmodels.stats import multitest
 from pandas import DataFrame, MultiIndex
-from collections import Counter, namedtuple
+from collections import Counter
 
-from genophenocorr.model import Cohort, FeatureType, VariantEffect
+from genophenocorr.model import Cohort, FeatureType, VariantEffect, Patient
 
-from .predicate import VariantEffectPredicate, HPOPresentPredicate, VariantPredicate, ExonPredicate, ProtFeatureTypePredicate, ProtFeaturePredicate
-from .predicate import HOMOZYGOUS, HETEROZYGOUS, NO_VARIANT
+from .predicate import PolyPredicate
+from .predicate.genotype import VariantEffectPredicate, VariantPredicate, ExonPredicate, ProtFeatureTypePredicate, ProtFeaturePredicate
+from .predicate.genotype import HOMOZYGOUS, HETEROZYGOUS, NO_VARIANT
 
+from ._api import CohortAnalysis
+from ._stats import run_recessive_fisher_exact, run_fisher_exact
+
+
+class BaseCohortAnalysis(CohortAnalysis):
+    """
+    `CohortAnalysis` supports running a palette of genotype-phenotype correlation analyses for a :class:`Cohort`.
+
+    The genotype-phenotype correlation is run only for the present HPO terms. The excluded phenotypic features are used
+    to calculate the total number of the patients investigated for a term.
+
+    .. note:
+
+      `missing_implies_excluded` option assumes that a patients were investigated for all cohort phenotypic features,
+      and unless the feature is not explicitly present, it is marked as excluded and used in calculating
+      of the total number of the patients investigated for the feature.
 
-class CohortAnalysis():
-    """This class can be used to run different analyses of a given Cohort. 
 
     Attributes:
         analysis_cohort (Cohort): The Cohort object given for this analysis
         analysis_patients (list): A subset of patients with no structural variants if they are being removed from the analysis. Otherwise, all patients in the cohort.
         analysis_transcript (string): The given transcript ID that will be used in this analysis
-        is_recessive (boolean): True if the variant is recessive. Default is False. 
-        include_unmeasured (boolean): True if we want to assume a patient without a specific phenotype listed does not have that phenotype. Otherwise, will only include those that explicitely say the phenotype was not observed. Defaults to True. 
+        is_recessive (boolean): True if the variant is recessive. Default is False.
+        missing_implies_excluded (boolean): True if we assume that a patient without a specific phenotype listed
+          *does not* have the phenotype. Otherwise, the only excluded phenotypes are those that are excluded explicitly.
+          Defaults to `False`.
         include_large_SV (boolean): True if we want to include large structural variants in the analysis. Defaults to True.
-        min_perc_patients_w_hpo (integer): Will only test phenotypes that are listed in at least this percent of patients. Defaults to 10 (10%).
+        min_perc_patients_w_hpo (int | float): Will only test phenotypes that are listed in at least this percent of patients.
+          Can be a `float` in range :math:`(0,1]` or a positive `int` with the minimum count (included). Defaults to `.1` (10%).
     Methods:
         compare_by_variant_type(var_type1:VariantEffect, var_type2:Optional VariantEffect): Runs Fisher Exact analysis, finds any correlation between given variant effects and phenotypes
         compare_by_variant_type(variant1:string, variant2:Optional string): Runs Fisher Exact analysis, finds any correlation between given variants and phenotypes
@@ -35,78 +50,90 @@ class CohortAnalysis():
         compare_by_variant_type(feature1:FeatureType, feature2:Optional FeatureType): Runs Fisher Exact analysis, finds any correlation between given protein feature types and phenotypes
         compare_by_variant_type(feature1:string, feature2:Optional string): Runs Fisher Exact analysis, finds any correlation between given protein features and phenotypes
     """
-    def __init__(self, cohort, transcript, hpo, recessive = False, include_unmeasured = True, p_val_correction = 'bonferroni', 
-                include_large_SV = True, min_perc_patients_w_hpo = 10) -> None:
-        """Constructs all necessary attributes for a CohortAnalysis object 
-
-        Args:
-            cohort (Cohort): The Cohort object that will be analyzed
-            transcript (string): The transcript ID that will be used in this analysis
-            recessive (boolean): True if the variant is recessive. Default is False. 
-            p_val_correction (string): Defaults to bonferroni. Options - bonferroni, sidak, holm-sidak, holm, simes-hochberg, hommel, fdr_bh, fdr_by, fdr_tsbh, fdr_tsbky.
-            include_unmeasured (boolean): True if we want to assume a patient without a specific phenotype listed does not have that phenotype. Otherwise, will only include those that explicitely say the phenotype was not observed. Defaults to True. 
-            include_large_SV (boolean): True if we want to include large structural variants in the analysis. Defaults to True.
-            min_perc_patients_w_hpo (integer): Will only test phenotypes that are listed in at least this percent of patients. Defaults to 10 (10%).
-        """
+    def __init__(self, cohort: Cohort, transcript: str, hpo: hpotk.MinimalOntology, recessive: bool = False,
+                 missing_implies_excluded: bool = False, p_val_correction: str = 'bonferroni',
+                include_large_SV: bool = True, min_perc_patients_w_hpo: typing.Union[float, int] = .1) -> None:
         if not isinstance(cohort, Cohort):
             raise ValueError(f"cohort must be type Cohort but was type {type(cohort)}")
         if transcript not in cohort.all_transcripts:
             raise ValueError(f"Transcript {transcript} not found in Cohort")
+        self._logger = logging.getLogger(__name__)
         self._cohort = cohort
         self._transcript = transcript
+        self._hpo = hpotk.util.validate_instance(hpo, hpotk.MinimalOntology, 'hpo')
         self._recessive = recessive
-        self._in_unmeasured = include_unmeasured ##TODO: Figure out how to move this out of analysis 
+        self._missing_implies_excluded = missing_implies_excluded ##TODO: Figure out how to move this out of analysis
         self._include_SV = include_large_SV
-        self._patient_list = cohort.all_patients if include_large_SV else [pat for pat in cohort.all_patients if not all([var.variant_coordinates.is_structural() for var in pat.variants])]
-        if not isinstance(min_perc_patients_w_hpo, int) or min_perc_patients_w_hpo <= 0 or min_perc_patients_w_hpo > 100:
-            raise ValueError(f"hpo_percent_of_patients must be an integer between 1 & 100 but was {min_perc_patients_w_hpo}")
-        self._percent_pats = min_perc_patients_w_hpo
-        self._hpo_present_test = HPOPresentPredicate(hpo)
+        self._patient_list = list(cohort.all_patients) if include_large_SV else [pat for pat in cohort.all_patients if not all(var.variant_coordinates.is_structural() for var in pat.variants)]
+        if len(self._patient_list) == 0:
+            raise ValueError('No patients left for analysis!')
+        self._percent_pats = self._check_min_perc_patients_w_hpo(min_perc_patients_w_hpo, len(self._patient_list))
+
+        # As of now, `self._testing_hpo_terms` is a Sequence[TermId] with terms that were *present* in >=1 cohort member.
         self._testing_hpo_terms = self._remove_low_hpo_terms()
-        self._patients_by_hpo = self._sort_patients_by_hpo()
+        self._patients_by_hpo = self._group_patients_by_hpo(self._testing_hpo_terms, self._patient_list,
+                                                            self._hpo, self._missing_implies_excluded)
         self._correction = p_val_correction
-        self._logger = logging.getLogger(__name__)
-        
+
     @property
-    def analysis_cohort(self):
+    def analysis_cohort(self) -> Cohort:
         """
         Returns:
             Cohort: The Cohort object being analyzed
         """
+        # TODO[0.3.0] - remove
+        # This is because we may actually filter out some patients in the constructor and, in result,
+        # not analyze the ENTIRE cohort. Exposing the cohort is, thus, superfluous and unless there are strong reasons
+        # for keeping the cohort, I think it should be removed.
+        warnings.warn('The `analysis_cohort` will be removed from the CohortAnalyzer API in 0.3.0',
+                      DeprecationWarning, stacklevel=2)
         return self._cohort
-    
+
     @property
-    def analysis_patients(self):
+    def analysis_patients(self) -> typing.Sequence[Patient]:
         """
         Returns:
-            list: A subset of patients with no structural variants if they are being removed from the analysis. Otherwise, all patients in the cohort.
+            Sequence: A sequence of patients that are subject to the analysis.
         """
         return self._patient_list
-    
+
     @property
-    def analysis_transcript(self):
+    def analysis_transcript(self) -> str:
         """
         Returns:
-            string: The given transcript ID that will be used in this analysis
+            string: The accession ID of the transcript that will be used in this analysis.
         """
         return self._transcript
 
     @property
-    def is_recessive(self):
+    def is_recessive(self) -> bool:
         """
         Returns:
-            boolean: True if the variant is recessive.
+            boolean: True if the analysis is assuming a disease with the Autosomal recessive mode of inheritance.
         """
         return self._recessive
 
     @property
-    def include_unmeasured(self):
+    def include_unmeasured(self) -> bool:
         """
         Returns:
-            boolean: True if we want to assume a patient without a specific phenotype listed does not have that phenotype. 
-            Otherwise, will only include those that explicitely say the phenotype was not observed. 
+            boolean: True if we want to assume a patient without a specific phenotype listed does not have that phenotype.
+            Otherwise, will only include those that explicitly say the phenotype was not observed.
         """
-        return self._in_unmeasured
+        # TODO[0.3.0] - remove
+        warnings.warn('The `include_unmeasured` will be removed from the CohortAnalyzer API in `v0.3.0`.'
+                      'Use `missing_implies_excluded` instead',
+                      DeprecationWarning, stacklevel=2)
+        return self.missing_implies_excluded
+
+    @property
+    def missing_implies_excluded(self) -> bool:
+        """
+        Returns:
+            boolean: True if we assume that patient without a specific phenotype listed *does not* have the phenotype.
+            Otherwise, the only excluded phenotypes are those that are excluded explicitly.
+        """
+        return self._missing_implies_excluded
 
     @property
     def include_large_SV(self):
@@ -120,18 +147,21 @@ def include_large_SV(self):
     def min_perc_patients_w_hpo(self):
         """
         Returns:
-            integer: Will only test phenotypes that are listed in at least this percent of patients. 
+            integer: Will only test phenotypes that are listed in at least this percent of patients.
         """
-        return self._percent_pats
+        # TODO[0.3.0] - remove
+        warnings.warn('The `min_perc_patients_w_hpo` will be removed from the CohortAnalyzer API in 0.3.0',
+                      DeprecationWarning, stacklevel=2)
+        return round(self._percent_pats * 100)
 
-    def _remove_low_hpo_terms(self):
+    def _remove_low_hpo_terms(self) -> typing.Sequence[hpotk.TermId]:
         final_hpo_list = set()
         all_hpo_counts = Counter()
         all_no_hpo_counts = Counter()
         for pat in self.analysis_patients:
-            all_hpo_counts.update([hpo for hpo in pat.phenotypes if hpo.observed == True])
+            all_hpo_counts.update(pat.present_phenotypes())
             if not self.include_unmeasured:
-                all_no_hpo_counts.update([hpo for hpo in pat.phenotypes if hpo.observed == False])
+                all_no_hpo_counts.update(pat.excluded_phenotypes())
         if self.include_unmeasured:
             for hpo, count in all_hpo_counts.items():
                 count = count if not None else 0
@@ -145,24 +175,10 @@ def _remove_low_hpo_terms(self):
                     final_hpo_list.add(hpo)
         if len(final_hpo_list) == 0:
             raise ValueError(f"No HPO terms found in over {self.min_perc_patients_w_hpo}% of patients.")
-        return final_hpo_list
-
-    def _sort_patients_by_hpo(self):
-        all_with_hpo = {}
-        all_without_hpo = {}
-        for hpo in self._testing_hpo_terms:
-            if self.include_unmeasured:
-                with_hpo = [pat for pat in self.analysis_patients if self._hpo_present_test.test(pat, hpo.identifier) == HPOPresentPredicate.OBSERVED]
-                not_hpo = [pat for pat in self.analysis_patients if self._hpo_present_test.test(pat, hpo.identifier) == HPOPresentPredicate.NOT_MEASURED] ## and not observed 
-            else:
-                with_hpo = [pat for pat in self.analysis_patients if self._hpo_present_test.test(pat, hpo.identifier) == HPOPresentPredicate.OBSERVED]
-                not_hpo = [pat for pat in self.analysis_patients if self._hpo_present_test.test(pat, hpo.identifier) == HPOPresentPredicate.NOT_OBSERVED]
-            all_with_hpo[hpo] = with_hpo
-            all_without_hpo[hpo] = not_hpo
-        patientsByHPO = namedtuple('patientByHPO', field_names=['all_with_hpo', 'all_without_hpo'])
-        return patientsByHPO(all_with_hpo, all_without_hpo)
-
-    def _run_dom_analysis(self, predicate, variable1, variable2):
+
+        return [p.identifier for p in final_hpo_list]
+
+    def _run_dom_analysis(self, predicate: PolyPredicate, variable1, variable2):
         final_dict = dict()
         all_pvals = []
         if not self.is_recessive:
@@ -176,11 +192,12 @@ def _run_dom_analysis(self, predicate, variable1, variable2):
                     with_hpo_var2_count = len([pat for pat in self._patients_by_hpo.all_with_hpo.get(hpo) if predicate.test(pat, variable2) in (HOMOZYGOUS, HETEROZYGOUS)])
                     not_hpo_var2_count = len([pat for pat in self._patients_by_hpo.all_without_hpo.get(hpo) if predicate.test(pat, variable2) in (HOMOZYGOUS, HETEROZYGOUS)])
                 if with_hpo_var1_count + not_hpo_var1_count == 0 or with_hpo_var2_count + not_hpo_var2_count == 0:
-                    self._logger.warning(f"Divide by 0 error with HPO {hpo.identifier.value}, not included in this analysis.")
+                    self._logger.warning(f"Divide by 0 error with HPO {hpo.value}, not included in this analysis.")
                     continue
-                p_val = self._run_fisher_exact([[with_hpo_var1_count, not_hpo_var1_count], [with_hpo_var2_count, not_hpo_var2_count]])
+                p_val = run_fisher_exact([[with_hpo_var1_count, not_hpo_var1_count], [with_hpo_var2_count, not_hpo_var2_count]])
                 all_pvals.append(p_val)
-                final_dict[f"{hpo.identifier.value} ({hpo.name})"] = [with_hpo_var1_count, "{:0.2f}%".format((with_hpo_var1_count/(with_hpo_var1_count + not_hpo_var1_count)) * 100), with_hpo_var2_count, "{:0.2f}%".format((with_hpo_var2_count/(with_hpo_var2_count + not_hpo_var2_count)) * 100), p_val]
+                term_name = self._hpo.get_term(hpo).name
+                final_dict[f"{hpo.value} ({term_name})"] = [with_hpo_var1_count, "{:0.2f}%".format((with_hpo_var1_count/(with_hpo_var1_count + not_hpo_var1_count)) * 100), with_hpo_var2_count, "{:0.2f}%".format((with_hpo_var2_count/(with_hpo_var2_count + not_hpo_var2_count)) * 100), p_val]
         else:
             return ValueError(f"Run a recessive analysis")
         corrected_pvals = multitest.multipletests(all_pvals, method = self._correction)[1]
@@ -206,27 +223,35 @@ def _run_rec_analysis(self, predicate, variable1, variable2):
                     with_hpo_var2_var2 = len([pat for pat in self._patients_by_hpo.all_with_hpo.get(hpo) if predicate.test(pat, variable2) == HOMOZYGOUS])
                     no_hpo_var2_var2 = len([pat for pat in self._patients_by_hpo.all_without_hpo.get(hpo) if predicate.test(pat, variable2) == HOMOZYGOUS])
                 if with_hpo_var1_var1 + no_hpo_var1_var1 == 0 or with_hpo_var1_var2 + no_hpo_var1_var2 == 0 or with_hpo_var2_var2 + no_hpo_var2_var2 == 0:
-                    self._logger.warning(f"Divide by 0 error with HPO {hpo.identifier.value}, not included in this analysis.")
+                    self._logger.warning(f"Divide by 0 error with HPO {hpo.value}, not included in this analysis.")
                     continue
-                p_val = self._run_recessive_fisher_exact([[with_hpo_var1_var1, no_hpo_var1_var1], [with_hpo_var1_var2, no_hpo_var1_var2], [with_hpo_var2_var2, no_hpo_var2_var2]] )
+                p_val = run_recessive_fisher_exact([[with_hpo_var1_var1, no_hpo_var1_var1], [with_hpo_var1_var2, no_hpo_var1_var2], [with_hpo_var2_var2, no_hpo_var2_var2]] )
                 all_pvals.append(p_val)
-                final_dict[f"{hpo.identifier.value} ({hpo.name})"] = [with_hpo_var1_var1, "{:0.2f}%".format((with_hpo_var1_var1/(with_hpo_var1_var1 + no_hpo_var1_var1)) * 100), with_hpo_var1_var2, "{:0.2f}%".format((with_hpo_var1_var2/(no_hpo_var1_var2 + with_hpo_var1_var2)) * 100), with_hpo_var2_var2, "{:0.2f}%".format((with_hpo_var2_var2/(with_hpo_var2_var2 + no_hpo_var2_var2)) * 100), p_val] 
+                term_name = self._hpo.get_term(hpo).name
+                final_dict[f"{hpo.value} ({term_name})"] = [with_hpo_var1_var1, "{:0.2f}%".format((with_hpo_var1_var1/(with_hpo_var1_var1 + no_hpo_var1_var1)) * 100), with_hpo_var1_var2, "{:0.2f}%".format((with_hpo_var1_var2/(no_hpo_var1_var2 + with_hpo_var1_var2)) * 100), with_hpo_var2_var2, "{:0.2f}%".format((with_hpo_var2_var2/(with_hpo_var2_var2 + no_hpo_var2_var2)) * 100), p_val]
         else:
             return ValueError("Run a dominant analysis")
         corrected_pvals = multitest.multipletests(all_pvals, method = self._correction)[1]
         return final_dict, corrected_pvals
 
+    def compare_by_variant_effect(self, effect: VariantEffect, tx_id: str):
+        # We need to cache the previous transcript
+        previous = self._transcript
+        self._transcript = tx_id
+        result = self.compare_by_variant_type(effect)
+        self._transcript = previous
+        return result
 
     def compare_by_variant_type(self, var_type1:VariantEffect, var_type2:VariantEffect = None):
         """Runs Fisher Exact analysis, finds any correlation between given variant effects across phenotypes.
-        
+
         Args:
-            var_type1 (VariantEffect): 
+            var_type1 (VariantEffect):
             var_type2 (VariantEffect, Optional): If None, we compare between those with var_type1 and those without var_type1
         Returns:
             DataFrame: A pandas DataFrame showing the results of the analysis
         """
-        var_effect_test = VariantEffectPredicate(self.analysis_transcript)
+        var_effect_test = VariantEffectPredicate(self._transcript, var_type1)
         if self.is_recessive:
             final_dict, corrected_pvals = self._run_rec_analysis(var_effect_test, var_type1, var_type2)
             if var_type2 is None:
@@ -249,11 +274,14 @@ def compare_by_variant_type(self, var_type1:VariantEffect, var_type2:VariantEffe
             final_df.insert(5, ('', "Corrected p-values"), corrected_pvals, True)
         return final_df.sort_values([('', 'Corrected p-values'), ('', 'p-value')])
 
+    def compare_by_variant_key(self, variant_key: str):
+        return self.compare_by_variant(variant_key)
+
     def compare_by_variant(self, variant1:str, variant2:str = None):
         """Runs Fisher Exact analysis, finds any correlation between given variants across phenotypes.
-        
+
         Args:
-            variant1 (string): 
+            variant1 (string):
             variant2 (string, Optional): If None, we compare between those with variant1 and those without variant1
         Returns:
             DataFrame: A pandas DataFrame showing the results of the analysis
@@ -271,7 +299,7 @@ def compare_by_variant(self, variant1:str, variant2:str = None):
             final_df = DataFrame.from_dict(final_dict, orient='index', columns=index.insert(6, ('', 'p-value')))
             final_df.insert(7, ('', "Corrected p-values"), corrected_pvals, True)
         else:
-            final_dict, corrected_pvals = self._run_dom_analysis(variant_test, variant1, variant2)  
+            final_dict, corrected_pvals = self._run_dom_analysis(variant_test, variant1, variant2)
             if variant2 is None:
                 col_name = f'Without {variant1}'
             else:
@@ -283,14 +311,14 @@ def compare_by_variant(self, variant1:str, variant2:str = None):
 
     def compare_by_exon(self, exon1:int, exon2:int = None):
         """Runs Fisher Exact analysis, finds any correlation between given exons across phenotypes.
-        
+
         Args:
-            exon1 (integer): 
+            exon1 (integer):
             exon2 (integer, Optional): If None, we compare between those within exon1 and those outside exon1
         Returns:
             DataFrame: A pandas DataFrame showing the results of the analysis
         """
-        exon_test = ExonPredicate(self.analysis_transcript)
+        exon_test = ExonPredicate(self.analysis_transcript, exon1)
         if self.is_recessive:
             final_dict, corrected_pvals = self._run_rec_analysis(exon_test, exon1, exon2)
             if exon2 is None:
@@ -302,7 +330,7 @@ def compare_by_exon(self, exon1:int, exon2:int = None):
             index = MultiIndex.from_product([[f'Homozygous {exon1}', col_name1, col_name2], ['Count', 'Percent']])
             final_df = DataFrame.from_dict(final_dict, orient='index', columns=index.insert(6, ('', 'p-value')))
             final_df.insert(7, ('', "Corrected p-values"), corrected_pvals, True)
-        else: 
+        else:
             final_dict, corrected_pvals = self._run_dom_analysis(exon_test, exon1, exon2)
             if exon2 is None:
                 col_name = f'Outside Exon {exon1}'
@@ -315,14 +343,14 @@ def compare_by_exon(self, exon1:int, exon2:int = None):
 
     def compare_by_protein_feature_type(self, feature1:FeatureType, feature2:FeatureType = None):
         """Runs Fisher Exact analysis, finds any correlation between given feature type across phenotypes.
-        
+
         Args:
-            feature1 (FeatureType): 
+            feature1 (FeatureType):
             feature2 (FeatureType, Optional): If None, we compare between those with feature1 and those without feature1
         Returns:
             DataFrame: A pandas DataFrame showing the results of the analysis
         """
-        feat_type_test = ProtFeatureTypePredicate(self.analysis_transcript)
+        feat_type_test = ProtFeatureTypePredicate(self.analysis_transcript, feature1)
         if self.is_recessive:
             final_dict, corrected_pvals = self._run_rec_analysis(feat_type_test, feature1, feature2)
             if feature2 is None:
@@ -348,14 +376,14 @@ def compare_by_protein_feature_type(self, feature1:FeatureType, feature2:Feature
 
     def compare_by_protein_feature(self, feature1:str, feature2:str = None):
         """Runs Fisher Exact analysis, finds any correlation between given feature and phenotypes.
-        
+
         Args:
-            feature1 (string): 
+            feature1 (string):
             feature2 (string, Optional): If None, we compare between those within feature1 and those outside feature1
         Returns:
             DataFrame: A pandas DataFrame showing the results of the analysis
         """
-        domain_test = ProtFeaturePredicate(self.analysis_transcript)
+        domain_test = ProtFeaturePredicate(self.analysis_transcript, feature1)
         if self.is_recessive:
             final_dict, corrected_pvals = self._run_rec_analysis(domain_test, feature1, feature2)
             if feature2 is None:
@@ -377,148 +405,7 @@ def compare_by_protein_feature(self, feature1:str, feature2:str = None):
             final_df = DataFrame.from_dict(final_dict, orient='index', columns=index.insert(4, ('', 'p-value')))
             final_df.insert(5, ('', "Corrected p-values"), corrected_pvals, True)
         return final_df.sort_values([('', 'Corrected p-values'), ('', 'p-value')])
-    
-
-    def _run_fisher_exact(self, two_by_two_table: typing.Sequence[typing.Sequence[int]]):
-        oddsr, p = stats.fisher_exact(two_by_two_table, alternative='two-sided')
-        return p
-
-    def _run_recessive_fisher_exact(self, two_by_three_table: typing.Sequence[typing.Sequence[int]]):
-        a = np.array(two_by_three_table, dtype=np.int64)
-        test_class = PythonMultiFisherExact()
-        val = test_class.calculate(a)
-        return val
-
 
-class MultiFisherExact(metaclass=abc.ABCMeta):
 
-    @abc.abstractmethod
-    def calculate(self, a: np.ndarray) -> float:
-        """
-        :param a: a 2x3 int array with counts
-        :returns: a p value calculated with Fisher's exact test
-        """
-        pass
-
-    @staticmethod
-    def _check_input(a: np.ndarray):
-        if not isinstance(a, np.ndarray):
-            raise ValueError(f'Expected a numpy array but got {type(a)}')
-        if not a.shape == (3, 2):
-            raise ValueError(f'Shape of the array must be (3, 2) but got {a.shape}')
-        if np.array_equal(a, np.zeros_like(a)):
-            raise ValueError(f'Array is all zeros, cannot run analysis')
-        
-
-
-class PythonMultiFisherExact(MultiFisherExact):
-
-    def calculate(self, a: np.ndarray) -> float:
-        MultiFisherExact._check_input(a)
-        return self._fisher_exact(a)
-
-    def _fisher_exact(self, table):
-        row_sum = []
-        col_sum = []
-
-        for i in range(len(table)):
-            temp = 0
-            for j in range(len(table[0])):
-                temp += table[i][j]
-            row_sum.append(temp)
-
-        for j in range(len(table[0])):
-            temp = 0
-            for i in range(len(table)):
-                temp += table[i][j]
-            col_sum.append(temp)
-
-        mat = [[0] * len(col_sum)] * len(row_sum)
-        pos = (0, 0)
-
-        p_0 = 1
-
-        for x in row_sum:
-            p_0 *= math.factorial(x)
-        for y in col_sum:
-            p_0 *= math.factorial(y)
-
-        n = 0
-        for x in row_sum:
-            n += x
-        p_0 /= Decimal(math.factorial(n))
-
-        for i in range(len(table)):
-            for j in range(len(table[0])):
-                p_0 /= Decimal(math.factorial(table[i][j]))
-
-        p = [0]
-        self._dfs(mat, pos, row_sum, col_sum, p_0, p)
-
-        return float(p[0])
-
-    def _dfs(self, mat, pos, r_sum, c_sum, p_0, p):
-
-        (xx, yy) = pos
-        (r, c) = (len(r_sum), len(c_sum))
-
-        mat_new = []
-
-        for i in range(len(mat)):
-            temp = []
-            for j in range(len(mat[0])):
-                temp.append(mat[i][j])
-            mat_new.append(temp)
-
-        if xx == -1 and yy == -1:
-            for i in range(r - 1):
-                temp = r_sum[i]
-                for j in range(c - 1):
-                    temp -= mat_new[i][j]
-                mat_new[i][c - 1] = temp
-            for j in range(c - 1):
-                temp = c_sum[j]
-                for i in range(r - 1):
-                    temp -= mat_new[i][j]
-                mat_new[r - 1][j] = temp
-            temp = r_sum[r - 1]
-            for j in range(c - 1):
-                temp -= mat_new[r - 1][j]
-            if temp < 0:
-                return
-            mat_new[r - 1][c - 1] = temp
-
-            p_1 = 1
-            for x in r_sum:
-                p_1 *= math.factorial(x)
-            for y in c_sum:
-                p_1 *= math.factorial(y)
-
-            n = 0
-            for x in r_sum:
-                n += x
-            p_1 /= Decimal(math.factorial(n))
-
-            for i in range(len(mat_new)):
-                for j in range(len(mat_new[0])):
-                    p_1 /= Decimal(math.factorial(mat_new[i][j]))
-            if p_1 <= p_0 + Decimal(0.00000001):
-                # print(mat_new)
-                # print(p_1)
-                p[0] += p_1
-        else:
-            max_1 = r_sum[xx]
-            max_2 = c_sum[yy]
-            for j in range(c):
-                max_1 -= mat_new[xx][j]
-            for i in range(r):
-                max_2 -= mat_new[i][yy]
-            for k in range(min(max_1, max_2) + 1):
-                mat_new[xx][yy] = k
-                if xx == r - 2 and yy == c - 2:
-                    pos_new = (-1, -1)
-                elif xx == r - 2:
-                    pos_new = (0, yy + 1)
-                else:
-                    pos_new = (xx + 1, yy)
-                self._dfs(mat_new, pos_new, r_sum, c_sum, p_0, p)
+    def compare_by_variant_keys(self, a: str, b: str):
+        raise NotImplementedError()
diff --git a/src/genophenocorr/analysis/_api.py b/src/genophenocorr/analysis/_api.py
new file mode 100644
index 00000000..4a6fc75a
--- /dev/null
+++ b/src/genophenocorr/analysis/_api.py
@@ -0,0 +1,273 @@
+import abc
+import typing
+from collections import namedtuple, defaultdict
+
+import hpotk
+import pandas as pd
+
+from genophenocorr.model import VariantEffect, Patient
+from .predicate import PatientCategory
+
+PatientsByHPO = namedtuple('PatientsByHPO', field_names=['all_with_hpo', 'all_without_hpo'])
+
+
+class GenotypePhenotypeAnalysisResult:
+    """
+    `GenotypePhenotypeAnalysisResult` summarizes results of genotype-phenotype correlation analysis of a cohort.
+    """
+
+    def __init__(self, n_usable: pd.Series,
+                 all_counts: pd.DataFrame,
+                 pvals: pd.Series,
+                 corrected_pvals: typing.Optional[pd.Series],
+                 phenotype_categories: typing.Iterable[PatientCategory],
+                 question: str):
+        self._n_usable = n_usable
+        self._all_counts = all_counts
+        self._pvals = pvals
+        self._corrected_pvals = corrected_pvals
+        self._phenotype_categories = tuple(phenotype_categories)
+        self._question = question
+
+    @property
+    def n_usable(self) -> pd.Series:
+        """
+        Get a :class:`pandas.Series` with mapping from HPO term and number of patients where the term was either
+         present or excluded, and are, thus, usable for genotype-phenotype correlation analysis.
+        """
+        return self._n_usable
+
+    @property
+    def all_counts(self) -> pd.DataFrame:
+        """
+        Get a :class:`pandas.DataFrame` with counts of patients in genotype and phenotype groups.
+
+        An example for a genotype predicate that bins into two categories (`Yes` and `No`) based on presence
+        of a missense variant in transcript, and phenotype predicate that checks presence/absence of a phenotype term::
+
+                                   Has MISSENSE_VARIANT in NM_123456.7
+                                   No       Yes
+            HPO         Present
+            HP:0001166  Yes        1        13
+            HP:0001166  No         7        5
+            HP:0003011  Yes        3        11
+            HP:0003011  No         7        16
+            HP:0012638  Yes        9        15
+            HP:0012638  No         3        4
+
+        """
+        return self._all_counts
+
+    @property
+    def pvals(self) -> pd.Series:
+        """
+        Get a :class:`pandas.Series` with p values for each tested HPO term.
+        """
+        return self._pvals
+
+    @property
+    def corrected_pvals(self) -> typing.Optional[pd.Series]:
+        """
+        Get an optional :class:`pandas.Series` with p values for each tested HPO term after multiple testing correction.
+        """
+        return self._corrected_pvals
+
+    @property
+    def phenotype_categories(self) -> typing.Sequence[PatientCategory]:
+        """
+        Get a sequence of phenotype patient categories that can be investigated.
+        """
+        return self._phenotype_categories
+
+    def summarize(self, hpo: hpotk.MinimalOntology,
+                  phenotype_category: PatientCategory) -> pd.DataFrame:
+        """
+        Create a data frame with summary of the genotype phenotype analysis.
+
+        The *rows* of the frame correspond to the analyzed HPO terms.
+
+        The columns of the data frame have `Count` and `Percentage` per used genotype predicate.
+        For instance, if we used :class:`genophenocorr.analysis.predicate.genotype.VariantEffectPredicate`
+        which compares phenotype with and without variants with certain :class:`genophenocorr.model.VariantEffect`,
+        we will have the following columns:
+
+        =====  =======  =====  =======
+        Yes             No
+        --------------  --------------
+        Count  Percent  Count  Percent
+        =====  =======  =====  =======
+
+        The final data frame may look something like this::
+
+            MISSENSE_VARIANT on `NM_1234.5`                       No                  Yes
+                                                                  Count    Percent    Count    Percent   p value   Corrected p value
+            Arachnodactyly [HP:0001166]                           1         3.8461    13       50.0000   0.000781  0.020299
+            Abnormality of the musculature [HP:0003011]           6        35.2941    11       64.7058   1.000000  1.000000
+            Abnormal nervous system physiology [HP:0012638]       9        37.5000    15       62.5000   1.000000  1.000000
+            ...                                                   ...      ...        ...      ...       ...       ...
+        """
+        if phenotype_category not in self._phenotype_categories:
+            raise ValueError(f'Unknown phenotype category: {phenotype_category}. Use one of {self._phenotype_categories}')
+
+        # Row index: a list of tested HPO terms
+        pheno_idx = pd.Index(self._n_usable.index)
+        # Column index: multiindex of counts and percentages for all genotype predicate groups
+        geno_idx = pd.MultiIndex.from_product(
+            iterables=(self._all_counts.columns, ('Count', 'Percent')),
+            names=(self._question, None))
+
+        # We'll fill this frame with data
+        df = pd.DataFrame(index=pheno_idx, columns=geno_idx)
+
+        # We must choose the counts with the `phenotype_category` of interest.
+        # This depends on the used phenotype predicate.
+        # Usually, we are interested patients who HAVE a feature,
+        # but we can also investigate those who do NOT have the feature.
+        # Hence, we need `phenotype_category`.
+        counts = self._all_counts.loc[(slice(None), phenotype_category), :]
+        counts.set_index(pheno_idx, inplace=True)
+
+        # Fill the frame cells
+        for col in geno_idx.levels[0]:
+            cnt = counts[col]
+            df[col, 'Count'] = cnt
+            df[col, 'Percent'] = cnt * 100 / self._n_usable
+
+        # Add columns with p values and corrected p values (if present)
+        df.insert(df.shape[1], ('', self._pvals.name), self._pvals)
+        if self._corrected_pvals is not None:
+            df.insert(df.shape[1], ('', self._corrected_pvals.name), self._corrected_pvals)
+
+        # Format the index values: `HP:0001250` -> `Seizure [HP:0001250]`
+        labeled_idx = df.index.map(lambda term_id: f'{hpo.get_term(term_id).name} [{term_id.value}]')
+
+        # Last, sort by corrected p value or just p value
+        df = df.set_index(labeled_idx)
+        if self._corrected_pvals is not None:
+            return df.sort_values(by=('', self._corrected_pvals.name))
+        else:
+            return df.sort_values(by=('', self._pvals.name))
+
+
+class CohortAnalysis(metaclass=abc.ABCMeta):
+    """
+    `CohortAnalysis` is a driver class for running genotype-phenotype correlation analyses.
+
+    The class provides various methods to test genotype-phenotype correlations. All methods wrap results
+    into :class:`GenotypePhenotypeAnalysisResult`.
+    """
+
+    @abc.abstractmethod
+    def compare_by_variant_effect(self, effect: VariantEffect, tx_id: str) -> GenotypePhenotypeAnalysisResult:
+        """
+        Compares variants with `effect` vs. variants with all other variant effects.
+
+        :param effect: variant effect
+        :param tx_id: the accession of the transcript of interest
+        """
+        pass
+
+    @abc.abstractmethod
+    def compare_by_variant_key(self, variant_key: str) -> GenotypePhenotypeAnalysisResult:
+        """
+        Compares variant with `variant_key` vs all the other variants.
+
+        .. seealso::
+
+          :attr:`genophenocorr.model.VariantCoordinates.variant_key`
+
+        :param variant_key: a `str` with the variant key, e.g. ``X_12345_12345_C_G``
+        """
+        pass
+
+    @abc.abstractmethod
+    def compare_by_exon(self, exon_number: int, tx_id: str) -> GenotypePhenotypeAnalysisResult:
+        """
+        Compares variants with `effect` vs. variants with all other variant effects.
+
+        .. note::
+
+          We use 1-based numbering to number the exons, not the usual 0-based numbering of the computer science.
+          Therefore, the first exon of the transcript has ``exon_number==1``, the second exon is ``2``, and so on ...
+
+        .. note::
+
+          We do not check if the `exon_number` spans beyond the number of exons of the given `transcript_id`!
+          Therefore, ``exon_number==10,000`` will effectively return :attr:`BooleanPredicate.FALSE` for *all* patients!!! 😱
+          Well, at least the patients of the *Homo sapiens sapiens* taxon...
+
+        :param exon_number: a positive `int` with exon number
+        :param tx_id: the accession of the transcript of interest
+        """
+        pass
+
+    @abc.abstractmethod
+    def compare_by_variant_keys(self, a: str, b: str) -> GenotypePhenotypeAnalysisResult:
+        """
+        Compare genotype-phenotype correlation between groups with variant `a` vs. variant `b`.
+
+        .. seealso::
+
+          :attr:`genophenocorr.model.VariantCoordinates.variant_key`
+
+        :param a: a `str` with variant key of variant `a`, e.g. ``X_12345_12345_C_G``
+        :param b: a `str` with variant key of variant `b`
+        """
+        pass
+
+    @staticmethod
+    def _check_min_perc_patients_w_hpo(min_perc_patients_w_hpo: typing.Union[int, float],
+                                       cohort_size: int) -> float:
+        """
+        Check if the input meets the requirements.
+        """
+        if isinstance(min_perc_patients_w_hpo, int):
+            if min_perc_patients_w_hpo > 0:
+                return min_perc_patients_w_hpo / cohort_size
+            else:
+                raise ValueError(f'`min_perc_patients_w_hpo` must be a positive `int` '
+                                 f'but got {min_perc_patients_w_hpo}')
+        elif isinstance(min_perc_patients_w_hpo, float):
+            if 0 < min_perc_patients_w_hpo <= 1:
+                return min_perc_patients_w_hpo
+            else:
+                raise ValueError(f'`min_perc_patients_w_hpo` must be a `float` in range (0, 1] '
+                                 f'but got {min_perc_patients_w_hpo}')
+        else:
+            raise ValueError(f'`min_perc_patients_w_hpo` must be a positive `int` or a `float` in range (0, 1] '
+                             f'but got {type(min_perc_patients_w_hpo)}')
+
+    @staticmethod
+    def _group_patients_by_hpo(phenotypic_features: typing.Iterable[hpotk.TermId],
+                               patients: typing.Iterable[Patient],
+                               hpo: hpotk.GraphAware,
+                               missing_implies_excluded: bool) -> PatientsByHPO:
+        all_with_hpo = defaultdict(list)
+        all_without_hpo = defaultdict(list)
+        for hpo_term in phenotypic_features:
+            for patient in patients:
+                found = False
+                for pf in patient.present_phenotypes():
+                    if hpo_term == pf.identifier or hpo.graph.is_ancestor_of(hpo_term, pf):
+                        # Patient is annotated with `hpo_term` because `pf` is equal to `hpo_term`
+                        # or it is a descendant of `hpo_term`.
+                        all_with_hpo[hpo_term].append(patient)
+
+                        # If one `pf` of the patient is found to be a descendant of `hpo`, we must break to prevent
+                        # adding the patient to `present_hpo` more than once due to another descendant!
+                        found = True
+                        break
+                if not found:
+                    # The patient is not annotated by the `hpo_term`.
+
+                    if missing_implies_excluded:
+                        # The `hpo_term` annotation is missing, hence implicitly excluded.
+                        all_without_hpo[hpo_term].append(patient)
+                    else:
+                        # The `hpo_term` must be explicitly excluded patient to be accounted for.
+                        for ef in patient.excluded_phenotypes():
+                            if hpo_term == ef.identifier or hpo.graph.is_descendant_of(hpo_term, ef):
+                                all_with_hpo[hpo_term].append(patient)
+                                break
+
+        return PatientsByHPO(all_with_hpo, all_without_hpo)
diff --git a/src/genophenocorr/analysis/_commie.py b/src/genophenocorr/analysis/_commie.py
new file mode 100644
index 00000000..d4627cf0
--- /dev/null
+++ b/src/genophenocorr/analysis/_commie.py
@@ -0,0 +1,207 @@
+import logging
+import typing
+from collections import Counter
+
+import hpotk
+import pandas as pd
+
+from statsmodels.stats import multitest
+
+from genophenocorr.model import Patient, Cohort, VariantEffect
+from .predicate import BooleanPredicate, PolyPredicate
+from .predicate.genotype import VariantEffectPredicate, VariantPredicate, ExonPredicate
+from .predicate.phenotype import PropagatingPhenotypeBooleanPredicateFactory, PhenotypePredicateFactory
+
+from ._api import CohortAnalysis, GenotypePhenotypeAnalysisResult
+from ._stats import run_fisher_exact
+
+
+def _filter_rare_phenotypes_using_hierarchy(patients: typing.Collection[Patient],
+                                            min_perc_patients_w_hpo: float,
+                                            hpo: hpotk.GraphAware,
+                                            missing_implies_excluded: bool = False) -> typing.Sequence[hpotk.TermId]:
+    """
+    We want to get a corpus of HPO terms that are present in at least certain fraction of the subjects.
+    We will get the IDs of the present features only.
+
+    :param min_perc_patients_w_hpo: the minimum fraction of patients that *must* have the feature.
+      The `float` must be in range :math:`(0,1]`.
+    """
+    # TODO - review the logic with Lauren. Should we also do annotation propagation here?
+    #  What is the meaning of _include_unmeasured in this context?
+    present_count = Counter()
+    excluded_count = Counter()
+
+    for patient in patients:
+        for pf in patient.phenotypes:
+            if pf.is_observed:
+                # A present phenotypic feature must be counted in.
+                present_count[pf.identifier] += 1
+                # implies presence of its ancestors.
+                for anc in hpo.graph.get_ancestors(pf):
+                    present_count[anc] += 1
+            else:
+                # An excluded phenotypic feature
+                excluded_count[pf.identifier] += 1
+                for desc in hpo.graph.get_descendants(pf):
+                    # implies exclusion of its descendants.
+                    excluded_count[desc] += 1
+
+    if missing_implies_excluded:
+        # We must do another pass to ensure that we account the missing features as excluded.
+        # `present_count` has all phenotypic features that were observed as present among the cohort members.
+        for pf in present_count:
+            for patient in patients:
+                pf_can_be_inferred_as_excluded_for_patient = False
+                for phenotype in patient.phenotypes:
+                    if pf == phenotype.identifier:
+                        # Patient is explicitly annotated with `pf`. No inference for this patient!
+                        pf_can_be_inferred_as_excluded_for_patient = False
+                        break
+
+                    if phenotype.is_observed and hpo.graph.is_ancestor_of(pf, phenotype):
+                        # The `pf` is implicitly observed in the patient. No inference for this patient!
+                        pf_can_be_inferred_as_excluded_for_patient = False
+                        break
+                    elif not phenotype.is_observed and hpo.graph.is_descendant_of(pf, phenotype):
+                        # The `pf` is implicitly excluded in the patient. No inference for this patient!
+                        pf_can_be_inferred_as_excluded_for_patient = False
+                        break
+                    else:
+                        # The `pf` is observed or excluded neither implicitly nor explicitly.
+                        # We can infer that it is excluded!
+                        pf_can_be_inferred_as_excluded_for_patient = True
+
+                if pf_can_be_inferred_as_excluded_for_patient:
+                    excluded_count[pf] += 1
+                    for desc in hpo.graph.get_descendants(pf):
+                        excluded_count[desc] += 1
+
+    total_count = Counter()
+    for term_id, count in present_count.items():
+        total_count[term_id] += count
+
+    for term_id, count in excluded_count.items():
+        total_count[term_id] += count
+
+    final_hpo = []
+    for term_id, n_present in present_count.items():
+        n_all = total_count[term_id]
+        ratio = n_present / n_all
+        if ratio >= min_perc_patients_w_hpo:
+            final_hpo.append(term_id)
+
+    if len(final_hpo) == 0:
+        raise ValueError(f"No HPO terms found in over {min_perc_patients_w_hpo * 100}% of patients.")
+
+    return final_hpo
+
+
+
+class CommunistCohortAnalysis(CohortAnalysis):
+
+    def __init__(self, cohort: Cohort,
+                 hpo: hpotk.MinimalOntology,
+                 missing_implies_excluded: bool = False,
+                 include_sv: bool = False,
+                 p_val_correction: typing.Optional[str] = None,
+                 min_perc_patients_w_hpo: typing.Union[float, int] = .1):
+        if not isinstance(cohort, Cohort):
+            raise ValueError(f"cohort must be type Cohort but was type {type(cohort)}")
+
+        self._logger = logging.getLogger(__name__)
+        self._hpo = hpotk.util.validate_instance(hpo, hpotk.MinimalOntology, 'hpo')
+        self._phenotype_predicate_factory = PropagatingPhenotypeBooleanPredicateFactory(self._hpo,
+                                                                                        missing_implies_excluded)
+        self._correction = p_val_correction
+        self._patient_list = list(cohort.all_patients) \
+            if include_sv \
+            else [pat for pat in cohort.all_patients if not all(var.variant_coordinates.is_structural() for var in pat.variants)]
+        if len(self._patient_list) == 0:
+            raise ValueError('No patients left for analysis!')
+        self._percent_pats = self._check_min_perc_patients_w_hpo(min_perc_patients_w_hpo, len(self._patient_list))
+
+        # As of now, `self._testing_hpo_terms` is a Sequence[TermId] with terms that were *present* in >=1 cohort member.
+        self._testing_hpo_terms = _filter_rare_phenotypes_using_hierarchy(self._patient_list, self._percent_pats,
+                                                                          self._hpo, missing_implies_excluded)
+
+        self._patients_by_hpo = self._group_patients_by_hpo(self._testing_hpo_terms, self._patient_list,
+                                                            self._hpo, missing_implies_excluded)
+
+    def compare_by_variant_effect(self, effect: VariantEffect, tx_id: str) -> GenotypePhenotypeAnalysisResult:
+        predicate = VariantEffectPredicate(tx_id, effect)
+        return self._apply_boolean_predicate(predicate)
+
+    def compare_by_variant_key(self, variant_key: str) -> GenotypePhenotypeAnalysisResult:
+        predicate = VariantPredicate(variant_key)
+        return self._apply_boolean_predicate(predicate)
+
+    def compare_by_exon(self, exon_number: int, tx_id: str) -> GenotypePhenotypeAnalysisResult:
+        predicate = ExonPredicate(tx_id, exon_number)
+        return self._apply_boolean_predicate(predicate)
+
+    def compare_by_variant_keys(self, a: str, b: str) -> GenotypePhenotypeAnalysisResult:
+        raise NotImplementedError()
+
+    def _apply_boolean_predicate(self, predicate: BooleanPredicate) -> GenotypePhenotypeAnalysisResult:
+        return self._run_gp_analysis(self._patient_list, self._testing_hpo_terms,
+                                     self._phenotype_predicate_factory, predicate)
+
+    def _apply_poly_predicate(self, predicate: PolyPredicate) -> GenotypePhenotypeAnalysisResult:
+        return self._run_gp_analysis(self._patient_list, self._testing_hpo_terms,
+                                     self._phenotype_predicate_factory, predicate)
+
+    def _run_gp_analysis(self, patients: typing.Iterable[Patient],
+                         phenotypic_features: typing.Iterable[hpotk.TermId],
+                         pheno_predicate_factory: PhenotypePredicateFactory,
+                         geno_predicate: PolyPredicate) -> GenotypePhenotypeAnalysisResult:
+        # 1) Count the patients
+        n_usable, all_counts = self._count_patients(patients, phenotypic_features,
+                                                    pheno_predicate_factory, geno_predicate)
+
+        # 2) Statistical tests
+        pvals_idx = pd.Index(phenotypic_features, name='p_val')
+        pvals = pd.Series(float('nan'), index=pvals_idx, name='p value')
+        for pf in phenotypic_features:
+            counts = all_counts.loc[pf]
+            # TODO - this is where we must fail unless we have the contingency table of the right size!
+            pvals[pf] = run_fisher_exact(counts)
+
+        # 3) Multiple correction
+        if self._correction is not None:
+            _, pvals_corrected, _, _ = multitest.multipletests(pvals, alpha=.05, method=self._correction)
+            corrected_idx = pd.Index(phenotypic_features, name='p_val_corrected')
+            corrected_pvals_series = pd.Series(data=pvals_corrected, index=corrected_idx,
+                                               name='Corrected p value')
+        else:
+            corrected_pvals_series = None
+
+        return GenotypePhenotypeAnalysisResult(n_usable, all_counts, pvals, corrected_pvals_series,
+                                               pheno_predicate_factory.get_categories(), geno_predicate.get_question())
+
+    @staticmethod
+    def _count_patients(patients: typing.Iterable[Patient],
+                        phenotypic_features: typing.Iterable[hpotk.TermId],
+                        pheno_predicate_factory: PhenotypePredicateFactory,
+                        geno_predicate: PolyPredicate) -> typing.Tuple[pd.Series, pd.DataFrame]:
+        row_idx = pd.MultiIndex.from_product(
+            iterables=(phenotypic_features, pheno_predicate_factory.get_categories()),
+            names=('phenotypic_feature', 'category')
+        )
+        col_idx = pd.Index(geno_predicate.get_categories(), name=geno_predicate.get_question())
+        counts = pd.DataFrame(data=0, index=row_idx, columns=col_idx)
+        n_usable_patients = pd.Series(data=0, index=pd.Index(phenotypic_features))
+
+        # Apply genotype and phenotype predicates
+        for pf in phenotypic_features:
+            pheno_predicate = pheno_predicate_factory.get_predicate(pf)
+
+            for patient in patients:
+                pheno_cat = pheno_predicate.test(patient)
+                geno_cat = geno_predicate.test(patient)
+
+                if pheno_cat is not None and geno_cat is not None:
+                    counts.loc[(pf, pheno_cat), geno_cat] += 1
+                    n_usable_patients[pf] += 1
+
+        return n_usable_patients, counts
diff --git a/src/genophenocorr/analysis/_config.py b/src/genophenocorr/analysis/_config.py
new file mode 100644
index 00000000..11d4e0c5
--- /dev/null
+++ b/src/genophenocorr/analysis/_config.py
@@ -0,0 +1,161 @@
+import logging
+import typing
+
+import hpotk
+
+from genophenocorr.model import Cohort
+
+from ._api import CohortAnalysis
+from ._commie import CommunistCohortAnalysis
+
+
+class CohortAnalysisConfiguration:
+    """
+    `CohortAnalysisConfiguration` is a value class for storing :class:`genophenocorr.analysis.CohortAnalysis`
+    configuration options.
+
+    Default values
+    ^^^^^^^^^^^^^^
+
+    ==============================  ===========  ====================
+        Option                        Type           Default value
+    ==============================  ===========  ====================
+     ``missing_implies_excluded``    `bool`      `False`
+     ``pval_correction``             `str`       `bonferroni`
+     ``min_perc_patients_w_hpo``     `float`     `0.1`
+     ``include_sv``                  `bool`      `False`
+    ==============================  ===========  ====================
+
+    """
+
+    def __init__(self, missing_implies_excluded: bool,
+                 pval_correction: typing.Optional[str],
+                 min_perc_patients_w_hpo: float,
+                 include_sv: bool):
+        self._missing_implies_excluded = missing_implies_excluded
+        self._pval_correction = pval_correction
+        self._min_perc_patients_w_hpo = min_perc_patients_w_hpo
+        self._include_sv = include_sv
+
+    @staticmethod
+    def builder():
+        """
+        Create a new :class:`CohortAnalysisConfigurationBuilder` initialized with default configuration options.
+        """
+        return CohortAnalysisConfigurationBuilder()
+
+    @property
+    def missing_implies_excluded(self) -> bool:
+        """
+        `True` if we assume that a patient without a specific phenotype listed *does not* have the phenotype.
+        Otherwise, the only excluded phenotypes are those that are excluded explicitly.
+        """
+        return self._missing_implies_excluded
+
+    @property
+    def pval_correction(self) -> typing.Optional[str]:
+        """
+        Get method for multiple testing p value correction. Default: `bonferroni`.
+        """
+        return self._pval_correction
+
+    @property
+    def min_perc_patients_w_hpo(self) -> float:
+        """
+        A threshold for removing rare HPO terms, only the terms that are observed in at least this fraction of patients
+        will be retained for the analysis.
+        """
+        return self._min_perc_patients_w_hpo
+
+    @property
+    def include_sv(self) -> bool:
+        """
+        `True` if we want to include structural variants in the analysis
+        (i.e. the variants that use symbolic VCF notation).
+        """
+        return self._include_sv
+
+
+class CohortAnalysisConfigurationBuilder:
+    """
+    `CohortAnalysisConfigurationBuilder` builds :class:`CohortAnalysisConfiguration`. The configuration options
+    can be customized or set with default values.
+
+    If an invalid option is passed, a warning is logged and the previous value is retained.
+    Consequently, it is *impossible* to mis-configure the analysis.
+    """
+
+    def __init__(self):
+        self._logger = logging.getLogger(__name__)
+        self._missing_implies_excluded = False
+        self._pval_correction = 'bonferroni'
+        self._min_perc_patients_w_hpo = .1
+        self._include_sv = False
+
+    def missing_implies_excluded(self, missing_implies_excluded: bool):
+        """
+        Set `missing_implies_excluded` option.
+        """
+        if isinstance(missing_implies_excluded, bool):
+            self._missing_implies_excluded = missing_implies_excluded
+        else:
+            self._logger.warning('Ignoring invalid `missing_implies_excluded` value %s',
+                                 missing_implies_excluded)
+
+        return self
+
+    def pval_correction(self, pval_correction: typing.Optional[str]):
+        """
+        Set `pval_correction` option.
+        """
+        # TODO - check admissible values
+        self._pval_correction = pval_correction
+        return self
+
+    def min_perc_patients_w_hpo(self, min_perc_patients_w_hpo: float):
+        """
+        Set `min_perc_patients_w_hpo` option.
+        """
+        # TODO - check float in range [0,1)
+        self._min_perc_patients_w_hpo = min_perc_patients_w_hpo
+        return self
+
+    def include_sv(self, include_sv: bool):
+        """
+        Set `include_sv` option.
+        """
+        if isinstance(include_sv, bool):
+            self._include_sv = include_sv
+        else:
+            self._logger.warning('Ignoring invalid `include_sv` value %s', include_sv)
+        return self
+
+    def build(self) -> CohortAnalysisConfiguration:
+        """
+        Build the configuration.
+        """
+        return CohortAnalysisConfiguration(self._missing_implies_excluded,
+                                           self._pval_correction,
+                                           self._min_perc_patients_w_hpo,
+                                           self._include_sv)
+
+
+def configure_cohort_analysis(cohort: Cohort,
+                              hpo: hpotk.MinimalOntology,
+                              config: typing.Optional[CohortAnalysisConfiguration] = None) -> CohortAnalysis:
+    """
+    Configure :class:`genophenocorr.analysis.CohortAnalysis` for given `cohort`.
+
+    :param cohort: a :class:`genophenocorr.model.Cohort` to analyze
+    :param hpo: a :class:`hpotk.MinimalOntology` with HPO to use in the analysis
+    :param config: an optional :class:`CohortAnalysisConfiguration` to parameterize the analysis.
+     The default parameters will be used if `None`.
+    """
+    if config is None:
+        config = CohortAnalysisConfiguration.builder().build()
+
+    return CommunistCohortAnalysis(cohort, hpo,
+                                   missing_implies_excluded=config.missing_implies_excluded,
+                                   include_sv=config.include_sv,
+                                   p_val_correction=config.pval_correction,
+                                   min_perc_patients_w_hpo=config.min_perc_patients_w_hpo)
diff --git a/src/genophenocorr/analysis/_stats.py b/src/genophenocorr/analysis/_stats.py
new file mode 100644
index 00000000..e5203388
--- /dev/null
+++ b/src/genophenocorr/analysis/_stats.py
@@ -0,0 +1,152 @@
+import abc
+import math
+import typing
+from decimal import Decimal
+
+import numpy as np
+import scipy
+
+
+class MultiFisherExact(metaclass=abc.ABCMeta):
+
+    @abc.abstractmethod
+    def calculate(self, a: np.ndarray) -> float:
+        """
+        :param a: a 2x3 int array with counts
+        :returns: a p value calculated with Fisher's exact test
+        """
+        pass
+
+    @staticmethod
+    def _check_input(a: np.ndarray):
+        if not isinstance(a, np.ndarray):
+            raise ValueError(f'Expected a numpy array but got {type(a)}')
+        if not a.shape == (3, 2):
+            raise ValueError(f'Shape of the array must be (3, 2) but got {a.shape}')
+        if np.array_equal(a, np.zeros_like(a)):
+            raise ValueError(f'Array is all zeros, cannot run analysis')
+
+
+class PythonMultiFisherExact(MultiFisherExact):
+
+    def calculate(self, a: np.ndarray) -> float:
+        MultiFisherExact._check_input(a)
+        return self._fisher_exact(a)
+
+    def _fisher_exact(self, table):
+        row_sum = []
+        col_sum = []
+
+        for i in range(len(table)):
+            temp = 0
+            for j in range(len(table[0])):
+                temp += table[i][j]
+            row_sum.append(temp)
+
+        for j in range(len(table[0])):
+            temp = 0
+            for i in range(len(table)):
+                temp += table[i][j]
+            col_sum.append(temp)
+
+        mat = [[0] * len(col_sum)] * len(row_sum)
+        pos = (0, 0)
+
+        p_0 = 1
+
+        for x in row_sum:
+            p_0 *= math.factorial(x)
+        for y in col_sum:
+            p_0 *= math.factorial(y)
+
+        n = 0
+        for x in row_sum:
+            n += x
+        p_0 /= Decimal(math.factorial(n))
+
+        for i in range(len(table)):
+            for j in range(len(table[0])):
+                p_0 /= Decimal(math.factorial(table[i][j]))
+
+        p = [0]
+        self._dfs(mat, pos, row_sum, col_sum, p_0, p)
+
+        return float(p[0])
+
+    def _dfs(self, mat, pos, r_sum, c_sum, p_0, p):
+
+        (xx, yy) = pos
+        (r, c) = (len(r_sum), len(c_sum))
+
+        mat_new = []
+
+        for i in range(len(mat)):
+            temp = []
+            for j in range(len(mat[0])):
+                temp.append(mat[i][j])
+            mat_new.append(temp)
+
+        if xx == -1 and yy == -1:
+            for i in range(r - 1):
+                temp = r_sum[i]
+                for j in range(c - 1):
+                    temp -= mat_new[i][j]
+                mat_new[i][c - 1] = temp
+            for j in range(c - 1):
+                temp = c_sum[j]
+                for i in range(r - 1):
+                    temp -= mat_new[i][j]
+                mat_new[r - 1][j] = temp
+            temp = r_sum[r - 1]
+            for j in range(c - 1):
+                temp -= mat_new[r - 1][j]
+            if temp < 0:
+                return
+            mat_new[r - 1][c - 1] = temp
+
+            p_1 = 1
+            for x in r_sum:
+                p_1 *= math.factorial(x)
+            for y in c_sum:
+                p_1 *= math.factorial(y)
+
+            n = 0
+            for x in r_sum:
+                n += x
+            p_1 /= Decimal(math.factorial(n))
+
+            for i in range(len(mat_new)):
+                for j in range(len(mat_new[0])):
+                    p_1 /= Decimal(math.factorial(mat_new[i][j]))
+            if p_1 <= p_0 + Decimal(0.00000001):
+                # print(mat_new)
+                # print(p_1)
+                p[0] += p_1
+        else:
+            max_1 = r_sum[xx]
+            max_2 = c_sum[yy]
+            for j in range(c):
+                max_1 -= mat_new[xx][j]
+            for i in range(r):
+                max_2 -= mat_new[i][yy]
+            for k in range(min(max_1, max_2) + 1):
+                mat_new[xx][yy] = k
+                if xx == r - 2 and yy == c - 2:
+                    pos_new = (-1, -1)
+                elif xx == r - 2:
+                    pos_new = (0, yy + 1)
+                else:
+                    pos_new = (xx + 1, yy)
+                self._dfs(mat_new, pos_new, r_sum, c_sum, p_0, p)
+
+
+def run_recessive_fisher_exact(two_by_three_table: typing.Sequence[typing.Sequence[int]]):
+    a = np.array(two_by_three_table, dtype=np.int64)
+    test_class = PythonMultiFisherExact()
+    val = test_class.calculate(a)
+    return val
+
+
+def run_fisher_exact(two_by_two_table: typing.Sequence[typing.Sequence[int]]):
+    oddsr, p = scipy.stats.fisher_exact(two_by_two_table, alternative='two-sided')
+    return p
diff --git a/src/genophenocorr/analysis/predicate/__init__.py b/src/genophenocorr/analysis/predicate/__init__.py
index 66593023..52ee0a37 100644
--- a/src/genophenocorr/analysis/predicate/__init__.py
+++ b/src/genophenocorr/analysis/predicate/__init__.py
@@ -1,9 +1,8 @@
-from ._api import SimplePredicate, PatientCategory, PolyPredicate
-from ._all_predicates import HPOPresentPredicate, VariantEffectPredicate, VariantPredicate, ExonPredicate, ProtFeaturePredicate, ProtFeatureTypePredicate
-from ._all_predicates import HOMOZYGOUS, HETEROZYGOUS, NO_VARIANT
+from . import genotype
+from . import phenotype
+
+from ._api import BooleanPredicate, PatientCategory, PolyPredicate
 
 __all__ = [
-    'SimplePredicate', 'PatientCategory', 'PolyPredicate',
-    'HPOPresentPredicate', 'VariantEffectPredicate', 'VariantPredicate', 'ExonPredicate', 'ProtFeaturePredicate', 'ProtFeatureTypePredicate',
-    'HOMOZYGOUS', 'HETEROZYGOUS', 'NO_VARIANT'
+    'PolyPredicate', 'PatientCategory', 'BooleanPredicate',
 ]
diff --git a/src/genophenocorr/analysis/predicate/_all_predicates.py b/src/genophenocorr/analysis/predicate/_all_predicates.py
deleted file mode 100644
index 3aa4b5c3..00000000
--- a/src/genophenocorr/analysis/predicate/_all_predicates.py
+++ /dev/null
@@ -1,273 +0,0 @@
-import typing
-
-import hpotk
-
-from genophenocorr.model import Patient, FeatureType, Genotype, VariantEffect
-from ._api import PolyPredicate, PatientCategory
-
-
-class HPOPresentPredicate(PolyPredicate[hpotk.TermId]):
-
-    OBSERVED = PatientCategory(cat_id=0,
-                               name='Observed',
-                               description="""
-                               The sample *is* annotated with the tested phenotype feature `q`.
-
-                               This is either because the sample is annotated with `q` (exact match),
-                               or because one of sample's annotations is a descendant `q` (annotation propagation).
-                               For instance, we tested for a Seizure and the sample *had* a Clonic seizure
-                               (a descendant of Seizure).
-                               """)
-
-    NOT_OBSERVED = PatientCategory(cat_id=1,
-                                   name='Not observed',
-                                   description="""
-                                   We are particular about the sample *not* having the tested feature `q`.
-
-                                   In other words, `q` was *excluded* in the sample or the sample is annotated with an excluded ancestor of `q`.
-
-                                   For instance, we tested for a Clonic seizure and the sample did *not* have any Seizure, which implies
-                                   *not* Clonic seizure.
-                                   """)
-
-    NOT_MEASURED = PatientCategory(cat_id=2,
-                                   name='Not measured',
-                                   description="""
-                                   We do not know if the sample has or has not the tested feature.
-                                   """)
-
-    def __init__(self,
-                 hpo: hpotk.MinimalOntology) -> None:
-        self._hpo = hpotk.util.validate_instance(hpo, hpotk.MinimalOntology, 'hpo')
-
-    def categories(self) -> typing.Sequence[PatientCategory]:
-        return self.OBSERVED, self.NOT_OBSERVED, self.NOT_MEASURED
-
-    def test(self, patient: Patient, query: hpotk.TermId) -> typing.Optional[PatientCategory]:
-        if not isinstance(patient, Patient):
-            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
-        query = hpotk.util.validate_instance(query, hpotk.TermId, 'query')
-
-        if len(patient.phenotypes) == 0:
-            return None
-
-        for pheno in patient.phenotypes:
-            if pheno.observed is None:
-                continue
-            if pheno.observed:
-                if any((query == ancestor
-                        for ancestor in self._hpo.graph.get_ancestors(pheno, include_source=True))):
-                    return self.OBSERVED
-            else:
-                if any((query == descendant
-                        for descendant in self._hpo.graph.get_descendants(pheno, include_source=True))):
-                    return self.NOT_OBSERVED
-
-        return self.NOT_MEASURED
-
-
-HETEROZYGOUS = PatientCategory(cat_id=0,
-                            name='Heterozygous',
-                            description="""
-                            This sample has the tested attribute on one allele.
-                            """)
-
-HOMOZYGOUS = PatientCategory(cat_id=1,
-                                name='Homozygous',
-                                description="""
-                                This sample has the tested attribute on both alleles.
-                                """)
-
-NO_VARIANT = PatientCategory(cat_id=2,
-                                name='No Variant',
-                                description="""
-                                The sample does not have the tested attribute.
-                                """)
-
-
-
-class VariantEffectPredicate(PolyPredicate[VariantEffect]):
-
-    def __init__(self, transcript:str) -> None:
-        self._transcript = transcript
-
-    def categories(self) -> typing.Sequence[PatientCategory]:
-        return HETEROZYGOUS, HOMOZYGOUS, NO_VARIANT
-
-    def test(self, patient: Patient, query: VariantEffect) -> typing.Optional[PatientCategory]:
-        if not isinstance(patient, Patient):
-            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
-        if not isinstance(query, VariantEffect):
-            raise ValueError(f"query must be type VariantEffect but was type {type(query)}")
-        vars = set()
-        for var in patient.variants:
-            for trans in var.tx_annotations:
-                if trans.transcript_id == self._transcript:
-                    for var_eff in trans.variant_effects:
-                        if var_eff == query:
-                            vars.add(var)
-        if len(vars) == 1:
-            for v in vars:
-                gt = v.genotype_for_sample(patient.patient_id)
-                if gt == Genotype.HETEROZYGOUS:
-                    return HETEROZYGOUS
-                elif gt == Genotype.HOMOZYGOUS_ALTERNATE:
-                    return HOMOZYGOUS
-                else:
-                    # TODO - is this really what we want to return here?
-                    return HETEROZYGOUS
-        elif len(vars) > 1:
-            return HOMOZYGOUS
-        else:
-            return NO_VARIANT
-
-class VariantPredicate(PolyPredicate[str]):
-
-    def __init__(self, transcript:str) -> None:
-        self._transcript = transcript
-
-    def categories(self) -> typing.Sequence[PatientCategory]:
-        return HETEROZYGOUS, HOMOZYGOUS, NO_VARIANT
-
-    def test(self, patient: Patient, query: str) -> typing.Optional[PatientCategory]:
-        if not isinstance(patient, Patient):
-            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
-        if not isinstance(query, str):
-            raise ValueError(f"query must be type string but was type {type(query)}")
-        vars = set()
-        for var in patient.variants:
-            if var.variant_coordinates.variant_key == query:
-                vars.add(var)
-        if len(vars) == 1:
-            for v in vars:
-                gt = v.genotype_for_sample(patient.patient_id)
-                if gt == Genotype.HETEROZYGOUS:
-                    return HETEROZYGOUS
-                elif gt == Genotype.HOMOZYGOUS_ALTERNATE:
-                    return HOMOZYGOUS
-                else:
-                    # TODO - is this really what we want to return here?
-                    return HETEROZYGOUS
-        elif len(vars) > 1:
-            return HOMOZYGOUS
-        else:
-            return NO_VARIANT
-
-class ExonPredicate(PolyPredicate[int]):
-
-    def __init__(self, transcript:str) -> None:
-        self._transcript = transcript
-
-    def categories(self) -> typing.Sequence[PatientCategory]:
-        return HETEROZYGOUS, HOMOZYGOUS, NO_VARIANT
-
-    def test(self, patient: Patient, query: int) -> typing.Optional[PatientCategory]:
-        if not isinstance(patient, Patient):
-            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
-        if not isinstance(query, int):
-            raise ValueError(f"query must be type integer but was type {type(query)}")
-        vars = set()
-        for var in patient.variants:
-            for trans in var.tx_annotations:
-                if trans.transcript_id == self._transcript:
-                    if trans.overlapping_exons is not None:
-                        if query in trans.overlapping_exons:
-                            vars.add(var)
-        if len(vars) == 1:
-            for v in vars:
-                gt = v.genotype_for_sample(patient.patient_id)
-                if gt == Genotype.HETEROZYGOUS:
-                    return HETEROZYGOUS
-                elif gt == Genotype.HOMOZYGOUS_ALTERNATE:
-                    return HOMOZYGOUS
-                else:
-                    # TODO - is this really what we want to return here?
-                    return HETEROZYGOUS
-        elif len(vars) > 1:
-            return HOMOZYGOUS
-        else:
-            return NO_VARIANT
-
-class ProtFeatureTypePredicate(PolyPredicate[FeatureType]):
-
-    def __init__(self, transcript:str) -> None:
-         self._transcript = transcript
-
-    def categories(self) -> typing.Sequence[PatientCategory]:
-        return HETEROZYGOUS, HOMOZYGOUS, NO_VARIANT
-
-    def test(self, patient: Patient, query:FeatureType) -> typing.Optional[PatientCategory]:
-        if not isinstance(patient, Patient):
-            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
-        if not isinstance(query, FeatureType):
-            raise ValueError(f"query must be type FeatureType but was type {type(query)}")
-        vars = set()
-        for var in patient.variants:
-            for trans in var.tx_annotations:
-                if trans.transcript_id == self._transcript:
-                    protein_location = trans.protein_effect_location
-                    if protein_location is not None:
-                        for prot in trans.protein_affected:
-                            for feat in prot.protein_features:
-                                if feat.feature_type == query:
-                                    if len(list(range(max(protein_location.start, feat.info.start), min(protein_location.end, feat.info.end) + 1))) > 0:
-                                        vars.add(var)
-        if len(vars) == 1:
-            for v in vars:
-                gt = v.genotype_for_sample(patient.patient_id)
-                if gt == Genotype.HETEROZYGOUS:
-                    return HETEROZYGOUS
-                elif gt == Genotype.HOMOZYGOUS_ALTERNATE:
-                    return HOMOZYGOUS
-                else:
-                    # TODO - is this really what we want to return here?
-                    return HETEROZYGOUS
-        elif len(vars) > 1:
-            return HOMOZYGOUS
-        else:
-            return NO_VARIANT
-
-class ProtFeaturePredicate(PolyPredicate[str]):
-    """
-    Test if the patient has a variant that overlaps with a protein feature.
-
-    The predicate is generic over a string that represents the name of the protein feature.
-    For instance, `EGF-like 2` for `FBN1 <https://www.uniprot.org/uniprotkb/P35555/entry#family_and_domains>`_
-    """
-
-    def __init__(self, transcript:str) -> None:
-         self._transcript = transcript
-
-    def categories(self) -> typing.Sequence[PatientCategory]:
-        return HETEROZYGOUS, HOMOZYGOUS, NO_VARIANT
-
-    def test(self, patient: Patient, query: str) -> typing.Optional[PatientCategory]:
-        if not isinstance(patient, Patient):
-            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
-        if not isinstance(query, str):
-            raise ValueError(f"query must be type string but was type {type(query)}")
-        vars = set()
-        for var in patient.variants:
-            for trans in var.tx_annotations:
-                if trans.transcript_id == self._transcript:
-                    protein_location = trans.protein_effect_location
-                    if protein_location is not None:
-                        for prot in trans.protein_affected:
-                            for feat in prot.protein_features:
-                                if feat.info.name == query:
-                                    if len(list(range(max(protein_location.start, feat.info.start), min(protein_location.end, feat.info.end) + 1))) > 0:
-                                        vars.add(var)
-        if len(vars) == 1:
-            for v in vars:
-                gt = v.genotype_for_sample(patient.patient_id)
-                if gt == Genotype.HETEROZYGOUS:
-                    return HETEROZYGOUS
-                elif gt == Genotype.HOMOZYGOUS_ALTERNATE:
-                    return HOMOZYGOUS
-                else:
-                    # TODO - is this really what we want to return here?
-                    return HETEROZYGOUS
-        elif len(vars) > 1:
-            return HOMOZYGOUS
-        else:
-            return NO_VARIANT
diff --git a/src/genophenocorr/analysis/predicate/_api.py b/src/genophenocorr/analysis/predicate/_api.py
index 75bacf3d..a63b4e97 100644
--- a/src/genophenocorr/analysis/predicate/_api.py
+++ b/src/genophenocorr/analysis/predicate/_api.py
@@ -5,22 +5,14 @@
 
 from genophenocorr.model import Patient
 
-T = typing.TypeVar('T')
 
-class SimplePredicate(metaclass=abc.ABCMeta):
-    """
-    Predicate for testing if a :class:`Patient` belongs to a group.
+class PatientCategory:
     """
-    def __init__(self) -> None:
-        pass
-        
-
-    @abc.abstractmethod
-    def test(self, patient: Patient) -> bool:
-        pass
+    `PatientCategory` represents one of several exclusive discrete groups.
 
-
-class PatientCategory:
+    Patient category has :attr:`cat_id`, a unique numeric identifier of the group,
+    :attr:`name` with human-readable group name, and :attr:`description` with an optional verbose description.
+    """
 
     def __init__(self, cat_id: int,
                  name: str,
@@ -31,23 +23,32 @@ def __init__(self, cat_id: int,
 
     @property
     def cat_id(self) -> int:
+        """
+        Get an `int` with the unique numeric identifier of the group.
+        """
         return self._cat_id
 
     @property
     def name(self) -> str:
+        """
+        Get a `str` with a human-readable name of the group.
+        """
         return self._name
 
     @property
     def description(self) -> typing.Optional[str]:
+        """
+        Get a `str` with an optional detailed group description.
+        """
         return self._description
 
-    def __str__(self) -> str:
+    def __repr__(self) -> str:
         return f"PatientCategory(cat_id={self.cat_id}, " \
                f"name={self.name}, " \
                f"description={self.description})"
 
-    def __repr__(self) -> str:
-        return str(self)
+    def __str__(self) -> str:
+        return self.name
 
     def __eq__(self, other) -> bool:
         return isinstance(other, PatientCategory) \
@@ -59,24 +60,78 @@ def __hash__(self) -> int:
         return hash((self.cat_id, self.name, self.description))
 
 
-class PolyPredicate(typing.Generic[T], metaclass=abc.ABCMeta):
+class PolyPredicate(metaclass=abc.ABCMeta):
     """
-    Predicate for categorizing a :class:`Patient` into one of discrete groups.
+    `PolyPredicate` bins a :class:`genophenocorr.model.Patient` into one of several discrete groups represented
+    by :class:`PatientCategory`.
+
+    The groups must be *exclusive* - the patient can be binned into one and only one group,
+    and *exhaustive* - the groups must cover all possible scenarios.
+
+    However, if the patient cannot be assigned into any meaningful category, `None` can be returned.
+
+    Note, that `PolyPredicate` must be used on a happy path - testing a patient must inherently make sense.
+    Predicate will *not* check if, for instance, the patient variants are compatible with a certain mode of inheritance.
     """
 
-    @property
+    @staticmethod
     @abc.abstractmethod
-    def categories(self) -> typing.Sequence[PatientCategory]:
+    def get_categories() -> typing.Sequence[PatientCategory]:
         """
         Get a sequence of all categories which the `PolyPredicate` can produce.
         """
         pass
 
     @abc.abstractmethod
-    def test(self, patient: Patient, query: T) -> typing.Optional[PatientCategory]:
+    def get_question(self) -> str:
+        """
+        Prepare a `str` with the question the predicate can answer.
         """
-        Assign patient into a category.
+        pass
+
+    @abc.abstractmethod
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        """
+        Assign a `patient` into a category.
 
         Return `None` if the patient cannot be assigned into any meaningful category.
         """
         pass
+
+    @classmethod
+    def n_categories(cls) -> int:
+        """
+        Get the number of categories the predicate can produce.
+        """
+        return len(cls.get_categories())
+
+    @staticmethod
+    def _check_patient(patient: Patient):
+        """
+        Check if the `patient` meets the predicate requirements.
+        """
+        if not isinstance(patient, Patient):
+            raise ValueError(f"patient must be type Patient but was type {type(patient)}")
+
+
+class BooleanPredicate(PolyPredicate, metaclass=abc.ABCMeta):
+    """
+    `BooleanPredicate` tests if a :class:`genophenocorr.model.Patient` belongs to a group and returns a boolean binning.
+    """
+
+    NO = PatientCategory(0, 'No', 'The patient does not belong to the group.')
+    """
+    Category for a patient who does *not* belong to the tested group.
+    """
+
+    YES = PatientCategory(1, 'Yes', 'The patient belongs to the group.')
+    """
+    Category for a patient who *belongs* to the tested group.
+    """
+
+    @staticmethod
+    def get_categories() -> typing.Sequence[PatientCategory]:
+        """
+        The predicate bins a patient into :class:`BooleanPredicate.NO` or :class:`BooleanPredicate.YES` category.
+        """
+        return BooleanPredicate.NO, BooleanPredicate.YES
diff --git a/src/genophenocorr/analysis/predicate/genotype/__init__.py b/src/genophenocorr/analysis/predicate/genotype/__init__.py
new file mode 100644
index 00000000..ad760f85
--- /dev/null
+++ b/src/genophenocorr/analysis/predicate/genotype/__init__.py
@@ -0,0 +1,9 @@
+from ._geno import ProtFeaturePredicate, ProtFeatureTypePredicate
+from ._geno import VariantEffectPredicate, VariantPredicate, ExonPredicate
+from ._geno import HOMOZYGOUS, HETEROZYGOUS, NO_VARIANT
+
+__all__ = [
+    'VariantEffectPredicate', 'VariantPredicate', 'ExonPredicate',
+    'ProtFeaturePredicate', 'ProtFeatureTypePredicate',
+    'HOMOZYGOUS', 'HETEROZYGOUS', 'NO_VARIANT'
+]
\ No newline at end of file
diff --git a/src/genophenocorr/analysis/predicate/genotype/_geno.py b/src/genophenocorr/analysis/predicate/genotype/_geno.py
new file mode 100644
index 00000000..3be957d8
--- /dev/null
+++ b/src/genophenocorr/analysis/predicate/genotype/_geno.py
@@ -0,0 +1,239 @@
+import typing
+
+import hpotk
+
+from genophenocorr.model import Patient, FeatureType, VariantEffect
+
+from .._api import PatientCategory, BooleanPredicate
+
+
+# TODO - should we remove these three?
+HETEROZYGOUS = PatientCategory(cat_id=0,
+                               name='Heterozygous',
+                               description="""
+                               This sample has the tested attribute on one allele.
+                               """)  #: :meta hide-value:
+
+HOMOZYGOUS = PatientCategory(cat_id=1,
+                             name='Homozygous',
+                             description="""
+                             This sample has the tested attribute on both alleles.
+                             """)  #: :meta hide-value:
+
+NO_VARIANT = PatientCategory(cat_id=2,
+                             name='No Variant',
+                             description="""
+                             The sample does not have the tested attribute.
+                             """)  #: :meta hide-value:
+
+
+
+class VariantEffectPredicate(BooleanPredicate):
+    """
+    `VariantEffectPredicate` tests if the `patient` has at least one variant that is predicted to have
+    the functional `effect` on the transcript of interest.
+
+    :param transcript_id: the accession of the transcript of interest.
+    :param effect: the tested variant effect.
+    """
+
+    def __init__(self, transcript_id: str,
+                 effect: VariantEffect) -> None:
+        self._tx_id = transcript_id
+        self._effect = hpotk.util.validate_instance(effect, VariantEffect, 'effect')
+
+    def get_question(self) -> str:
+        return f'{self._effect.name} on {self._tx_id}'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.variants) == 0:
+            return None
+
+        for variant in patient.variants:
+            for ann in variant.tx_annotations:
+                if ann.transcript_id == self._tx_id:
+                    for var_eff in ann.variant_effects:
+                        if var_eff == self._effect:
+                            return BooleanPredicate.YES
+
+        return BooleanPredicate.NO
+
+    def __str__(self):
+        return repr(self)
+
+    def __repr__(self):
+        return f'VariantEffectPredicate(transcript_id={self._tx_id}, effect={self._effect})'
+
+
+class VariantPredicate(BooleanPredicate):
+    """
+    `VariantPredicate` tests if the `patient` has ar least one allele of the variant described by the `variant_key`.
+
+    .. note::
+
+      You can get the variant key by calling :class:`genophenocorr.model.VariantCoordinates.variant_key`.
+
+    :param variant_key: a `str` with the variant key.
+    """
+
+    def __init__(self, variant_key: str) -> None:
+        self._variant_key = hpotk.util.validate_instance(variant_key, str, 'variant_key')
+
+    def get_question(self) -> str:
+        return f'>=1 allele of the variant {self._variant_key}'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.variants) == 0:
+            return None
+
+        for variant in patient.variants:
+            if variant.variant_coordinates.variant_key == self._variant_key:
+                return BooleanPredicate.YES
+
+        return BooleanPredicate.NO
+
+    def __str__(self):
+        return repr(self)
+
+    def __repr__(self):
+        return f'VariantPredicate(variant_key={self._variant_key})'
+
+
+class ExonPredicate(BooleanPredicate):
+    """
+    `ExonPredicate` tests if the `patient` has a variant that affects *n*-th exon of the transcript of interest.
+
+    .. warning::
+
+      We use 1-based numbering to number the exons, not the usual 0-based numbering of the computer science.
+      Therefore, the first exon of the transcript has ``exon_number==1``, the second exon is ``2``, and so on ...
+
+    .. warning::
+
+      We do not check if the `exon_number` spans beyond the number of exons of the given `transcript_id`!
+      Therefore, ``exon_number==10,000`` will effectively return :attr:`BooleanPredicate.FALSE` for *all* patients!!! 😱
+      Well, at least the patients of the *Homo sapiens sapiens* taxon...
+
+    :param transcript_id: the accession of the transcript of interest.
+    :param exon_number: a positive `int` of the target exon.
+    """
+
+    def __init__(self, transcript_id: str,
+                 exon_number: int) -> None:
+        self._tx_id = transcript_id
+        self._exon_number = hpotk.util.validate_instance(exon_number, int, 'exon_number')
+        if self._exon_number <= 0:
+            raise ValueError(f'`exon_number` must be a positive `int` but got {self._exon_number}')
+
+    def get_question(self) -> str:
+        return f'Variant in exon {self._exon_number} on {self._tx_id}'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.variants) == 0:
+            return None
+
+        for variant in patient.variants:
+            for ann in variant.tx_annotations:
+                if ann.transcript_id == self._tx_id:
+                    if ann.overlapping_exons is not None:
+                        if self._exon_number in ann.overlapping_exons:
+                            return BooleanPredicate.YES
+
+        return BooleanPredicate.NO
+
+    def __str__(self):
+        return repr(self)
+
+    def __repr__(self):
+        return f'ExonPredicate(tx_id={self._tx_id}, exon_number={self._exon_number})'
+
+
+class ProtFeatureTypePredicate(BooleanPredicate):
+    """
+    `ExonPredicate` tests if the `patient` has a variant that affects a :class:`FeatureType`
+    in the transcript of interest.
+
+    :param transcript_id: the accession of the transcript of interest.
+    :param feature_type: an instance of the target :class:`FeatureType`.
+    """
+
+    def __init__(self, transcript_id: str, feature_type: FeatureType) -> None:
+        self._tx_id = transcript_id
+        self._feature_type = hpotk.util.validate_instance(feature_type, FeatureType, 'feature_type')
+
+    def get_question(self) -> str:
+        return f'Variant that affecting {self._feature_type.name} protein feature on {self._tx_id}'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.variants) == 0:
+            return None
+
+        for variant in patient.variants:
+            for ann in variant.tx_annotations:
+                if ann.transcript_id == self._tx_id:
+                    prot_loc = ann.protein_effect_location
+                    for prot in ann.protein_affected:
+                        for feat in prot.protein_features:
+                            if feat.feature_type == self._feature_type:
+                                if len(list(range(max(prot_loc.start, feat.info.start), min(prot_loc.end, feat.info.end) + 1))) > 0:
+                                    return BooleanPredicate.YES
+
+        return BooleanPredicate.NO
+
+    def __str__(self):
+        return repr(self)
+
+    def __repr__(self):
+        return f'ProtFeatureTypePredicate(tx_id={self._tx_id}, feature_type={self._feature_type})'
+
+
+class ProtFeaturePredicate(BooleanPredicate):
+    """
+    `ProtFeaturePredicate` tests if the `patient` has a variant that overlaps with a protein feature.
+
+    The predicate needs the name of the protein feature.
+    For instance, `EGF-like 2` for `FBN1 <https://www.uniprot.org/uniprotkb/P35555/entry#family_and_domains>`_
+
+    :param transcript_id: the accession of the transcript of interest.
+    :param protein_feature_name: a `str` with the name of the protein feature.
+    """
+
+    def __init__(self, transcript_id: str, protein_feature_name: str) -> None:
+        self._tx_id = transcript_id
+        self._pf_name = hpotk.util.validate_instance(protein_feature_name, str, 'protein_feature_name')
+
+    def get_question(self) -> str:
+        return f'Variant in {self._pf_name} of {self._tx_id}'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.variants) == 0:
+            return None
+
+        for variant in patient.variants:
+            for ann in variant.tx_annotations:
+                if ann.transcript_id == self._tx_id:
+                    prot_loc = ann.protein_effect_location
+                    for prot in ann.protein_affected:
+                        for feat in prot.protein_features:
+                            if feat.info.name == self._pf_name:
+                                if len(list(range(max(prot_loc.start, feat.info.start), min(prot_loc.end, feat.info.end) + 1))) > 0:
+                                    return BooleanPredicate.YES
+
+        return BooleanPredicate.NO
+
+
+    def __str__(self):
+        return repr(self)
+
+    def __repr__(self):
+        return f'ProtFeaturePredicate(tx_id={self._tx_id}, exon_number={self._pf_name})'
diff --git a/src/genophenocorr/analysis/predicate/phenotype/__init__.py b/src/genophenocorr/analysis/predicate/phenotype/__init__.py
new file mode 100644
index 00000000..c1d44e15
--- /dev/null
+++ b/src/genophenocorr/analysis/predicate/phenotype/__init__.py
@@ -0,0 +1,7 @@
+from ._pheno import PhenotypePredicateFactory, PropagatingPhenotypePredicate, PropagatingPhenotypeBooleanPredicate
+from ._pheno import PropagatingPhenotypeBooleanPredicateFactory
+
+__all__ = [
+    'PropagatingPhenotypePredicate', 'PropagatingPhenotypeBooleanPredicate',
+    'PhenotypePredicateFactory', 'PropagatingPhenotypeBooleanPredicateFactory'
+]
diff --git a/src/genophenocorr/analysis/predicate/phenotype/_pheno.py b/src/genophenocorr/analysis/predicate/phenotype/_pheno.py
new file mode 100644
index 00000000..94d8a437
--- /dev/null
+++ b/src/genophenocorr/analysis/predicate/phenotype/_pheno.py
@@ -0,0 +1,166 @@
+import abc
+import typing
+
+import hpotk
+
+from genophenocorr.model import Patient
+
+from .._api import PolyPredicate, BooleanPredicate, PatientCategory
+
+
+class PhenotypePolyPredicate(PolyPredicate, metaclass=abc.ABCMeta):
+    pass
+
+
+class BasePhenotypeBooleanPredicate(PhenotypePolyPredicate, BooleanPredicate, metaclass=abc.ABCMeta):
+    pass
+
+
+class PropagatingPhenotypeBooleanPredicate(BasePhenotypeBooleanPredicate):
+
+    def __init__(self, hpo: hpotk.MinimalOntology,
+                 phenotypic_feature: hpotk.TermId,
+                 missing_implies_excluded: bool = False) -> None:
+        self._hpo = hpotk.util.validate_instance(hpo, hpotk.MinimalOntology, 'hpo')
+        self._query = hpotk.util.validate_instance(phenotypic_feature, hpotk.TermId, 'phenotypic_feature')
+        self._missing_implies_excluded = hpotk.util.validate_instance(missing_implies_excluded, bool,
+                                                                      'missing_implies_excluded')
+
+    def get_question(self) -> str:
+        query_label = self._hpo.get_term(self._query).name
+        return f'Is \'{query_label}\' present in the patient?'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.phenotypes) == 0:
+            return None
+
+        for phenotype in patient.phenotypes:
+            if phenotype.is_observed:
+                if any(self._query == anc for anc in self._hpo.graph.get_ancestors(phenotype, include_source=True)):
+                    return BooleanPredicate.YES
+            else:
+                if self._missing_implies_excluded:
+                    return BooleanPredicate.NO
+                else:
+                    if any(self._query == desc for desc in
+                           self._hpo.graph.get_descendants(phenotype, include_source=True)):
+                        return BooleanPredicate.NO
+
+        return None
+
+
+class PropagatingPhenotypePredicate(PhenotypePolyPredicate):
+    """
+    `PropagatingPhenotypePredicate` tests if the `patient` is annotated with a `query` HPO term.
+
+    The predicate returns the following results:
+
+    * :attr:`PropagatingPhenotypePredicate.PRESENT` if the patient is annotated with the `query` term or its descendant
+    * :attr:`PropagatingPhenotypePredicate.EXCLUDED` presence of the `query` term or its ancestor was specifically
+      excluded in the patient
+    * :attr:`PropagatingPhenotypePredicate.NOT_MEASURED` if the patient is not annotated with the `query` and presence of `query`
+      was *not* excluded
+    """
+
+    PRESENT = PatientCategory(cat_id=0,
+                              name='Present',
+                              description="""
+                              The sample *is* annotated with the tested phenotype feature `q`.
+
+                              This is either because the sample is annotated with `q` (exact match),
+                              or because one of sample's annotations is a descendant `q` (annotation propagation).
+                              For instance, we tested for a Seizure and the sample *had* a Clonic seizure 
+                              (a descendant of Seizure).
+                              """)  #: :meta hide-value:
+
+    EXCLUDED = PatientCategory(cat_id=1,
+                               name='Excluded',
+                               description="""
+                               We are particular about the sample *not* having the tested feature `q`.
+
+                               In other words, `q` was *excluded* in the sample or the sample is annotated with an excluded ancestor of `q`.
+
+                               For instance, we tested for a Clonic seizure and the sample did *not* have any Seizure, which implies 
+                               *not* Clonic seizure.  
+                               """)  #: :meta hide-value:
+
+    NOT_MEASURED = PatientCategory(cat_id=2,
+                                   name='Not measured',
+                                   description="""
+                                   We do not know if the sample has or has not the tested feature.
+                                   """)  #: :meta hide-value:
+
+    def __init__(self, hpo: hpotk.MinimalOntology,
+                 phenotypic_feature: hpotk.TermId) -> None:
+        self._hpo = hpotk.util.validate_instance(hpo, hpotk.MinimalOntology, 'hpo')
+        self._query = hpotk.util.validate_instance(phenotypic_feature, hpotk.TermId, 'phenotypic_feature')
+
+    @staticmethod
+    def get_categories() -> typing.Sequence[PatientCategory]:
+        return PropagatingPhenotypePredicate.PRESENT, PropagatingPhenotypePredicate.EXCLUDED, PropagatingPhenotypePredicate.NOT_MEASURED
+
+    def get_question(self) -> str:
+        query_label = self._hpo.get_term(self._query).name
+        return f'Is \'{query_label}\' present in the patient?'
+
+    def test(self, patient: Patient) -> typing.Optional[PatientCategory]:
+        self._check_patient(patient)
+
+        if len(patient.phenotypes) == 0:
+            return None
+
+        for phenotype in patient.phenotypes:
+            if phenotype.is_observed:
+                if any(self._query == anc for anc in self._hpo.graph.get_ancestors(phenotype, include_source=True)):
+                    return PropagatingPhenotypePredicate.PRESENT
+            else:
+                if any(self._query == desc for desc in self._hpo.graph.get_descendants(phenotype, include_source=True)):
+                    return self.EXCLUDED
+
+        return self.NOT_MEASURED
+
+    def __str__(self):
+        return repr(self)
+
+    def __repr__(self):
+        return f'PropagatingPhenotypePredicate(query={self._query})'
+
+
+class PhenotypePredicateFactory(metaclass=abc.ABCMeta):
+    """
+    `PhenotypePredicateFactory` creates a :class:`PhenotypePolyPredicate` to assess a phenotype query.
+    """
+
+    @staticmethod
+    @abc.abstractmethod
+    def get_categories() -> typing.Sequence[PatientCategory]:
+        """
+        Get a sequence of all categories which the `PhenotypePolyPredicate` provided by this factory can produce.
+        """
+        pass
+
+    @abc.abstractmethod
+    def get_predicate(self, query: hpotk.TermId) -> PhenotypePolyPredicate:
+        """
+        Get an instance of `PhenotypePolyPredicate` for the phenotype `query`.
+
+        :param query: a term ID of the query phenotype.
+        """
+        pass
+
+
+class PropagatingPhenotypeBooleanPredicateFactory(PhenotypePredicateFactory):
+
+    def __init__(self, hpo: hpotk.MinimalOntology,
+                 missing_implies_excluded: bool = False):
+        self._hpo = hpotk.util.validate_instance(hpo, hpotk.MinimalOntology, 'hpo')
+        self._missing_implies_excluded = missing_implies_excluded
+
+    @staticmethod
+    def get_categories() -> typing.Sequence[PatientCategory]:
+        return PropagatingPhenotypeBooleanPredicate.get_categories()
+
+    def get_predicate(self, query: hpotk.TermId) -> PhenotypePolyPredicate:
+        return PropagatingPhenotypeBooleanPredicate(self._hpo, query, self._missing_implies_excluded)
diff --git a/src/genophenocorr/model/_cohort.py b/src/genophenocorr/model/_cohort.py
index 500a32b4..6c3cebb6 100644
--- a/src/genophenocorr/model/_cohort.py
+++ b/src/genophenocorr/model/_cohort.py
@@ -65,6 +65,18 @@ def proteins(self) -> typing.Sequence[ProteinMetadata]:
         """
         return self._proteins
 
+    def present_phenotypes(self) -> typing.Iterator[Phenotype]:
+        """
+        Get an iterator over *present* phenotypes of the patient.
+        """
+        return filter(lambda p: p.is_observed, self._phenotypes)
+
+    def excluded_phenotypes(self) -> typing.Iterator[Phenotype]:
+        """
+        Get an iterator over *excluded* phenotypes of the patient.
+        """
+        return filter(lambda p: p.is_excluded, self._phenotypes)
+
     def __str__(self) -> str:
         return (f"Patient("
                 f"patient_id:{self.patient_id}, "
@@ -107,7 +119,7 @@ def from_patients(members: typing.Sequence[Patient], include_patients_with_no_HP
             cohort_phenotypes.update(patient.phenotypes)
             cohort_variants.update(patient.variants)
             var_counts.update([var.variant_coordinates.variant_key for var in patient.variants])
-            pheno_count.update([pheno.identifier.value for pheno in patient.phenotypes if pheno.observed == True])
+            pheno_count.update(pheno.identifier.value for pheno in patient.present_phenotypes())
             cohort_proteins.update(patient.proteins)
             prot_counts.update([prot.protein_id for prot in patient.proteins])
         all_counts = {'patients': len(members), 'variants': var_counts, 'phenotypes': pheno_count,
@@ -115,7 +127,7 @@ def from_patients(members: typing.Sequence[Patient], include_patients_with_no_HP
         return Cohort(members, cohort_phenotypes, cohort_variants, cohort_proteins,
                       all_counts, excluded_members)  # cohort_proteins, all_counts
 
-    """This class creates a collection of patients and makes it easier to determine overlapping diseases, 
+    """This class creates a collection of patients and makes it easier to determine overlapping diseases,
     phenotypes, variants, and proteins among the patients. If a list of JSON files is given, it will
     add each file as a patient into the grouping.
 
@@ -127,10 +139,10 @@ def from_patients(members: typing.Sequence[Patient], include_patients_with_no_HP
         all_transcripts (Sequence[string]): A set of all transcript IDs referenced in all the Variant objects
         total_patient_count (integer): The total number of Patient objects
     Methods:
-        list_all_patients(): A list of all patient IDs 
-        list_all_phenotypes(top:Optional[integer]): A list of all the top phenotype IDs (or all IDs if top is None) and the count of how many patients have it. 
-        list_all_variants(top:Optional[integer]): A list of all the top variants (or all variants if top is None) and the count of how many patients have it. 
-        list_all_proteins(top:Optional[integer]): A list of all the top protein IDs (or all IDs if top is None) and the count of how many patients have it. 
+        list_all_patients(): A list of all patient IDs
+        list_all_phenotypes(top:Optional[integer]): A list of all the top phenotype IDs (or all IDs if top is None) and the count of how many patients have it.
+        list_all_variants(top:Optional[integer]): A list of all the top variants (or all variants if top is None) and the count of how many patients have it.
+        list_all_proteins(top:Optional[integer]): A list of all the top protein IDs (or all IDs if top is None) and the count of how many patients have it.
         list_data_by_tx(transcript:Optional[string]): A list and count of all the variants effects found for all transcripts or a given transcript if transcript is not None.
     """
 
@@ -175,7 +187,7 @@ def all_phenotypes(self):
         return self._phenotype_set
 
     @property
-    def all_variants(self):
+    def all_variants(self) -> typing.Collection[Variant]:
         """
         Returns:
             set: A set of all the Variant objects in the Cohort
@@ -276,25 +288,30 @@ def get_excluded_ids(self):
     def get_excluded_count(self):
         return len(self.all_excluded_patients)
 
-    def get_protein_features_affected(self, transcript):
-        all_features = Counter()
+    def get_protein_features_affected(self, tx_id: str):
+        # TODO - add documentation, type annotations.
         protein_set = set()
-        var_coords = []
+        protein_locations = []
         for var in self.all_variants:
             for tx in var.tx_annotations:
-                if tx.transcript_id == transcript:
+                if tx.transcript_id == tx_id:
+                    # TODO - Here we are adding a sequence of `ProteinMetadata` into the `protein_set`.
+                    #  However, it looks fishy. Shouldn't we be really adding the individual `ProteinMetadata`
+                    #  into the set?
                     protein_set.add(tx.protein_affected)
-                    if tx.protein_effect_location is None or tx.protein_effect_location[0] is None or tx.protein_effect_location[1] is None:
-                        continue
-                    else:
-                        var_coords.append(tx.protein_effect_location)
+                    if tx.protein_effect_location is not None:
+                        protein_locations.append(tx.protein_effect_location)
+
         if len(protein_set) != 1:
             raise ValueError(f"Found more than 1 protein: {protein_set}")
         else:
             protein = list(protein_set)[0][0]
-        for pair in var_coords:
-            all_features.update(list(protein.get_features_variant_overlaps(pair[0], pair[1])))
+
+        all_features = Counter()
+        for location in protein_locations:
+            all_features.update(protein.get_features_variant_overlaps(location))
+
         return all_features
-        
+
     def __len__(self) -> int:
         return len(self._patient_set)
diff --git a/src/genophenocorr/model/_phenotype.py b/src/genophenocorr/model/_phenotype.py
index 7439caae..67891ec5 100644
--- a/src/genophenocorr/model/_phenotype.py
+++ b/src/genophenocorr/model/_phenotype.py
@@ -1,36 +1,28 @@
 import typing
+import warnings
 
 import hpotk
 
+
 class Phenotype(hpotk.model.Identified, hpotk.model.Named):
     """A class that represents an HPO verified phenotype
 
     Attributes:
-        identifier (hpotk.model.Named): The HPO ID associated with this phenotype
+        term_id (hpotk.model.Named): The HPO ID associated with this phenotype
         name (str): The official HPO name for this phenotype
-        observed (bool): Is True if this phenotype was observed in the respective patient
+        is_observed (bool): Is True if this phenotype was observed in the respective patient
     """
 
     @staticmethod
-    def from_term(term: hpotk.model.MinimalTerm, observed: typing.Optional[bool]):
-        return Phenotype(term.identifier, term.name, observed)
+    def from_term(term: hpotk.model.MinimalTerm, is_observed: bool):
+        return Phenotype(term.identifier, term.name, is_observed)
 
     def __init__(self, term_id: hpotk.TermId,
                  name: str,
-                 observed: typing.Optional[bool]) -> None:
-        """Constructs all necessary attributes for a Phenotype object
-
-        Args:
-            term (hpotk.model.MinimalTerm): An hpotk object that has an HPO ID and Name
-            observed (bool): Is True if this phenotype was observed in the respective patient
-        """
+                 is_observed: bool) -> None:
         self._term_id = hpotk.util.validate_instance(term_id, hpotk.TermId, 'term_id')
         self._name = hpotk.util.validate_instance(name, str, 'name')
-
-        if observed is not None:
-            if not isinstance(observed, bool):
-                raise ValueError(f"observed variable must be type boolean, but is type {type(observed)}")
-        self._observed = observed
+        self._observed = hpotk.util.validate_instance(is_observed, bool, 'is_observed')
 
     @property
     def identifier(self) -> hpotk.TermId:
@@ -58,22 +50,38 @@ def observed(self) -> typing.Optional[bool]:
         Returns:
             boolean: True if this phenotype was observed in the respective patient.
         """
+        warnings.warn('`observed` property was deprecated and will be removed in `v0.3.0`. '
+                      'Use `is_observed` instead', DeprecationWarning, stacklevel=2)
+        return self.is_observed
+
+    @property
+    def is_observed(self) -> bool:
+        """
+        Returns `True` if the phenotype was *present* in the respective patient.
+        """
         return self._observed
 
+    @property
+    def is_excluded(self) -> bool:
+        """
+        Returns `True` if the phenotype presence was *explicitly excluded* in the respective patient.
+        """
+        return not self.is_observed
+
     def __eq__(self, other):
         return isinstance(other, Phenotype) \
             and self.identifier == other.identifier \
             and self.name == other.name \
-            and self.observed == other.observed
+            and self.is_observed == other.is_observed
 
     def __hash__(self):
-        return hash((self.identifier, self.name, self.observed))
+        return hash((self.identifier, self.name, self.is_observed))
 
     def __str__(self):
         return f"Phenotype(" \
                f"identifier={self.identifier}, " \
                f"name={self.name}, " \
-               f"observed={self._observed})"
+               f"is_observed={self._observed})"
 
     def __repr__(self):
         return str(self)
diff --git a/src/genophenocorr/model/_variant.py b/src/genophenocorr/model/_variant.py
index d9a8b95e..4d18658a 100644
--- a/src/genophenocorr/model/_variant.py
+++ b/src/genophenocorr/model/_variant.py
@@ -239,9 +239,10 @@ def alt(self) -> str:
     @property
     def change_length(self) -> int:
         """
-        Returns:
-            integer: The change between the ref and alt alleles due to the variant presence. SNVs lead to change length of zero,
-            deletions and insertions/duplications lead to negative and positive change lengths, respectively.
+        Get the change of length between the `ref` and `alt` alleles due to the variant presence.
+
+        SNVs lead to change length of zero, deletions and insertions/duplications lead to negative
+        and positive change lengths, respectively.
         """
         return self._change_length
 
@@ -250,8 +251,11 @@ def variant_key(self) -> str:
         """
         Get a readable representation of the variant's coordinates.
 
-        For instance, `X_12345_12345_C_G` for sequence variant or `22_10001_20000_INV`
-        Note that both start and end coordinates use 1-based (included) coordinate system.
+        For instance, ``X_12345_12345_C_G`` for a sequence variant or ``22_10001_20000_INV`` for a symbolic variant.
+
+        .. note::
+
+          Both *start* and *end* coordinates use 1-based (included) coordinate system.
         """
         if self.is_structural():
             return f'{self.chrom}_{self.start + 1}_{self.end}_{self.alt[1:-1]}'
@@ -261,8 +265,7 @@ def variant_key(self) -> str:
     @property
     def variant_class(self) -> str:
         """
-        Returns:
-            string: The variant class. (e.g. `DUP`, `SNV`, `INS`, `MNV`, `INV`, ...)
+        Get a `str` with the variant class. (e.g. `DUP`, `SNV`, `INS`, `MNV`, `INV`, ...).
         """
 
         if self.is_structural():
@@ -281,12 +284,19 @@ def variant_class(self) -> str:
                     return 'MNV'
 
     def is_structural(self) -> bool:
-        """Checks if the variant coordinates use structural variant notation
-        (e.g. `chr5  101 . N <DEL> .  .  SVTYPE=DEL;END=120;SVLEN=-10`)
-        as opposed to the sequence/literal notation (`chr5  101 . NACGTACGTAC N`).
+        """Checks if the variant coordinates use structural variant notation as described by Variant Call Format
+        (`VCF <https://en.wikipedia.org/wiki/Variant_Call_Format>`_).
 
-        Returns:
-            boolean: True if the variant coordinates use structural variant notation
+        Ane example of *structural* variant notation::
+
+          chr5  101 . N <DEL> .  .  SVTYPE=DEL;END=120;SVLEN=-10
+
+
+        as opposed to the *sequence* (literal) notation::
+
+          chr5  101 . NACGTACGTAC N
+
+        :return: `True` if the variant coordinates use structural variant notation.
         """
         return len(self._alt) != 0 and self._alt.startswith('<') and self._alt.endswith('>')
 


Test patch:
diff --git a/src/genophenocorr/analysis/_test_fisherExact.py b/src/genophenocorr/analysis/_test_fisherExact.py
index 6a059608..e6035e12 100644
--- a/src/genophenocorr/analysis/_test_fisherExact.py
+++ b/src/genophenocorr/analysis/_test_fisherExact.py
@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
 from contextlib import nullcontext as does_not_raise
-from ._analyzers import PythonMultiFisherExact
+from ._stats import PythonMultiFisherExact
 
 
 @pytest.fixture
diff --git a/src/genophenocorr/preprocessing/_test_variant.py b/src/genophenocorr/preprocessing/_test_variant.py
index 253dc1f1..c790daf2 100644
--- a/src/genophenocorr/preprocessing/_test_variant.py
+++ b/src/genophenocorr/preprocessing/_test_variant.py
@@ -34,6 +34,7 @@ def pp_vc_finder(build: GenomeBuild,
     return PhenopacketVariantCoordinateFinder(build, hgvs_vc_finder)
 
 
+@pytest.mark.online
 @pytest.mark.parametrize("pp_path, expected",
                          [('test_data/deletion_test.json', '16_89284129_89284134_CTTTTT_C'),
                           ('test_data/insertion_test.json', '16_89280829_89280829_C_CA'),
diff --git a/tests/test_analysis.py b/tests/test_analysis.py
new file mode 100644
index 00000000..37c239d7
--- /dev/null
+++ b/tests/test_analysis.py
@@ -0,0 +1,28 @@
+import hpotk
+import pytest
+
+import pandas as pd
+
+from genophenocorr.analysis import configure_cohort_analysis
+from genophenocorr.analysis.predicate import BooleanPredicate
+from genophenocorr.data import get_toy_cohort
+from genophenocorr.model import Cohort, VariantEffect
+
+
+@pytest.mark.skip(reason='Disabled unless explicitly enabled')
+class TestCommunistCohortAnalysis:
+
+    @pytest.fixture
+    def toy_cohort(self) -> Cohort:
+        return get_toy_cohort()
+
+    @pytest.fixture
+    def hpo(self) -> hpotk.MinimalOntology:
+        return hpotk.load_minimal_ontology('/home/ielis/data/ontologies/hpo/2023-07-21/hp.2023-07-21.json')
+
+    def test_compare_by_variant_effect(self, toy_cohort: Cohort, hpo: hpotk.MinimalOntology):
+        pd.set_option('expand_frame_repr', False)
+        cohort_analysis = configure_cohort_analysis(toy_cohort, hpo)
+        results = cohort_analysis.compare_by_variant_effect(VariantEffect.MISSENSE_VARIANT, 'NM_1234.5')
+        summary = results.summarize(hpo, BooleanPredicate.YES)
+        print(summary)
diff --git a/tests/test_predicates.py b/tests/test_predicates.py
index 0db1680c..7e3fb12a 100644
--- a/tests/test_predicates.py
+++ b/tests/test_predicates.py
@@ -3,8 +3,11 @@
 import hpotk
 import pytest
 
-from genophenocorr.analysis.predicate import *
+from genophenocorr.analysis.predicate import BooleanPredicate, PatientCategory
+from genophenocorr.analysis.predicate.genotype import *
+from genophenocorr.analysis.predicate.phenotype import PropagatingPhenotypePredicate
 from genophenocorr.model import Cohort, Patient, FeatureType, VariantEffect
+
 from .fixtures import toy_hpo, toy_cohort
 
 
@@ -23,18 +26,18 @@ def find_patient(pat_id, cohort) -> typing.Optional[Patient]:
                              # Test exact match
                              ('HP:0001166',  # Arachnodactyly
                              'HetSingleVar',
-                              HPOPresentPredicate.OBSERVED),
+                              PropagatingPhenotypePredicate.PRESENT),
                              # Test inferred annotations
                              ('HP:0001250',  # Seizure
                               'HetSingleVar',
-                              HPOPresentPredicate.OBSERVED),
+                              PropagatingPhenotypePredicate.PRESENT),
                              # Test excluded feature
                              ('HP:0001257',  # Spasticity
                               'HetSingleVar',
-                              HPOPresentPredicate.NOT_OBSERVED),
+                              PropagatingPhenotypePredicate.EXCLUDED),
                             ('HP:0006280',  # Chronic pancreatitis
                               'HetSingleVar',
-                              HPOPresentPredicate.NOT_MEASURED),
+                             PropagatingPhenotypePredicate.NOT_MEASURED),
                          ])
 def test_HPOPresentPredicate(toy_cohort: Cohort,
                              toy_hpo: hpotk.Ontology,
@@ -42,88 +45,94 @@ def test_HPOPresentPredicate(toy_cohort: Cohort,
                              patient_id: str,
                              expected: PatientCategory):
     patient = find_patient(patient_id, toy_cohort)
-    predicate = HPOPresentPredicate(hpo=toy_hpo)
-    actual = predicate.test(patient, query=hpotk.TermId.from_curie(query))
+    predicate = PropagatingPhenotypePredicate(hpo=toy_hpo, phenotypic_feature=hpotk.TermId.from_curie(query))
+    actual = predicate.test(patient)
     assert actual == expected
 
 
-@pytest.mark.parametrize('patient_id, variantEffect, expected_result',
-                        (['HetSingleVar', VariantEffect.FRAMESHIFT_VARIANT, HETEROZYGOUS],
-                        ['HetSingleVar', VariantEffect.MISSENSE_VARIANT, NO_VARIANT],
-                        ['HetDoubleVar1', VariantEffect.STOP_GAINED, HETEROZYGOUS],
-                        ['HomoVar', VariantEffect.FRAMESHIFT_VARIANT, HOMOZYGOUS],
-                        ['LargeCNV', VariantEffect.STOP_LOST, HETEROZYGOUS],
-                        ['LargeCNV', VariantEffect.FEATURE_TRUNCATION, HETEROZYGOUS]))
-def test_VariantEffectPredicate(patient_id: str, 
-                                variantEffect: VariantEffect,
+@pytest.mark.parametrize('patient_id, variant_effect, expected_result',
+                         (['HetSingleVar', VariantEffect.FRAMESHIFT_VARIANT, BooleanPredicate.YES],
+                        ['HetSingleVar', VariantEffect.MISSENSE_VARIANT, BooleanPredicate.NO],
+                        ['HetDoubleVar1', VariantEffect.STOP_GAINED, BooleanPredicate.YES],
+                        ['HomoVar', VariantEffect.FRAMESHIFT_VARIANT, BooleanPredicate.YES],
+                        ['LargeCNV', VariantEffect.STOP_LOST, BooleanPredicate.YES],
+                        ['LargeCNV', VariantEffect.FEATURE_TRUNCATION, BooleanPredicate.YES]))
+def test_VariantEffectPredicate(patient_id: str,
+                                variant_effect: VariantEffect,
                                 expected_result: PatientCategory,
                                 toy_cohort: Cohort):
     patient = find_patient(patient_id, toy_cohort)
-    predicate = VariantEffectPredicate('NM_013275.6')
-    result = predicate.test(patient, variantEffect)
+    predicate = VariantEffectPredicate('NM_013275.6', effect=variant_effect)
+    result = predicate.test(patient)
     assert result == expected_result
 
 
 @pytest.mark.parametrize('patient_id, variant, hasVarResult',
-                        (['HetSingleVar', '16_89279850_89279850_G_GC', HETEROZYGOUS],
+                         (['HetSingleVar', '16_89279850_89279850_G_GC', BooleanPredicate.YES],
                         # the `HetSingleVar` patient does NOT have the variant.
-                        ['HetSingleVar', '16_89279708_89279725_AGTGTTCGGGGCGGGGCC_A', NO_VARIANT],
+                        ['HetSingleVar', '16_89279708_89279725_AGTGTTCGGGGCGGGGCC_A', BooleanPredicate.NO],
                         # but `HetDoubleVar2` below DOES have the variant.
-                        ['HetDoubleVar2', '16_89279708_89279725_AGTGTTCGGGGCGGGGCC_A', HETEROZYGOUS],
-                        ['HetDoubleVar1', '16_89284601_89284602_GG_A', HETEROZYGOUS],
-                        ['HetDoubleVar1', '16_89280752_89280752_G_T', HETEROZYGOUS],
+                        ['HetDoubleVar2', '16_89279708_89279725_AGTGTTCGGGGCGGGGCC_A', BooleanPredicate.YES],
+                        ['HetDoubleVar1', '16_89284601_89284602_GG_A', BooleanPredicate.YES],
+                        ['HetDoubleVar1', '16_89280752_89280752_G_T', BooleanPredicate.YES],
                         # the `HomoVar` patient does NOT have the variant
-                        ['HomoVar', '16_89280752_89280752_G_T', NO_VARIANT],
-                        ['HomoVar', '16_89279458_89279459_TG_T', HOMOZYGOUS],
-                        ['LargeCNV', '16_89190071_89439815_DEL', HETEROZYGOUS]))
+                        ['HomoVar', '16_89280752_89280752_G_T', BooleanPredicate.NO],
+                        ['HomoVar', '16_89279458_89279459_TG_T', BooleanPredicate.YES],
+                        ['LargeCNV', '16_89190071_89439815_DEL', BooleanPredicate.YES]))
 def test_VariantPredicate(patient_id, variant, hasVarResult, toy_cohort):
-    predicate = VariantPredicate('NM_013275.6')
+    predicate = VariantPredicate(variant_key=variant)
     patient = find_patient(patient_id, toy_cohort)
-    result = predicate.test(patient, variant)
+    result = predicate.test(patient)
     assert result == hasVarResult
 
 
 @pytest.mark.parametrize('patient_id, exon, hasVarResult',
-                        (['HetSingleVar', 9, HETEROZYGOUS],
-                        ['HetSingleVar', 13, NO_VARIANT],
-                        ['HetDoubleVar1', 9, HOMOZYGOUS],
-                        ['HetDoubleVar2', 10, HETEROZYGOUS],
-                        ['HetDoubleVar2', 9, HETEROZYGOUS],
-                        ['HomoVar', 10, NO_VARIANT],
-                        ['HomoVar', 9, HOMOZYGOUS],
-                        ['LargeCNV', 1, NO_VARIANT],
-                        ['LargeCNV', 13, HETEROZYGOUS]))
+                         (['HetSingleVar', 9, BooleanPredicate.YES],
+                        ['HetSingleVar', 13, BooleanPredicate.NO],
+                        ['HetDoubleVar1', 9, BooleanPredicate.YES],
+                        ['HetDoubleVar2', 10, BooleanPredicate.YES],
+                        ['HetDoubleVar2', 9, BooleanPredicate.YES],
+                        ['HomoVar', 10, BooleanPredicate.NO],
+                        ['HomoVar', 9, BooleanPredicate.YES],
+                        ['LargeCNV', 1, BooleanPredicate.NO],
+                        ['LargeCNV', 13, BooleanPredicate.YES]))
 def test_ExonPredicate(patient_id, exon, hasVarResult, toy_cohort):
     patient = find_patient(patient_id, toy_cohort)
-    predicate = ExonPredicate('NM_013275.6')
-    result = predicate.test(patient, exon)
+    predicate = ExonPredicate('NM_013275.6', exon_number=exon)
+    result = predicate.test(patient)
     assert result == hasVarResult
 
 
-@pytest.mark.parametrize('patient_id, featureType, hasVarResult',
-                        (['HetDoubleVar2', FeatureType.REGION, HOMOZYGOUS],
-                        ['HetDoubleVar2', FeatureType.REPEAT, NO_VARIANT],
-                        ['HetSingleVar', FeatureType.REGION, HETEROZYGOUS],
-                        ['HomoVar', FeatureType.REGION, HOMOZYGOUS],
-                        ['HetDoubleVar1', FeatureType.REPEAT, NO_VARIANT]))
+@pytest.mark.parametrize('patient_id, feature_type, hasVarResult',
+                         (['HetDoubleVar2', FeatureType.REGION, BooleanPredicate.YES],
+                        ['HetDoubleVar2', FeatureType.REPEAT, BooleanPredicate.NO],
+                        ['HetSingleVar', FeatureType.REGION, BooleanPredicate.YES],
+                        ['HomoVar', FeatureType.REGION, BooleanPredicate.YES],
+                        ['HetDoubleVar1', FeatureType.REPEAT, BooleanPredicate.NO]))
                         ## TODO Why do CNV not show as affecting a feature?
                         ##['LargeCNV', FeatureType.REGION , HETEROZYGOUS]))
-def test_ProteinFeatureTypePredicate(patient_id, featureType, hasVarResult, toy_cohort):
+def test_ProteinFeatureTypePredicate(patient_id, feature_type, hasVarResult, toy_cohort):
     patient = find_patient(patient_id, toy_cohort)
-    predicate = ProtFeatureTypePredicate('NM_013275.6')
-    result = predicate.test(patient, featureType)
+    predicate = ProtFeatureTypePredicate('NM_013275.6', feature_type=feature_type)
+    result = predicate.test(patient)
     assert result == hasVarResult
 
 
 @pytest.mark.parametrize('patient_id, feature, hasVarResult',
-                        (['HetDoubleVar2', 'Disordered', HETEROZYGOUS],
-                        ['HetDoubleVar2', 'BadFeature', NO_VARIANT],
-                        ['HetSingleVar', 'Disordered', HETEROZYGOUS],
-                        ['HomoVar', 'Disordered', HOMOZYGOUS],
-                        ['HetDoubleVar1', 'Disordered', HETEROZYGOUS]))
+                         (['HetDoubleVar2', 'Disordered', BooleanPredicate.YES],
+                        ['HetDoubleVar2', 'BadFeature', BooleanPredicate.NO],
+                        ['HetSingleVar', 'Disordered', BooleanPredicate.YES],
+                        ['HomoVar', 'Disordered', BooleanPredicate.YES],
+                        ['HetDoubleVar1', 'Disordered', BooleanPredicate.YES]))
 def test_ProteinFeaturePredicate(patient_id, feature, hasVarResult, toy_cohort):
-    predicate = ProtFeaturePredicate('NM_013275.6')
+    predicate = ProtFeaturePredicate('NM_013275.6', protein_feature_name=feature)
     patient = find_patient(patient_id, toy_cohort)
-    result = predicate.test(patient, feature)
+    result = predicate.test(patient)
     assert result == hasVarResult
 
+
+class TestGenericPredicateMethods:
+
+    def test_get_categories(self):
+        assert len(PropagatingPhenotypePredicate.get_categories()) == PropagatingPhenotypePredicate.n_categories()
+        assert len(BooleanPredicate.get_categories()) == BooleanPredicate.n_categories()