--- a
+++ b/openomics/database/ontology.py
@@ -0,0 +1,917 @@
+import os
+import warnings
+from collections.abc import Iterable
+from io import TextIOWrapper, StringIO
+from typing import Tuple, List, Dict, Union, Callable, Optional
+
+import dask.dataframe as dd
+import networkx as nx
+import numpy as np
+import obonet
+import pandas as pd
+import scipy.sparse as ssp
+from logzero import logger
+from networkx import NetworkXError
+from pandas import DataFrame
+
+from openomics.io.read_gaf import read_gaf
+from openomics.transforms.adj import slice_adj
+from openomics.transforms.agg import get_agg_func
+from .base import Database
+
+__all__ = ['GeneOntology', 'UniProtGOA', 'InterPro', 'HumanPhenotypeOntology', ]
+
+
+class Ontology(Database):
+    annotations: pd.DataFrame
+
+    def __init__(self,
+                 path,
+                 file_resources=None,
+                 blocksize=None,
+                 **kwargs):
+        """
+        Manages dataset input processing from tables and construct an ontology network from .obo file. There ontology
+        network is G(V,E) where there exists e_ij for child i to parent j to present "node i is_a node j".
+
+        Args:
+            path:
+            file_resources:
+            blocksize:
+            verbose:
+        """
+        super().__init__(path=path, file_resources=file_resources, blocksize=blocksize, **kwargs)
+
+        self.network, self.node_list = self.load_network(self.file_resources)
+        self.annotations = self.load_annotation(self.file_resources, self.blocksize)
+
+        self.close()
+
+    def load_network(self, file_resources) -> Tuple[nx.MultiDiGraph, List[str]]:
+        raise NotImplementedError()
+
+    def load_annotation(self, file_resources, blocksize) -> Union[pd.DataFrame, dd.DataFrame]:
+        pass
+
+    def filter_network(self, namespace) -> None:
+        """
+        Filter the subgraph node_list to only `namespace` terms.
+        Args:
+            namespace: one of {"biological_process", "cellular_component", "molecular_function"}
+        """
+        terms = self.data[self.data["namespace"] == namespace]["go_id"].unique()
+        print("{} terms: {}".format(namespace,
+                                    len(terms))) if self.verbose else None
+        self.network = self.network.subgraph(nodes=list(terms))
+        self.node_list = np.array(list(terms))
+
+    def adj(self, node_list):
+        adj_mtx = nx.adj_matrix(self.network, nodelist=node_list)
+
+        if node_list is None or list(node_list) == list(self.node_list):
+            return adj_mtx
+        elif set(node_list) < set(self.node_list):
+            return slice_adj(adj_mtx, list(self.node_list), node_list,
+                             None)
+        elif not (set(node_list) < set(self.node_list)):
+            raise Exception("A node in node_list is not in self.node_list.")
+
+        return adj_mtx
+
+    def filter_annotation(self, annotation: pd.Series):
+        go_terms = set(self.node_list)
+        filtered_annotation = annotation.map(lambda x: list(set(x) & go_terms)
+                                             if isinstance(x, list) else [])
+
+        return filtered_annotation
+
+    def get_child_nodes(self):
+        adj = self.adj(self.node_list)
+        leaf_terms = self.node_list[np.nonzero(adj.sum(axis=0) == 0)[1]]
+        return leaf_terms
+
+    def get_root_nodes(self):
+        adj = self.adj(self.node_list)
+        parent_terms = self.node_list[np.nonzero(adj.sum(axis=1) == 0)[0]]
+        return parent_terms
+
+    def get_dfs_paths(self, root_nodes: list, filter_duplicates=False):
+        """
+        Return all depth-first search paths from root node(s) to children node by traversing the ontology directed graph.
+        Args:
+            root_nodes (list): ["GO:0008150"] if biological processes, ["GO:0003674"] if molecular_function, or ["GO:0005575"] if cellular_component
+            filter_duplicates (bool): whether to remove duplicated paths that end up at the same leaf nodes
+
+        Returns: pd.DataFrame of all paths starting from the root nodes.
+        """
+        if not isinstance(root_nodes, list):
+            root_nodes = list(root_nodes)
+
+        paths = list(dfs_path(self.network, root_nodes))
+        paths = list(flatten_list(paths))
+        paths_df = pd.DataFrame(paths)
+
+        if filter_duplicates:
+            paths_df = paths_df[~paths_df.duplicated(keep="first")]
+            paths_df = filter_dfs_paths(paths_df)
+
+        return paths_df
+
+    def remove_predecessor_terms(self, annotation: pd.Series, sep="\||;"):
+        # leaf_terms = self.get_child_nodes()
+        # if not annotation.map(lambda x: isinstance(x, (list, np.ndarray))).any() and sep:
+        #     annotation = annotation.str.split(sep)
+        #
+        # parent_terms = annotation.map(lambda x: list(
+        #     set(x) & set(leaf_terms)) if isinstance(x, (list, np.ndarray)) else None)
+        # return parent_terms
+        raise NotImplementedError
+
+    def get_subgraph(self, edge_types: Union[str, List[str]]) -> Union[nx.MultiDiGraph, nx.DiGraph]:
+        if not hasattr(self, "_subgraphs"):
+            self._subgraphs = {}
+        elif edge_types in self._subgraphs:
+            return self._subgraphs[edge_types]
+
+        if edge_types and isinstance(self.network, (nx.MultiGraph, nx.MultiDiGraph)):
+            # Needed to create new nx.Graph because .edge_subgraph is too slow to iterate on (idk why)
+            g = nx.from_edgelist([(u, v) for u, v, k in self.network.edges if k in edge_types],
+                                 create_using=nx.DiGraph if self.network.is_directed() else nx.Graph)
+        else:
+            raise Exception("Must provide `edge_types` keys for a nx.MultiGraph type.")
+
+        self._subgraphs[edge_types] = g
+
+        return g
+
+    def add_predecessor_terms(self, anns: pd.Series, edge_type: Union[str, List[str]] = 'is_a', sep="\||;"):
+        anns_w_parents = anns.map(lambda x: [] if not isinstance(x, (list, np.ndarray)) else x) + \
+                         get_predecessor_terms(anns, self.get_subgraph(edge_type))
+
+        return anns_w_parents
+
+    @staticmethod
+    def get_node_color(file="~/Bioinformatics_ExternalData/GeneOntology/go_colors_biological.csv", ):
+        go_colors = pd.read_csv(file)
+
+        def selectgo(x):
+            terms = [term for term in x if isinstance(term, str)]
+            if len(terms) > 0:
+                return terms[-1]
+            else:
+                return None
+
+        go_colors["node"] = go_colors[[
+            col for col in go_colors.columns if col.isdigit()
+        ]].apply(selectgo, axis=1)
+        go_id_colors = go_colors[go_colors["node"].notnull()].set_index("node")["HCL.color"]
+        go_id_colors = go_id_colors[~go_id_colors.index.duplicated(keep="first")]
+
+        print(go_id_colors.unique().shape, go_colors["HCL.color"].unique().shape)
+        return go_id_colors
+
+    def split_annotations(self, src_node_col="gene_name", dst_node_col="go_id", groupby: List[str] = ["Qualifier"],
+                          train_date="2017-06-15", valid_date="2017-11-15", test_date="2021-12-31",
+                          query: Optional[str] = "Evidence in ['EXP', 'IDA', 'IPI', 'IMP', 'IGI', 'IEP', 'TAS', 'IC']",
+                          filter_src_nodes: pd.Index = None, filter_dst_nodes: pd.Index = None,
+                          agg: Union[Callable, str] = "unique") -> Tuple[DataFrame, DataFrame, DataFrame]:
+        """
+
+        Args:
+            src_node_col (str): Name of column containg the the src node types.
+            dst_node_col (str): Name of column containg the the dst node types.
+            train_date (str): A date before which the annotations belongs in the training set.
+            valid_date (str): A date before which the annotations belongs in the validation set.
+            test_date (str): A date before which the annotations belongs in the testing set.
+            groupby (str): A list of strings to groupby annotations on, default [`src_node_col`, "Qualifier"].
+            query (str, optional): A pandas query string to filter annotations. Default, only select ['EXP', 'IDA', 'IPI', 'IMP', 'IGI', 'IEP', 'TAS', 'IC'] annotations.
+            filter_src_nodes (pd.Index): A subset annotations by these values on `src_node_col`.
+            filter_dst_nodes (pd.Index): A subset annotations by these values on `dst_node_col`.
+            agg (str): Either "unique" or "add_parent", or a callable function, or a dd.Aggregation() for aggregating on the `dst_node_col` column after groupby on `groupby`.
+        """
+        raise NotImplementedError
+
+
+class GeneOntology(Ontology):
+    """Loads the GeneOntology database from http://geneontology.org .
+
+    Default path: "http://geneontology.org/gene-associations/".
+
+    Default file_resources: {
+        "go-basic.obo": "http://purl.obolibrary.org/obo/go/go-basic.obo",
+        "goa_human.gaf": "goa_human.gaf.gz",
+        "goa_human_rna.gaf": "goa_human_rna.gaf.gz",
+        "goa_human_isoform.gaf": "goa_human_isoform.gaf.gz",
+    }
+    """
+    COLUMNS_RENAME_DICT = {
+        "DB_Object_Symbol": "gene_name",
+        "DB_Object_ID": "gene_id",
+        "GO_ID": "go_id",
+        "Taxon_ID": 'species_id',
+    }
+
+    DROP_COLS = {'DB:Reference', 'With', 'Annotation_Extension', 'Gene_Product_Form_ID'}
+
+    def __init__(
+        self,
+        path="http://geneontology.org/gene-associations/",
+        species='human',
+        file_resources=None,
+        index_col='DB_Object_Symbol',
+        keys=None,
+        col_rename=COLUMNS_RENAME_DICT,
+        blocksize=None,
+        **kwargs
+    ):
+        """
+        Loads the GeneOntology database from http://geneontology.org .
+
+            Default path: "http://geneontology.org/gene-associations/" .
+            Default file_resources: {
+                "go-basic.obo": "http://purl.obolibrary.org/obo/go/go-basic.obo",
+                "goa_human.gaf": "goa_human.gaf.gz",
+                "goa_human_rna.gaf": "goa_human_rna.gaf.gz",
+                "goa_human_isoform.gaf": "goa_human_isoform.gaf.gz",
+            }
+
+        Data for GO term annotations in .gpi files are already included in .obo file, so this module doesn't maker use of .gpi files.
+
+        Handles downloading the latest Gene Ontology obo and annotation data, preprocesses them. It provide
+        functionalities to create a directed acyclic graph of GO terms, filter terms, and filter annotations.
+        """
+        if species and not hasattr(self, 'species'):
+            self.species = species.lower()
+        elif species is None:
+            self.species = 'uniprot'
+
+        if file_resources is None:
+            file_resources = {
+                f"goa_{self.species}.gaf.gz": f"goa_{self.species}.gaf.gz",
+            }
+            if species != 'uniprot':
+                file_resources[f"goa_{self.species}_rna.gaf.gz"] = f"goa_{self.species}_rna.gaf.gz"
+                file_resources[f"goa_{self.species}_isoform.gaf.gz"] = f"goa_{self.species}_isoform.gaf.gz"
+
+        if not any('.obo' in file for file in file_resources):
+            warnings.warn(
+                f'No .obo file provided in `file_resources`, so automatically adding "http://purl.obolibrary.org/obo/go/go-basic.obo"')
+            file_resources["go-basic.obo"] = "http://purl.obolibrary.org/obo/go/go-basic.obo"
+
+        super().__init__(path, file_resources, index_col=index_col, keys=keys, col_rename=col_rename,
+                         blocksize=blocksize, **kwargs)
+
+        # By default, the __init__ constructor run load_dataframe() before load_network(), but for GeneOntology,
+        # we get node data from the nx.Graph, so we must run load_dataframe() again when self.network is not None.
+        self.data = self.load_dataframe(self.file_resources, self.blocksize)
+
+    def info(self):
+        print("network {}".format(nx.info(self.network)))
+
+    def load_dataframe(self, file_resources: Dict[str, TextIOWrapper], blocksize=None) -> DataFrame:
+        if hasattr(self, 'network') and self.network is not None:
+            # Annotations for each GO term from nodes in the NetworkX graph created by the .obo file
+            go_terms = pd.DataFrame.from_dict(dict(self.network.nodes(data=True)), orient='index')
+            go_terms["def"] = go_terms["def"].apply(
+                lambda x: x.split('"')[1] if isinstance(x, str) else None)
+            go_terms.index.name = "go_id"
+
+            # Get depth of each term node in its ontology
+            hierarchy = nx.subgraph_view(self.network, filter_edge=lambda u, v, e: e == 'is_a')
+            for namespace in go_terms['namespace'].unique():
+                root_term = go_terms.query(f'name == "{namespace}"').index.item()
+                namespace_terms = go_terms.query(f'namespace == "{namespace}"').index
+                shortest_paths = nx.shortest_path_length(hierarchy.subgraph(namespace_terms), root_term)
+                go_terms.loc[namespace_terms, 'depth'] = namespace_terms.map(shortest_paths)
+            go_terms['depth'] = go_terms['depth'].fillna(0).astype(int)
+        else:
+            go_terms = None
+
+        return go_terms
+
+    def load_network(self, file_resources) -> Tuple[nx.Graph, np.ndarray]:
+        network, node_list = None, None
+        fn = next((fn for fn in file_resources if fn.endswith(".obo")), None)
+        if fn:
+            network: nx.MultiDiGraph = obonet.read_obo(file_resources[fn])
+            network = network.reverse(copy=True)
+            node_list = np.array(network.nodes)
+
+        return network, node_list
+
+    def load_annotation(self, file_resources, blocksize=None) -> Union[pd.DataFrame, dd.DataFrame]:
+        # Handle .gaf annotation files
+        dfs = {}
+        for filename, filepath_or_buffer in file_resources.items():
+            gaf_name = filename.split(".")[0]
+            # Ensure no duplicate GAF file (if having files uncompressed with same prefix)
+            if gaf_name in dfs: continue
+
+            if blocksize and isinstance(filepath_or_buffer, str):
+                if filename.endswith(".processed.parquet"):
+                    # Parsed and filtered gaf file
+                    dfs[gaf_name] = dd.read_parquet(filepath_or_buffer, chunksize=blocksize)
+                    if dfs[gaf_name].index.name != self.index_col and self.index_col in dfs[gaf_name].columns:
+                        dfs[gaf_name] = dfs[gaf_name].set_index(self.index_col, sorted=True)
+                    if not dfs[gaf_name].known_divisions:
+                        dfs[gaf_name].divisions = dfs[gaf_name].compute_current_divisions()
+
+                elif (filename.endswith(".parquet") or filename.endswith(".gaf")):
+                    # .parquet from .gaf.gz file, unfiltered, with raw str values
+                    dfs[gaf_name] = read_gaf(filepath_or_buffer, blocksize=blocksize, index_col=self.index_col,
+                                             keys=self.keys, usecols=self.usecols)
+
+                elif filename.endswith(".gaf.gz"):
+                    # Compressed .gaf file downloaded
+                    dfs[gaf_name] = read_gaf(filepath_or_buffer, blocksize=blocksize, index_col=self.index_col,
+                                             keys=self.keys, usecols=self.usecols, compression='gzip')
+
+            else:
+                if filename.endswith(".processed.parquet"):
+                    dfs[gaf_name] = pd.read_parquet(filepath_or_buffer)
+                if filename.endswith(".gaf"):
+                    dfs[gaf_name] = read_gaf(filepath_or_buffer, index_col=self.index_col, keys=self.keys,
+                                             usecols=self.usecols)
+
+        if len(dfs):
+            annotations = dd.concat(list(dfs.values()), interleave_partitions=True) \
+                if blocksize else pd.concat(dfs.values())
+
+            annotations = annotations.rename(columns=UniProtGOA.COLUMNS_RENAME_DICT)
+            if annotations.index.name in UniProtGOA.COLUMNS_RENAME_DICT:
+                annotations.index = annotations.index.rename(
+                    UniProtGOA.COLUMNS_RENAME_DICT[annotations.index.name])
+        else:
+            annotations = None
+
+        return annotations
+
+    def split_annotations(self, src_node_col="gene_name", dst_node_col="go_id", groupby: List[str] = ["Qualifier"],
+                          train_date="2017-06-15", valid_date="2017-11-15", test_date="2021-12-31",
+                          query: str = "Evidence in ['EXP', 'IDA', 'IPI', 'IMP', 'IGI', 'IEP', 'TAS', 'IC']",
+                          filter_src_nodes: pd.Index = None, filter_dst_nodes: pd.Index = None,
+                          agg: Union[str, Callable, dd.Aggregation] = "unique") \
+        -> Tuple[DataFrame, DataFrame, DataFrame]:
+        assert isinstance(groupby, list) and groupby, f"`groupby` must be a nonempty list of strings. Got {groupby}"
+
+        # Set the source column (i.e. protein_id or gene_name), to be the first in groupby
+        if src_node_col not in groupby:
+            groupby = [src_node_col] + groupby
+        if "Qualifier" not in groupby and "Qualifier" in self.annotations.columns:
+            groupby.append("Qualifier")
+
+        # Aggregator function
+        if agg == "add_parent":
+            subgraph = self.get_subgraph(edge_types="is_a")
+            node_ancestors = {node: nx.ancestors(subgraph, node) for node in subgraph.nodes}
+
+            if isinstance(self.annotations, dd.DataFrame):
+                agg = dd.Aggregation(name='_unique_add_parent',
+                                     chunk=lambda s: s.unique(),
+                                     agg=lambda s0: s0.apply(get_predecessor_terms, node_ancestors, keep_terms=True),
+                                     finalize=lambda s1: s1.apply(lambda li: np.hstack(li) if li else None))
+            else:
+                agg = lambda s: get_predecessor_terms(s, g=node_ancestors, join_groups=True, keep_terms=True)
+
+        elif agg == 'unique' and isinstance(self.annotations, dd.DataFrame):
+            agg = get_agg_func('unique', use_dask=True)
+
+        elif isinstance(self.annotations, dd.DataFrame) and not isinstance(agg, dd.Aggregation):
+            raise Exception("`agg` must be a dd.Aggregation for groupby.agg() on columns of a dask DataFrame")
+
+        def _remove_dup_neg_go_id(s: pd.Series) -> pd.Series:
+            if s.isna().any():
+                return s
+            elif isinstance(s[neg_dst_col], Iterable) and isinstance(s[dst_node_col], Iterable):
+                rm_dups_go_id = [go_id for go_id in s[neg_dst_col] if go_id not in s[dst_node_col]]
+                if len(rm_dups_go_id) == 0:
+                    rm_dups_go_id = None
+                s[neg_dst_col] = rm_dups_go_id
+            return s
+
+        neg_dst_col = f"neg_{dst_node_col}"
+
+        # Filter annotations
+        annotations = self.annotations
+        if query:
+            annotations = annotations.query(query)
+        if filter_src_nodes is not None:
+            annotations = annotations[annotations[src_node_col].isin(filter_src_nodes)]
+        if filter_dst_nodes is not None:
+            annotations = annotations[annotations[dst_node_col].isin(filter_dst_nodes)]
+        if annotations.index.name in groupby:
+            annotations = annotations.reset_index()
+
+        # Split train/valid/test annotations
+        train_anns = annotations.loc[annotations["Date"] <= pd.to_datetime(train_date)]
+        valid_anns = annotations.loc[(annotations["Date"] <= pd.to_datetime(valid_date)) & \
+                                     (annotations["Date"] > pd.to_datetime(train_date))]
+        test_anns = annotations.loc[(annotations["Date"] <= pd.to_datetime(test_date)) & \
+                                    (annotations["Date"] > pd.to_datetime(valid_date))]
+
+        outputs = []
+        for anns in [train_anns, valid_anns, test_anns]:
+            # Keep track of which annotation has a "NOT" Qualifier
+            is_neg_ann = anns.loc[:, "Qualifier"].map(lambda li: "NOT" in li)
+
+            # Convert `Qualifiers` entries of list of strings to string
+            args = dict(meta=pd.Series([""])) if isinstance(anns, dd.DataFrame) else {}
+            anns.loc[:, 'Qualifier'] = anns.loc[:, 'Qualifier'].apply(
+                lambda li: "".join([i for i in li if i != "NOT"]), **args)
+
+            # Aggregate gene-GO annotations
+            if isinstance(anns, pd.DataFrame) and len(anns.index):
+                pos_anns = anns[~is_neg_ann].groupby(groupby).agg({dst_node_col: agg})
+                neg_anns = anns[is_neg_ann].groupby(groupby).agg(**{neg_dst_col: (dst_node_col, agg)})
+                pos_neg_anns = pd.concat([pos_anns, neg_anns], axis=1)
+
+            elif isinstance(anns, dd.DataFrame) and len(anns.index) and dst_node_col in anns.columns:
+                pos_anns = anns[~is_neg_ann].groupby(groupby).agg({dst_node_col: agg})
+                if False and len(is_neg_ann.index):
+                    neg_anns = anns[is_neg_ann].groupby(groupby).agg({dst_node_col: agg})
+                    neg_anns.columns = [neg_dst_col]
+                    pos_neg_anns = dd.concat([pos_anns, neg_anns], axis=1)
+                else:
+                    pos_neg_anns = pos_anns
+                    pos_neg_anns[neg_dst_col] = None
+
+            else:
+                pos_neg_anns = pd.DataFrame(
+                    columns=[dst_node_col, neg_dst_col],
+                    index=pd.MultiIndex(levels=[[] for i in range(len(groupby))],
+                                        codes=[[] for i in range(len(groupby))], names=groupby))
+                outputs.append(pos_neg_anns)
+                continue
+
+            if isinstance(pos_neg_anns, pd.DataFrame):
+                pos_neg_anns = pos_neg_anns.drop([""], axis='index', errors="ignore")
+
+            # Remove "GO:0005515" (protein binding) annotations for a gene if it's the gene's only annotation
+            _exclude_single_fn = lambda li: None \
+                if isinstance(li, Iterable) and len(li) == 1 and "GO:0005515" in li else li
+            args = dict(meta=pd.Series([list()])) if isinstance(anns, dd.DataFrame) else {}
+            pos_neg_anns.loc[:, dst_node_col] = pos_neg_anns[dst_node_col].apply(_exclude_single_fn, **args)
+
+            # Drop rows with all nan values
+            if isinstance(pos_neg_anns, pd.DataFrame):
+                pos_neg_anns = pos_neg_anns.drop(pos_neg_anns.index[pos_neg_anns.isna().all(1)], axis='index')
+
+            # Ensure no negative terms duplicates positive annotations
+            if len(is_neg_ann.index):
+                args = dict(meta=pd.DataFrame({dst_node_col: [], neg_dst_col: []})) \
+                    if isinstance(anns, dd.DataFrame) else {}
+                pos_neg_anns = pos_neg_anns.apply(_remove_dup_neg_go_id, axis=1, **args)
+
+            outputs.append(pos_neg_anns)
+
+        return tuple(outputs)
+
+
+def get_predecessor_terms(anns: Union[pd.Series, Iterable], g: Union[Dict[str, List[str]], nx.MultiDiGraph],
+                          join_groups=False, keep_terms=True, exclude={'GO:0005575', 'GO:0008150', 'GO:0003674'}) \
+    -> Union[pd.Series, List[str]]:
+    """
+
+    Args:
+        anns ():
+        g (nx.MultiDiGraph, Dict[str,Set[str]]): Either a NetworkX DAG or a precomputed lookup table of node to ancestors
+        join_groups (): whether to concatenate multiple
+        keep_terms ():
+        exclude ():
+
+    Returns:
+
+    """
+    if exclude is None:
+        exclude = {}
+
+    def _get_ancestors(terms: Iterable):
+        try:
+            if isinstance(terms, Iterable):
+                if isinstance(g, dict):
+                    parents = {parent \
+                               for term in terms if term in g \
+                               for parent in g[term] if parent not in exclude}
+
+                elif isinstance(g, nx.Graph):
+                    parents = {parent \
+                               for term in terms if term in g.nodes \
+                               for parent in nx.ancestors(g, term) if parent not in exclude}
+                else:
+                    raise Exception("Provided `g` arg must be either an nx.Graph or a Dict")
+            else:
+                parents = []
+
+            if keep_terms and isinstance(terms, list):
+                terms.extend(parents)
+                out = terms
+            else:
+                out = list(parents)
+
+        except (NetworkXError, KeyError) as nxe:
+            if "foo" not in nxe.__str__():
+                logger.error(f"{nxe.__class__.__name__} get_predecessor_terms._get_ancestors: {nxe}")
+            out = terms if keep_terms else []
+
+        return out
+
+    if isinstance(anns, pd.Series):
+        if (anns.map(type) == str).any():
+            anns = anns.map(lambda s: [s])
+
+        parent_terms = anns.map(_get_ancestors)
+        if join_groups:
+            parent_terms = sum(parent_terms, [])
+
+    elif isinstance(anns, Iterable):
+        parent_terms = _get_ancestors(anns)
+
+    else:
+        parent_terms = []
+
+    return parent_terms
+
+class UniProtGOA(GeneOntology):
+    """Loads the GeneOntology database from https://www.ebi.ac.uk/GOA/ .
+
+    Default path: "ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/UNIPROT/" .
+    Default file_resources: {
+        "goa_uniprot_all.gaf": "goa_uniprot_all.gaf.gz",
+    }
+    """
+
+    COLUMNS_RENAME_DICT = {
+        "DB_Object_ID": "protein_id",
+        "DB_Object_Symbol": "gene_name",
+        "GO_ID": "go_id",
+        "Taxon_ID": 'species_id',
+    }
+    def __init__(
+        self,
+        path="ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/",
+        species="HUMAN",
+        file_resources=None,
+        index_col='DB_Object_ID', keys=None,
+        col_rename=COLUMNS_RENAME_DICT,
+        blocksize=None,
+        **kwargs,
+    ):
+        """
+        Loads the UniProtGOA database from https://www.ebi.ac.uk/GOA/ .
+
+            Default path: "ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/UNIPROT/" .
+            Default file_resources: {
+                "goa_uniprot_all.gaf.gz": "goa_uniprot_all.gaf.gz",
+                "go.obo": "http://current.geneontology.org/ontology/go.obo",
+            }
+
+        Handles downloading the latest Gene Ontology obo and annotation data, preprocesses them. It provides
+        functionalities to create a directed acyclic graph of GO terms, filter terms, and filter annotations.
+
+        Args:
+            path ():
+            species ():
+            file_resources ():
+            index_col ():
+            keys ():
+            col_rename ():
+            blocksize ():
+            **kwargs ():
+        """
+        if species is None:
+            self.species = species = 'UNIPROT'
+            substr = 'uniprot_all'
+        else:
+            self.species = species.upper()
+            substr = species.lower()
+
+        if file_resources is None:
+            file_resources = {
+                "go.obo": "http://current.geneontology.org/ontology/go.obo",
+                f"goa_{self.species.lower()}.gaf.gz": os.path.join(species, f"goa_{substr}.gaf.gz"),
+                # f"goa_{self.species.lower()}_isoform.gaf.gz": os.path.join(species, f"goa_{substr}_isoform.gaf.gz"),
+                # f"goa_{self.species.lower()}_complex.gaf.gz": os.path.join(species, f"goa_{substr}_complex.gaf.gz"),
+            }
+
+        if not any('.obo' in file for file in file_resources):
+            warnings.warn(f'No .obo file provided in `file_resources`, '
+                          f'so automatically adding "http://purl.obolibrary.org/obo/go/go-basic.obo"')
+            file_resources["go-basic.obo"] = "http://purl.obolibrary.org/obo/go/go-basic.obo"
+
+        super().__init__(path=path, file_resources=file_resources, index_col=index_col, keys=keys,
+                         col_rename=col_rename,
+                         blocksize=blocksize, **kwargs)
+
+
+class InterPro(Ontology):
+    """
+    Default parameters
+    path="https://ftp.ebi.ac.uk/pub/databases/interpro/current_release/"
+    file_resources = {
+        "entry.list": "entry.list",
+        "protein2ipr.dat.gz": "protein2ipr.dat.gz",
+        "interpro2go": "interpro2go",
+        "ParentChildTreeFile.txt": "ParentChildTreeFile.txt",
+    }
+    """
+
+    def __init__(self, path="https://ftp.ebi.ac.uk/pub/databases/interpro/current_release/", index_col='UniProtKB-AC',
+                 keys=None, file_resources=None, col_rename=None, **kwargs):
+        """
+
+        Args:
+            path ():
+            index_col (str): Default 'UniProtKB-AC'.
+            keys ():
+            file_resources ():
+            col_rename ():
+            **kwargs ():
+        """
+        assert index_col is not None
+
+        if file_resources is None:
+            file_resources = {}
+            file_resources["entry.list"] = os.path.join(path, "entry.list")
+            file_resources["protein2ipr.dat.gz"] = os.path.join(path, "protein2ipr.dat.gz")
+            file_resources["interpro2go"] = os.path.join(path, "interpro2go")
+            file_resources["ParentChildTreeFile.txt"] = os.path.join(path, "ParentChildTreeFile.txt")
+
+        if any('protein2ipr' in fn for fn in file_resources):
+            assert keys is not None, "Processing protein2ipr requires user to provide `keys` as a list of " \
+                                     "`UniProtKB-AC` values"
+
+        super().__init__(path=path, file_resources=file_resources, index_col=index_col, keys=keys,
+                         col_rename=col_rename, **kwargs)
+
+    def load_dataframe(self, file_resources: Dict[str, TextIOWrapper], blocksize=None):
+        ipr_entries = pd.read_table(file_resources["entry.list"], index_col="ENTRY_AC")
+
+        ipr2go_fn = next((fn for fn in file_resources if 'interpro2go' in fn), None)
+        if ipr2go_fn:
+            ipr2go = self.parse_interpro2go(file_resources[ipr2go_fn])
+            if ipr2go is not None:
+                ipr_entries = ipr_entries.join(ipr2go.groupby('ENTRY_AC')["go_id"].unique(), on="ENTRY_AC")
+
+        ipr2go_fn = next((fn for fn in file_resources if 'interpro.xml' in fn), None)
+        if ipr2go_fn:
+            ipr_xml = pd.read_xml(file_resources[ipr2go_fn], xpath='//interpro',
+                                  compression='gzip' if ipr2go_fn.endswith('.gz') else None) \
+                .dropna(axis=1, how='all') \
+                .rename(columns={'id': 'ENTRY_AC'}) \
+                .set_index('ENTRY_AC')
+            ipr_entries = ipr_entries.join(ipr_xml, on="ENTRY_AC")
+
+        return ipr_entries
+
+    def load_network(self, file_resources) -> Tuple[nx.Graph, np.ndarray]:
+        network, node_list = None, None
+        for filename in file_resources:
+            if 'ParentChildTreeFile' in filename and isinstance(file_resources[filename], str):
+                network: nx.MultiDiGraph = self.parse_ipr_treefile(file_resources[filename])
+                node_list = np.array(network.nodes)
+
+        return network, node_list
+
+    def load_annotation(self, file_resources, blocksize=None):
+        if not any('protein2ipr' in fn for fn in file_resources):
+            return None
+
+        ipr_entries = self.data
+
+        # Use Dask
+        args = dict(names=['UniProtKB-AC', 'ENTRY_AC', 'ENTRY_NAME', 'accession', 'start', 'stop'],
+                    usecols=['UniProtKB-AC', 'ENTRY_AC', 'start', 'stop'],
+                    dtype={'UniProtKB-AC': 'category', 'ENTRY_AC': 'category', 'start': 'int8', 'stop': 'int8'},
+                    low_memory=True,
+                    blocksize=None if isinstance(blocksize, bool) else blocksize)
+        if 'protein2ipr.parquet' in file_resources:
+            annotations = dd.read_parquet(file_resources["protein2ipr.parquet"])
+        else:
+            annotations = dd.read_table(file_resources["protein2ipr.dat"], **args)
+        if self.keys is not None and self.index_col in annotations.columns:
+            annotations = annotations.loc[annotations[self.index_col].isin(self.keys)]
+        elif self.keys is not None and self.index_col == annotations.index.name:
+            annotations = annotations.loc[annotations.index.isin(self.keys)]
+        # if annotations.index.name != self.index_col:
+        #     annotations = annotations.set_index(self.index_col, sorted=True)
+        # if not annotations.known_divisions:
+        #     annotations.divisions = annotations.compute_current_divisions()
+        # Set ordering for rows and columns
+        row_order = self.keys
+        col_order = ipr_entries.index
+        row2idx = {node: i for i, node in enumerate(row_order)}
+        col2idx = {node: i for i, node in enumerate(col_order)}
+
+        def edgelist2coo(edgelist_df: DataFrame, source='UniProtKB-AC', target='ENTRY_AC') -> Optional[ssp.coo_matrix]:
+            if edgelist_df.shape[0] == 1 and edgelist_df.iloc[0, 0] == 'foo':
+                return None
+
+            if edgelist_df.index.name == source:
+                source_nodes = edgelist_df.index
+            else:
+                source_nodes = edgelist_df[source]
+
+            edgelist_df = edgelist_df.assign(row=source_nodes.map(row2idx).astype('int'),
+                                             col=edgelist_df[target].map(col2idx).astype('int'))
+
+            edgelist_df = edgelist_df.dropna(subset=['row', 'col'])
+            if edgelist_df.shape[0] == 0:
+                return None
+
+            values = np.ones(edgelist_df.index.size)
+            coo = ssp.coo_matrix((values, (edgelist_df['row'], edgelist_df['col'])),
+                                 shape=(len(row2idx), ipr_entries.index.size))
+            return coo
+
+        # Create a sparse adjacency matrix each partition, then combine them
+        adj = annotations.reduction(chunk=edgelist2coo,
+                                    aggregate=lambda x: x.dropna().sum() if not x.isna().all() else None,
+                                    meta=pd.Series([ssp.coo_matrix])).compute()
+        assert len(adj) == 1, f"len(adj) = {len(adj)}"
+        # Create a sparse matrix of UniProtKB-AC x ENTRY_AC
+        annotations = pd.DataFrame.sparse.from_spmatrix(adj[0], index=row_order, columns=col_order)
+
+        return annotations
+
+    def parse_interpro2go(self, file: str) -> pd.DataFrame:
+        def _process_line(line: str) -> Tuple[str, str, str]:
+            pos = line.find('> GO')
+            interpro_terms, go_term = line[:pos], line[pos:]
+            interpro_id, interpro_name = interpro_terms.strip().split(' ', 1)
+            go_name, go_id = go_term.split(';')
+            go_desc = go_name.strip('> GO:')
+
+            return (interpro_id.strip().split(':')[1], go_id.strip(), go_desc)
+
+        if isinstance(file, str):
+            with open(os.path.expanduser(file), 'r') as file:
+                tuples = [_process_line(line.strip()) for line in file if line[0] != '!']
+
+            ipr2go = pd.DataFrame(tuples, columns=['ENTRY_AC', "go_id", "go_desc"])
+            return ipr2go
+
+    def parse_ipr_treefile(self, lines: Union[List[str], StringIO]) -> nx.MultiDiGraph:
+        """Parse the InterPro Tree from the given file.
+        Args:
+            lines: A readable file or file-like
+        """
+        if isinstance(lines, str):
+            lines = open(os.path.expanduser(lines), 'r')
+
+        graph = nx.MultiDiGraph()
+        previous_depth, previous_name = 0, None
+        stack = [previous_name]
+
+        def count_front(s: str) -> int:
+            """Count the number of leading dashes on a string."""
+            for position, element in enumerate(s):
+                if element != '-':
+                    return position
+
+        for line in lines:
+            depth = count_front(line)
+            interpro_id, name, *_ = line[depth:].split('::')
+
+            if depth == 0:
+                stack.clear()
+                stack.append(interpro_id)
+
+                graph.add_node(interpro_id, interpro_id=interpro_id, name=name)
+
+            else:
+                if depth > previous_depth:
+                    stack.append(previous_name)
+
+                elif depth < previous_depth:
+                    del stack[-1]
+
+                parent = stack[-1]
+
+                graph.add_node(interpro_id, interpro_id=interpro_id, parent=parent, name=name)
+                graph.add_edge(parent, interpro_id, key="is_a")
+
+            previous_depth, previous_name = depth, interpro_id
+
+        lines.close()
+        return graph
+
+
+class HumanPhenotypeOntology(Ontology):
+    """Loads the Human Phenotype Ontology database from https://hpo.jax.org/app/ .
+
+        Default path: "http://geneontology.org/gene-associations/" .
+        Default file_resources: {
+            "hp.obo": "http://purl.obolibrary.org/obo/hp.obo",
+        }
+        """
+
+    COLUMNS_RENAME_DICT = {}
+
+    def __init__(
+        self,
+        path="https://hpo.jax.org/",
+        file_resources=None,
+        col_rename=COLUMNS_RENAME_DICT,
+        blocksize=0,
+        verbose=False,
+    ):
+        """
+        Handles downloading the latest Human Phenotype Ontology obo and annotation data, preprocesses them. It provide
+        functionalities to create a directed acyclic graph of Ontology terms, filter terms, and filter annotations.
+        """
+        if file_resources is None:
+            file_resources = {
+                "hp.obo": "http://purl.obolibrary.org/obo/hp.obo",
+            }
+        super().__init__(
+            path,
+            file_resources,
+            col_rename=col_rename,
+            blocksize=blocksize,
+            verbose=verbose,
+        )
+
+    def info(self):
+        print("network {}".format(nx.info(self.network)))
+    def load_network(self, file_resources):
+        for file in file_resources:
+            if ".obo" in file:
+                network = obonet.read_obo(file_resources[file])
+                network = network.reverse(copy=True)
+                node_list = np.array(network.nodes)
+        return network, node_list
+
+
+
+
+def traverse_predecessors(network, seed_node, type=["is_a", "part_of"]):
+    """
+    Returns all successor terms from seed_node by traversing the ontology network with edges == `type`.
+    Args:
+        seed_node: seed node of the traversal
+        type: the ontology type to include
+    Returns:
+        generator of list of lists for each dfs branches.
+    """
+    parents = dict(network.pred[seed_node])
+    for parent, v in parents.items():
+        if list(v.keys())[0] in type:
+            yield [parent] + list(traverse_predecessors(network, parent, type))
+
+
+def flatten(lst):
+    return sum(([x] if not isinstance(x, list) else flatten(x) for x in lst),
+               [])
+
+
+def dfs_path(graph, path):
+    node = path[-1]
+    successors = list(graph.successors(node))
+    if len(successors) > 0:
+        for child in successors:
+            yield list(dfs_path(graph, path + [child]))
+    else:
+        yield path
+
+
+def flatten_list(list_in):
+    if isinstance(list_in, list):
+        for l in list_in:
+            if isinstance(list_in[0], list):
+                for y in flatten_list(l):
+                    yield y
+            elif isinstance(list_in[0], str):
+                yield list_in
+    else:
+        yield list_in
+
+
+def filter_dfs_paths(paths_df: pd.DataFrame):
+    idx = {}
+    for col in sorted(paths_df.columns[:-1], reverse=True):
+        idx[col] = ~(paths_df[col].notnull()
+                     & paths_df[col].duplicated(keep="first")
+                     & paths_df[col + 1].isnull())
+
+    idx = pd.DataFrame(idx)
+
+    paths_df = paths_df[idx.all(axis=1)]
+    return paths_df
+
+
+def write_taxonomy(network, root_nodes, file_path):
+    """
+
+    Args:
+        network: A network with edge(i, j) where i is a node and j is a child of i.
+        root_nodes (list): a list of node names
+        file_path (str):
+    """
+    file = open(file_path, "a")
+    file.write("Root\t" + "\t".join(root_nodes) + "\n")
+
+    for root_node in root_nodes:
+        for node, children in nx.traversal.bfs_successors(network, root_node):
+            if len(children) > 0:
+                file.write(node + "\t" + "\t".join(children) + "\n")
+    file.close()