Biopython reads my tree eternally long

Question

I have a nexus tree (1332 taxa) with a lot of additional data. When I tried to read it through tree = Phylo.read(treepath, "nexus"), my kernel got eternally loaded. If I abort the process, I get the following message:

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [95], in <cell line: 1>()
----> 1 tree = Phylo.read(r"D:\Downloads_hard_disk\Job\NAs_Hepatitis_A\Results\experiments\test_right\test_right_ExpPop_3\test_right_ExpPop_3.biopy.tree", "nexus")

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:60, in read(file, format, **kwargs)
     58 try:
     59     tree_gen = parse(file, format, **kwargs)
---> 60     tree = next(tree_gen)
     61 except StopIteration:
     62     raise ValueError("There are no trees in this file.") from None

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:49, in parse(file, format, **kwargs)
     34 """Parse a file iteratively, and yield each of the trees it contains.
     35 
     36 If a file only contains one tree, this still returns an iterable object that
   (...)
     46 
     47 """
     48 with File.as_handle(file) as fp:
---> 49     yield from getattr(supported_formats[format], "parse")(fp, **kwargs)

File ~\miniconda3\lib\site-packages\Bio\Phylo\NexusIO.py:40, in parse(handle)
     32 def parse(handle):
     33     """Parse the trees in a Nexus file.
     34 
     35     Uses the old Nexus.Trees parser to extract the trees, converts them back to
   (...)
     38     eventually change Nexus to use the new NewickIO parser directly.)
     39     """
---> 40     nex = Nexus.Nexus(handle)
     42     # NB: Once Nexus.Trees is modified to use Tree.Newick objects, do this:
     43     # return iter(nex.trees)
     44     # Until then, convert the Nexus.Trees.Tree object hierarchy:
     45     def node2clade(nxtree, node):

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:668, in Nexus.__init__(self, input)
    665 self.options["gapmode"] = "missing"
    667 if input:
--> 668     self.read(input)
    669 else:
    670     self.read(DEFAULTNEXUS)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:718, in Nexus.read(self, input)
    716     break
    717 if title in KNOWN_NEXUS_BLOCKS:
--> 718     self._parse_nexus_block(title, contents)
    719 else:
    720     self._unknown_nexus_block(title, contents)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:759, in Nexus._parse_nexus_block(self, title, contents)
    757 for line in block.commandlines:
    758     try:
--> 759         getattr(self, "_" + line.command)(line.options)
    760     except AttributeError:
    761         raise NexusError("Unknown command: %s " % line.command) from None

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:1181, in Nexus._tree(self, options)
   1179     elif special == "W":
   1180         weight = float(value)
-> 1181 tree = Tree(name=name, weight=weight, rooted=rooted, tree=opts.rest().strip())
   1182 # if there's an active translation table, translate
   1183 if self.translate:

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:82, in Tree.__init__(self, tree, weight, rooted, name, data, values_are_support, max_support)
     80 # there's discrepancy whether newick allows semicolons et the end
     81 tree = tree.rstrip(";")
---> 82 subtree_info, base_info = self._parse(tree)
     83 root.data = self._add_nodedata(root.data, [[], base_info])
     84 self._add_subtree(parent_id=root.id, tree=subtree_info)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:123, in Tree._parse(self, tree)
    121 elif not incomment and tree[p] == ")":
    122     plevel -= 1
--> 123 elif tree[p:].startswith(NODECOMMENT_START):
    124     incomment = True
    125 elif incomment and tree[p] == NODECOMMENT_END:

KeyboardInterrupt: 

I am interested what is the specific reason of this: either the ineffective algorithm or a bug in Biopython? Is there a way to read the tree as Bio.Phylo.Newick.Tree object? Or should I completely use another package? In the last case, could you recommend effective one?

Answer 1

The reason is the RAM has been exceeded for Biopython's capabilities and you are encountering a 'RAM-bottleneck'. I've listed the options from 1 or 3 and recommend option 3c.

The options are use:

1. machine with much greater RAM
2. acceleration coding techniques in Python (such parallelisation)
3. Beast tools written in Java

I assume this calculation is being performed on your laptop because in theory <1500 isn't much, but agreed there are ~1500 taxa.

Option 1 I personally don't think you are very far into the bottle-neck. Certainly a Cloud computer (e.g. GCP) could easily perform this calculation and maybe your the headnode of your local cluster. Usually the headnode has the most RAM.

Option 2 Code acceleration is worth keeping in mind because this achieves >> 1 log performance. Parallisation will solve this using an optimal chunk-size. However, large amounts of code can be needed. Code optimisation - Python is very sensitive to this because it uses C-libraries - for example, use map and reduce and avoid ANY regex in your Python code because IMO its a key weakness of Python. Pandas iterrows is famously slow. If the Biopython library is the problem code efficiency isn't an option (unless you want to rewrite Biopython). Label compression or header recoding techniques work, I'm not sure whether this has been done in Beast output already(?).

Option 3 The tools for dealing with this which are written and optimised to achieve rapid output are:

1. Tracer - does lots of analysis on MCMC output. Its really geared towards assessing convergence, but it can do a lot of stuff.
2. TreeAnnotator - makes a summary tree from the MCMC output you have

3. TreeStat2 <- I think this is what you are looking for.