v2.0.2 - updating code for hiding taxa. now takes a file with species as in tree (gtdb issue). will now take root level in species tree by default.

Author: alex-wave

README.md

@@ -120,8 +120,8 @@ Required arguments: ``--db``, ``--oma_path``
 | [``--min_fam_size``](#markdown-header--min_fam_size)|6|Only root-HOGs with a protein count passing this threshold are used.
 | [``--min_fam_completeness``](#markdown-header--min_fam_completeness)|0.5|Only root-HOGs passing this threshold are used. The completeness of a HOG is defined as the number of observed species divided by the expected number of species at the HOG taxonomic level.
 | [``--logic``](#markdown-header--logic)|OR|Logic used between the two above arguments to filter root-HOGs. Options are "AND" or "OR".
-| [``--root_taxon``](#markdown-header--root_taxon)|LUCA|HOGs defined at, or descending from, this taxon are uses as root-HOGs.
-| [``--hidden_taxa``](#markdown-header--hidden_taxa)||The proteins from these taxa are removed before the database computation. Usage: a list of comma-separated taxa (scientific name) with underscore replacing spaces (_e.g._ Bacteria,Homo_sapiens).
+| [``--root_taxon``](#markdown-header--root_taxon)|root of speciestree.nwk|HOGs defined at, or descending from, this taxon are uses as root-HOGs.
+| [``--hidden_taxa``](#markdown-header--hidden_taxa)||The proteins from these taxa are removed before the database computation. Usage: a file containing taxa on seperate lines (scientific name). These must match EXACTLY with the node name in the tree given.
 | [``--reduced_alphabet``](#markdown-header--reduced_alphabet)||Use reduced alphabet from Linclust paper
 | [``--k``](#markdown-header--k)|6|_k_-mer length
 | [``--oma_path``](#markdown-header--oma_path)||Path to a directory with both OmaServer.h5 and speciestree.nwk

omamer/__init__.py

@@ -23,7 +23,7 @@
 from datetime import date
 
 __packagename__ = "omamer"
-__version__ = "2.0.1"
+__version__ = "2.0.2"
 __copyright__ = "(C) 2019-{:d} Victor Rossier <[email protected]> and Alex Warwick Vesztrocy <[email protected]>".format(
     date.today().year
 )

omamer/_runners.py

@@ -26,6 +26,8 @@
 def mkdb_oma(args):
     from .database import DatabaseFromOMA
     from .index import Index
+    
+    from ete3 import Tree
     import os
 
     _ensure_db_build_dependencies_available()
@@ -35,32 +37,53 @@ def mkdb_oma(args):
     LOG.info("arguments for build:")
     for (k, v) in vars(args).items():
         LOG.info(" - {}: {}".format(k, v))
+    
+    oma_db_fn = os.path.join(args.oma_path, "OmaServer.h5")
+    nwk_fn = os.path.join(args.oma_path, "speciestree.nwk")
+
+    # check if root_taxon in tree
+    t = Tree(nwk_fn, format=1, quoted_node_names=True)
+    if args.root_taxon is not None:
+        pruned_t = t.search_nodes(name=args.root_taxon)
+        if len(pruned_t) == 0:
+            raise ValueError("Unable to find {} as root taxon".format(args.root_taxon))
+        elif len(pruned_t) > 1:
+            raise ValueError("Ambiguous root taxon {} ({})".format(args.root_taxon, ", ".join(map(lambda x: x.name, pruned_t))))
+        t = pruned_t[0]
+    root_taxon = t.name
+
+    if args.hidden_taxa:
+        # load and check for hidden taxa
+        hidden_taxa = list(map(lambda x: x.rstrip(), args.hidden_taxa))
+       
+        for ht in hidden_taxa:
+            r = t.search_nodes(name=ht)
+            if len(r) == 0:
+                raise ValueError("Unable to find {} as hidden taxon".format(ht))
+            elif len(r) > 1:
+                raise ValueError("Ambiguous hidden taxon {} ({})".format(ht, ", ".join(map(lambda x: x.name, r))))
+    else:
+        hidden_taxa = []
 
     db = DatabaseFromOMA(
         args.db,
-        root_taxon=args.root_taxon,
+        root_taxon=root_taxon,
         min_fam_size=args.min_fam_size,
         logic=args.logic,
         min_fam_completeness=args.min_fam_completeness,
         include_younger_fams=True,
         mode="w",
     )
 
-    oma_db_fn = os.path.join(args.oma_path, "OmaServer.h5")
-    nwk_fn = os.path.join(args.oma_path, "speciestree.nwk")
-
     # add sequences from database
-    LOG.info("Adding sequences")
+    LOG.info("Loading sequences")
     seq_buff = db.build_database(oma_db_fn, nwk_fn)
 
-    hidden_taxa = []
-    if args.hidden_taxa:
-        hidden_taxa = [" ".join(x.split("_")) for x in args.hidden_taxa.split(",")]
-
     LOG.info("Building index")
     db.ki = Index(
         db, k=args.k, reduced_alphabet=args.reduced_alphabet, hidden_taxa=hidden_taxa
     )
+    db.ki.sp_filter
     db.ki.build_kmer_table(seq_buff)
     db.add_metadata()
     db.add_md5_hash()
@@ -207,7 +230,7 @@ def search(args):
 
 def _ensure_db_build_dependencies_available():
     try:
-        from pysais import sais
+        from PySAIS import sais
     except ImportError:
         LOG.error("To build OMAmer databases, pysais must be installed. Please ensure you installed omamer with the 'build' extra, e.g. `pip install omamer[build]`")
         import sys

omamer/database.py

@@ -696,6 +696,7 @@ def get_metadata(self):
         meta["k-mer length"] = self.ki.k
         meta["alphabet size"] = self.ki.alphabet.n
         meta["nr species"] = len(self._db_Species)
+        meta["nr indexed species"] = (len(self._db_Species) - np.sum(self.ki.sp_filter))
         meta["hidden taxa"] = self.ki.hidden_taxa
         return meta
 

omamer/index.py

@@ -20,6 +20,7 @@
     You should have received a copy of the GNU Lesser General Public License
     along with OMAmer. If not, see <http://www.gnu.org/licenses/>.
 """
+from property_manager import lazy_property
 import numba
 import numpy as np
 
@@ -46,26 +47,36 @@ def __init__(self, db, k=6, reduced_alphabet=False, hidden_taxa=[]):
 
         self.alphabet = Alphabet(n=alphabet_n)
 
-    @property
+    @lazy_property
     def sp_filter(self):
-        tax_tab = self.db._db_Taxonomy[:]
         sp_filter = np.full((len(self.db._db_Species),), False)
-        for hidden_taxon in self.hidden_taxa:
-            descendant_species = get_leaves(
-                np.searchsorted(tax_tab["ID"], hidden_taxon.encode("ascii")),
-                tax_tab,
-                self.db._db_ChildrenTax[:],
-            )
-            # case where hidden taxon is species
-            if descendant_species.size == 0:
-                sp_filter[
-                    np.searchsorted(
-                        self.db._db_Species.col("ID"), hidden_taxon.encode("ascii")
-                    )
-                ] = True
-            else:
-                for sp_off in descendant_species:
-                    sp_filter[tax_tab["SpeOff"][sp_off]] = True
+        if len(self.hidden_taxa) > 0:
+            tax_tab = self.db._db_Taxonomy[:]
+            child_tax = self.db._db_ChildrenTax[:]
+
+            if len(self.hidden_taxa) > 0:
+                LOG.debug(' - creating species filter to hide declared taxa')
+
+            for hidden_taxon in self.hidden_taxa:
+                LOG.debug('   - identifying: {}'.format(hidden_taxon))
+                taxon = np.argwhere(tax_tab["ID"] == hidden_taxon.encode("ascii")).flatten()
+                if len(taxon) == 0:
+                    raise ValueError("Can't find {} in taxonomy.".format(hidden_taxon))
+                elif len(taxon) > 1:
+                    raise ValueError("Ambiguous taxon {}.".format(hidden_taxon))
+
+                tax_ii = taxon[0]
+                sp_ii = tax_tab[tax_ii]['SpeOff']
+
+                if sp_ii >= 0:
+                    # leaf (i.e., extant species listed)
+                    LOG.debug('     - hiding {}'.format(self.db._db_Species[sp_ii]['ID'].decode('ascii')))
+                    sp_filter[sp_ii] = True
+                else:
+                    # filter all leaves below declared taxon
+                    for sp_jj in tax_tab["SpeOff"][get_leaves(tax_ii, tax_tab, child_tax)]:
+                        LOG.debug('     - hiding {}'.format(self.db._db_Species[sp_jj]['ID'].decode('ascii')))
+                        sp_filter[sp_jj] = True
         return sp_filter
 
     @property

omamer/main.py

@@ -128,14 +128,14 @@ def get_thread_count():
     )
     mkdb_parser.add_argument(
         "--root_taxon",
-        default="LUCA",
-        help="HOGs defined at, or descending from, this taxon are uses as root-HOGs.",
+        required=False,
+        help="HOGs defined at, or descending from, this taxon are uses as root-HOGs. Default is the top level in species tree.",
     )
     mkdb_parser.add_argument(
         "--hidden_taxa",
-        default="",
-        help="The proteins from these taxa are removed before the database computation. Usage: a list of comma-separated taxa (scientific name) with underscore replacing spaces (e.g. Bacteria,Homo_sapiens).",
-        type=str,
+        required=False,
+        help="Optional -- path to a file giving a list of taxa to hide the proteins during index creation only. HOGs will still exist in the database, but it will not be possible to place to them. Names must match EXACTLY to those in the given newick species tree.",
+        type=FileType("r"),
     )
     # mkdb_parser.add_argument(
     #    "--species", default="", help="Alternatively to --hidden_taxa, provide a file with species offsets in sp_tab (tmp option for scaling experiment)", type=str