codonlogo: corebio/seq_io/_nexus/Trees.py comparison

comparison corebio/seq_io/_nexus/Trees.py @ 0:c55bdc2fb9fa

Uploaded

author	davidmurphy
date	Thu, 27 Oct 2011 12:09:09 -0400
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:c55bdc2fb9fa
+#
+# Trees.py
+#
+# Copyright 2005 by Frank Kauff & Cymon J. Cox. All rights reserved.
+# This code is part of the Biopython distribution and governed by its
+# license. Please see the LICENSE file that should have been included
+# as part of this package.
+#
+# Tree class handles phylogenetic trees. Provides a set of methods to read and write newick-format tree
+# descriptions, get information about trees (monphyly of taxon sets, congruence between trees, common ancestors,...)
+# and to manipulate trees (reroot trees, split terminal nodes).
+#
+# Bug reports welcome: fkauff@duke.edu
+#
+import sys, random, sets
+import Nodes
+PRECISION_BRANCHLENGTH=6
+PRECISION_SUPPORT=6
+class TreeError(Exception): pass
+class NodeData:
+"""Stores tree-relevant data associated with nodes (e.g. branches or otus)."""
+def __init__(self,taxon=None,branchlength=0.0,support=None):
+self.taxon=taxon
+self.branchlength=branchlength
+self.support=support
+class Tree(Nodes.Chain):
+"""Represents a tree using a chain of nodes with on predecessor (=ancestor)
+and multiple successors (=subclades).
+"""
+# A newick tree is parsed into nested list and then converted to a node list in two stages
+# mostly due to historical reasons. This could be done in one swoop). Note: parentheses ( ) and
+# colon : are not allowed in taxon names. This is against NEXUS standard, but makes life much
+# easier when parsing trees.
+## NOTE: Tree should store its data class in something like self.dataclass=data,
+## so that nodes that are generated have easy access to the data class
+## Some routines use automatically NodeData, this needs to be more concise
+def __init__(self,tree=None,weight=1.0,rooted=False,name='',data=NodeData,values_are_support=False,max_support=1.0):
+"""Ntree(self,tree)."""
+Nodes.Chain.__init__(self)
+self.dataclass=data
+self.__values_are_support=values_are_support
+self.max_support=max_support
+self.weight=weight
+self.rooted=rooted
+self.name=name
+root=Nodes.Node(data())
+self.add(root)
+self.root=root.id
+if tree:    # use the tree we have
+# if Tree is called from outside Nexus parser, we need to get rid of linebreaks, etc
+tree=tree.strip().replace('\n','').replace('\r','')
+# there's discrepancy whether newick allows semicolons et the end
+tree=tree.rstrip(';')
+self._add_subtree(parent_id=root.id,tree=self._parse(tree)[0])
+def _parse(self,tree):
+"""Parses (a,b,c...)[[[xx]:]yy] into subcomponents and travels down recursively."""
+if tree.count('(')!=tree.count(')'):
+raise TreeError, 'Parentheses do not match in (sub)tree: '+tree
+if tree.count('(')==0: # a leaf
+colon=tree.rfind(':')
+if colon>-1:
+return [tree[:colon],self._get_values(tree[colon+1:])]
+else:
+return [tree,[None]]
+else:
+closing=tree.rfind(')')
+val=self._get_values(tree[closing+1:])
+if not val:
+val=[None]
+subtrees=[]
+plevel=0
+prev=1
+for p in range(1,closing):
+if tree[p]=='(':
+plevel+=1
+elif tree[p]==')':
+plevel-=1
+elif tree[p]==',' and plevel==0:
+subtrees.append(tree[prev:p])
+prev=p+1
+subtrees.append(tree[prev:closing])
+subclades=[self._parse(subtree) for subtree in subtrees]
+return [subclades,val]
+def _add_subtree(self,parent_id=None,tree=None):
+"""Adds leaf or tree (in newick format) to a parent_id. (self,parent_id,tree)."""
+if parent_id is None:
+raise TreeError('Need node_id to connect to.')
+for st in tree:
+if type(st[0])==list: # it's a subtree
+nd=self.dataclass()
+if len(st[1])>=2: # if there's two values, support comes first. Is that always so?
+nd.support=st[1][0]
+if st[1][1] is not None:
+nd.branchlength=st[1][1]
+elif len(st[1])==1: # otherwise it could be real branchlengths or support as branchlengths
+if not self.__values_are_support: # default
+if st[1][0] is not None:
+nd.branchlength=st[1][0]
+else:
+nd.support=st[1][0]
+sn=Nodes.Node(nd)
+self.add(sn,parent_id)
+self._add_subtree(sn.id,st[0])
+else: # it's a leaf
+nd=self.dataclass()
+nd.taxon=st[0]
+if len(st)>1:
+if len(st[1])>=2: # if there's two values, support comes first. Is that always so?
+nd.support=st[1][0]
+if st[1][1] is not None:
+nd.branchlength=st[1][1]
+elif len(st[1])==1: # otherwise it could be real branchlengths or support as branchlengths
+if not self.__values_are_support: # default
+if st[1][0] is not None:
+nd.branchlength=st[1][0]
+else:
+nd.support=st[1][0]
+leaf=Nodes.Node(nd)
+self.add(leaf,parent_id)
+def _get_values(self, text):
+"""Extracts values (support/branchlength) from xx[:yyy], xx."""
+if text=='':
+return None
+return [float(t) for t in text.split(':') if t.strip()]
+def _walk(self,node=None):
+"""Return all node_ids downwards from a node."""
+if node is None:
+node=self.root
+for n in self.node(node).succ:
+yield n
+for sn in self._walk(n):
+yield sn
+def node(self,node_id):
+"""Return the instance of node_id.
+node = node(self,node_id)
+"""
+if node_id not in self.chain:
+raise TreeError('Unknown node_id: %d' % node_id)
+return self.chain[node_id]
+def split(self,parent_id=None,n=2,branchlength=1.0):
+"""Speciation: generates n (default two) descendants of a node.
+[new ids] = split(self,parent_id=None,n=2,branchlength=1.0):
+"""
+if parent_id is None:
+raise TreeError('Missing node_id.')
+ids=[]
+parent_data=self.chain[parent_id].data
+for i in range(n):
+node=Nodes.Node()
+if parent_data:
+node.data=self.dataclass()
+# each node has taxon and branchlength attribute
+if parent_data.taxon:
+node.data.taxon=parent_data.taxon+str(i)
+node.data.branchlength=branchlength
+ids.append(self.add(node,parent_id))
+return ids
+def search_taxon(self,taxon):
+"""Returns the first matching taxon in self.data.taxon. Not restricted to terminal nodes.
+node_id = search_taxon(self,taxon)
+"""
+for id,node in self.chain.items():
+if node.data.taxon==taxon:
+return id
+return None
+def prune(self,taxon):
+"""Prunes a terminal taxon from the tree.
+id_of_previous_node = prune(self,taxon)
+If taxon is from a bifurcation, the connectiong node will be collapsed
+and its branchlength added to remaining terminal node. This might be no
+longer a meaningful value'
+"""
+id=self.search_taxon(taxon)
+if id is None:
+raise TreeError('Taxon not found: %s' % taxon)
+elif id not in self.get_terminals():
+raise TreeError('Not a terminal taxon: %s' % taxon)
+else:
+prev=self.unlink(id)
+self.kill(id)
+if not prev==self.root and len(self.node(prev).succ)==1:
+succ=self.node(prev).succ[0]
+new_bl=self.node(prev).data.branchlength+self.node(succ).data.branchlength
+self.collapse(prev)
+self.node(succ).data.branchlength=new_bl
+return prev
+def get_taxa(self,node_id=None):
+"""Return a list of all otus downwards from a node (self, node_id).
+nodes = get_taxa(self,node_id=None)
+"""
+if node_id is None:
+node_id=self.root
+if node_id not in self.chain:
+raise TreeError('Unknown node_id: %d.' % node_id)
+if self.chain[node_id].succ==[]:
+if self.chain[node_id].data:
+return [self.chain[node_id].data.taxon]
+else:
+return None
+else:
+list=[]
+for succ in self.chain[node_id].succ:
+list.extend(self.get_taxa(succ))
+return list
+def get_terminals(self):
+"""Return a list of all terminal nodes."""
+return [i for i in self.all_ids() if self.node(i).succ==[]]
+def sum_branchlength(self,root=None,node=None):
+"""Adds up the branchlengths from root (default self.root) to node.
+sum = sum_branchlength(self,root=None,node=None)
+"""
+if root is None:
+root=self.root
+if node is None:
+raise TreeError('Missing node id.')
+blen=0.0
+while node is not None and node is not root:
+blen+=self.node(node).data.branchlength
+node=self.node(node).prev
+return blen
+def set_subtree(self,node):
+"""Return subtree as a set of nested sets.
+sets = set_subtree(self,node)
+"""
+if self.node(node).succ==[]:
+return self.node(node).data.taxon
+else:
+return sets.Set([self.set_subtree(n) for n in self.node(node).succ])
+def is_identical(self,tree2):
+"""Compare tree and tree2 for identity.
+result = is_identical(self,tree2)
+"""
+return self.set_subtree(self.root)==tree2.set_subtree(tree2.root)
+def is_compatible(self,tree2,threshold,strict=True):
+"""Compares branches with support>threshold for compatibility.
+result = is_compatible(self,tree2,threshold)
+"""
+# check if both trees have the same set of taxa. strict=True enforces this.
+missing2=sets.Set(self.get_taxa())-sets.Set(tree2.get_taxa())
+missing1=sets.Set(tree2.get_taxa())-sets.Set(self.get_taxa())
+if strict and (missing1 or missing2):
+if missing1:
+print 'Taxon/taxa %s is/are missing in tree %s' % (','.join(missing1) , self.name)
+if missing2:
+print 'Taxon/taxa %s is/are missing in tree %s' % (','.join(missing2) , tree2.name)
+raise TreeError, 'Can\'t compare trees with different taxon compositions.'
+t1=[(sets.Set(self.get_taxa(n)),self.node(n).data.support) for n in self.all_ids() if \
+self.node(n).succ and\
+(self.node(n).data and self.node(n).data.support and self.node(n).data.support>=threshold)]
+t2=[(sets.Set(tree2.get_taxa(n)),tree2.node(n).data.support) for n in tree2.all_ids() if \
+tree2.node(n).succ and\
+(tree2.node(n).data and tree2.node(n).data.support and tree2.node(n).data.support>=threshold)]
+conflict=[]
+for (st1,sup1) in t1:
+for (st2,sup2) in t2:
+if not st1.issubset(st2) and not st2.issubset(st1):                     # don't hiccup on upstream nodes
+intersect,notin1,notin2=st1 & st2, st2-st1, st1-st2                 # all three are non-empty sets
+# if notin1==missing1 or notin2==missing2  <==> st1.issubset(st2) or st2.issubset(st1) ???
+if intersect and not (notin1.issubset(missing1) or notin2.issubset(missing2)):         # omit conflicts due to missing taxa
+conflict.append((st1,sup1,st2,sup2,intersect,notin1,notin2))
+return conflict
+def common_ancestor(self,node1,node2):
+"""Return the common ancestor that connects to nodes.
+node_id = common_ancestor(self,node1,node2)
+"""
+l1=[self.root]+self.trace(self.root,node1)
+l2=[self.root]+self.trace(self.root,node2)
+return [n for n in l1 if n in l2][-1]
+def distance(self,node1,node2):
+"""Add and return the sum of the branchlengths between two nodes.
+dist = distance(self,node1,node2)
+"""
+ca=self.common_ancestor(node1,node2)
+return self.sum_branchlength(ca,node1)+self.sum_branchlength(ca,node2)
+def is_monophyletic(self,taxon_list):
+"""Return node_id of common ancestor if taxon_list is monophyletic, -1 otherwise.
+result = is_monophyletic(self,taxon_list)
+"""
+if isinstance(taxon_list,str):
+taxon_set=sets.Set([taxon_list])
+else:
+taxon_set=sets.Set(taxon_list)
+node_id=self.root
+while 1:
+subclade_taxa=sets.Set(self.get_taxa(node_id))
+if subclade_taxa==taxon_set:                                        # are we there?
+return node_id
+else:                                                               # check subnodes
+for subnode in self.chain[node_id].succ:
+if sets.Set(self.get_taxa(subnode)).issuperset(taxon_set):  # taxon_set is downstream
+node_id=subnode
+break   # out of for loop
+else:
+return -1   # taxon set was not with successors, for loop exhausted
+def is_bifurcating(self,node=None):
+"""Return True if tree downstream of node is strictly bifurcating."""
+if not node:
+node=self.root
+if node==self.root and len(self.node(node).succ)==3: #root can be trifurcating, because it has no ancestor
+return self.is_bifurcating(self.node(node).succ[0]) and \
+self.is_bifurcating(self.node(node).succ[1]) and \
+self.is_bifurcating(self.node(node).succ[2])
+if len(self.node(node).succ)==2:
+return self.is_bifurcating(self.node(node).succ[0]) and self.is_bifurcating(self.node(node).succ[1])
+elif len(self.node(node).succ)==0:
+return True
+else:
+return False
+def branchlength2support(self):
+"""Move values stored in data.branchlength to data.support, and set branchlength to 0.0
+This is necessary when support has been stored as branchlength (e.g. paup), and has thus
+been read in as branchlength.
+"""
+for n in self.chain.keys():
+self.node(n).data.support=self.node(n).data.branchlength
+self.node(n).data.branchlength=0.0
+def convert_absolute_support(self,nrep):
+"""Convert absolute support (clade-count) to rel. frequencies.
+Some software (e.g. PHYLIP consense) just calculate how often clades appear, instead of
+calculating relative frequencies."""
+for n in self._walk():
+if self.node(n).data.support:
+self.node(n).data.support/=float(nrep)
+def randomize(self,ntax=None,taxon_list=None,branchlength=1.0,branchlength_sd=None,bifurcate=True):
+"""Generates a random tree with ntax taxa and/or taxa from taxlabels.
+new_tree = randomize(self,ntax=None,taxon_list=None,branchlength=1.0,branchlength_sd=None,bifurcate=True)
+Trees are bifurcating by default. (Polytomies not yet supported).
+"""
+if not ntax and taxon_list:
+ntax=len(taxon_list)
+elif not taxon_list and ntax:
+taxon_list=['taxon'+str(i+1) for i in range(ntax)]
+elif not ntax and not taxon_list:
+raise TreeError('Either numer of taxa or list of taxa must be specified.')
+elif ntax<>len(taxon_list):
+raise TreeError('Length of taxon list must correspond to ntax.')
+# initiate self with empty root
+self.__init__()
+terminals=self.get_terminals()
+# bifurcate randomly at terminal nodes until ntax is reached
+while len(terminals)<ntax:
+newsplit=random.choice(terminals)
+new_terminals=self.split(parent_id=newsplit,branchlength=branchlength)
+# if desired, give some variation to the branch length
+if branchlength_sd:
+for nt in new_terminals:
+bl=random.gauss(branchlength,branchlength_sd)
+if bl<0:
+bl=0
+self.node(nt).data.branchlength=bl
+terminals.extend(new_terminals)
+terminals.remove(newsplit)
+# distribute taxon labels randomly
+random.shuffle(taxon_list)
+for (node,name) in zip(terminals,taxon_list):
+self.node(node).data.taxon=name
+def display(self):
+"""Quick and dirty lists of all nodes."""
+table=[('#','taxon','prev','succ','brlen','blen (sum)','support')]
+for i in self.all_ids():
+n=self.node(i)
+if not n.data:
+table.append((str(i),'-',str(n.prev),str(n.succ),'-','-','-'))
+else:
+tx=n.data.taxon
+if not tx:
+tx='-'
+blength=n.data.branchlength
+if blength is None:
+blength='-'
+sum_blength='-'
+else:
+sum_blength=self.sum_branchlength(node=i)
+support=n.data.support
+if support is None:
+support='-'
+table.append((str(i),tx,str(n.prev),str(n.succ),blength,sum_blength,support))
+print '\n'.join(['%3s %32s %15s %15s %8s %10s %8s' % l for l in table])
+print '\nRoot: ',self.root
+def to_string(self,support_as_branchlengths=False,branchlengths_only=False,plain=True,plain_newick=False):
+"""Return a paup compatible tree line.
+to_string(self,support_as_branchlengths=False,branchlengths_only=False,plain=True)
+"""
+# if there's a conflict in the arguments, we override plain=True
+if support_as_branchlengths or branchlengths_only:
+plain=False
+self.support_as_branchlengths=support_as_branchlengths
+self.branchlengths_only=branchlengths_only
+self.plain=plain
+def make_info_string(data,terminal=False):
+"""Creates nicely formatted support/branchlengths."""
+# CHECK FORMATTING
+if self.plain: # plain tree only. That's easy.
+return ''
+elif self.support_as_branchlengths: # support as branchlengths (eg. PAUP), ignore actual branchlengths
+if terminal:    # terminal branches have 100% support
+return ':%1.2f' % self.max_support
+else:
+return ':%1.2f' % (data.support)
+elif self.branchlengths_only: # write only branchlengths, ignore support
+return ':%1.5f' % (data.branchlength)
+else:   # write suport and branchlengths (e.g. .con tree of mrbayes)
+if terminal:
+return ':%1.5f' % (data.branchlength)
+else:
+if data.branchlength is not None and data.support is not None:  # we have blen and suppport
+return '%1.2f:%1.5f' % (data.support,data.branchlength)
+elif data.branchlength is not None:                             # we have only blen
+return '0.00000:%1.5f' % (data.branchlength)
+elif data.support is not None:                                  # we have only support
+return '%1.2f:0.00000' % (data.support)
+else:
+return '0.00:0.00000'
+def newickize(node):
+"""Convert a node tree to a newick tree recursively."""
+if not self.node(node).succ:    #terminal
+return self.node(node).data.taxon+make_info_string(self.node(node).data,terminal=True)
+else:
+return '(%s)%s' % (','.join(map(newickize,self.node(node).succ)),make_info_string(self.node(node).data))
+return subtree
+treeline=['tree']
+if self.name:
+treeline.append(self.name)
+else:
+treeline.append('a_tree')
+treeline.append('=')
+if self.weight<>1:
+treeline.append('[&W%s]' % str(round(float(self.weight),3)))
+if self.rooted:
+treeline.append('[&R]')
+treeline.append('(%s)' % ','.join(map(newickize,self.node(self.root).succ)))
+treeline.append(';')
+if plain_newick:
+return treeline[-2]
+else:
+return ' '.join(treeline)
+def __str__(self):
+"""Short version of to_string(), gives plain tree"""
+return self.to_string(plain=True)
+def unroot(self):
+"""Defines a unrooted Tree structure, using data of a rooted Tree."""
+# travel down the rooted tree structure and save all branches and the nodes they connect
+def _get_branches(node):
+branches=[]
+for b in self.node(node).succ:
+branches.append([node,b,self.node(b).data.branchlength,self.node(b).data.support])
+branches.extend(_get_branches(b))
+return branches
+self.unrooted=_get_branches(self.root)
+# if root is bifurcating, then it is eliminated
+if len(self.node(self.root).succ)==2:
+# find the two branches that connect to root
+rootbranches=[b for b in self.unrooted if self.root in b[:2]]
+b1=self.unrooted.pop(self.unrooted.index(rootbranches[0]))
+b2=self.unrooted.pop(self.unrooted.index(rootbranches[1]))
+# Connect them two each other. If both have support, it should be identical (or one set to None?).
+# If both have branchlengths, they will be added
+newbranch=[b1[1],b2[1],b1[2]+b2[2]]
+if b1[3] is None:
+newbranch.append(b2[3]) # either None (both rootbranches are unsupported) or some support
+elif b2[3] is None:
+newbranch.append(b1[3]) # dito
+elif b1[3]==b2[3]:
+newbranch.append(b1[3]) # identical support
+elif b1[3]==0 or b2[3]==0:
+newbranch.append(b1[3]+b2[3]) # one is 0, take the other
+else:
+raise TreeError, 'Support mismatch in bifurcating root: %f, %f' % (float(b1[3]),float(b2[3]))
+self.unrooted.append(newbranch)
+def root_with_outgroup(self,outgroup=None):
+def _connect_subtree(parent,child):
+"""Hook subtree starting with node child to parent."""
+for i,branch in enumerate(self.unrooted):
+if parent in branch[:2] and child in branch[:2]:
+branch=self.unrooted.pop(i)
+break
+else:
+raise TreeError, 'Unable to connect nodes for rooting: nodes %d and %d are not connected' % (parent,child)
+self.link(parent,child)
+self.node(child).data.branchlength=branch[2]
+self.node(child).data.support=branch[3]
+#now check if there are more branches connected to the child, and if so, connect them
+child_branches=[b for b in self.unrooted if child in b[:2]]
+for b in child_branches:
+if child==b[0]:
+succ=b[1]
+else:
+succ=b[0]
+_connect_subtree(child,succ)
+# check the outgroup we're supposed to root with
+if outgroup is None:
+return self.root
+outgroup_node=self.is_monophyletic(outgroup)
+if outgroup_node==-1:
+return -1
+# if tree is already rooted with outgroup on a bifurcating root,
+# or the outgroup includes all taxa on the tree, then we're fine
+if (len(self.node(self.root).succ)==2 and outgroup_node in self.node(self.root).succ) or outgroup_node==self.root:
+return self.root
+self.unroot()
+# now we find the branch that connects outgroup and ingroup
+#print self.node(outgroup_node).prev
+for i,b in enumerate(self.unrooted):
+if outgroup_node in b[:2] and self.node(outgroup_node).prev in b[:2]:
+root_branch=self.unrooted.pop(i)
+break
+else:
+raise TreeError, 'Unrooted and rooted Tree do not match'
+if outgroup_node==root_branch[1]:
+ingroup_node=root_branch[0]
+else:
+ingroup_node=root_branch[1]
+# now we destroy the old tree structure, but keep node data. Nodes will be reconnected according to new outgroup
+for n in self.all_ids():
+self.node(n).prev=None
+self.node(n).succ=[]
+# now we just add both subtrees (outgroup and ingroup) branch for branch
+root=Nodes.Node(data=NodeData())            # new root
+self.add(root)                              # add to tree description
+self.root=root.id                           # set as root
+self.unrooted.append([root.id,ingroup_node,root_branch[2],root_branch[3]])  # add branch to ingroup to unrooted tree
+self.unrooted.append([root.id,outgroup_node,0.0,0.0])   # add branch to outgroup to unrooted tree
+_connect_subtree(root.id,ingroup_node)      # add ingroup
+_connect_subtree(root.id,outgroup_node)     # add outgroup
+# if theres still a lonely node in self.chain, then it's the old root, and we delete it
+oldroot=[i for i in self.all_ids() if self.node(i).prev is None and i!=self.root]
+if len(oldroot)>1:
+raise TreeError, 'Isolated nodes in tree description: %s' % ','.join(oldroot)
+elif len(oldroot)==1:
+self.kill(oldroot[0])
+return self.root
+def consensus(trees, threshold=0.5,outgroup=None):
+"""Compute a majority rule consensus tree of all clades with relative frequency>=threshold from a list of trees."""
+total=len(trees)
+if total==0:
+return None
+# shouldn't we make sure that it's NodeData or subclass??
+dataclass=trees[0].dataclass
+max_support=trees[0].max_support
+clades={}
+#countclades={}
+alltaxa=sets.Set(trees[0].get_taxa())
+# calculate calde frequencies
+c=0
+for t in trees:
+c+=1
+#if c%50==0:
+#    print c
+if alltaxa!=sets.Set(t.get_taxa()):
+raise TreeError, 'Trees for consensus must contain the same taxa'
+t.root_with_outgroup(outgroup=outgroup)
+for st_node in t._walk(t.root):
+subclade_taxa=t.get_taxa(st_node)
+subclade_taxa.sort()
+subclade_taxa=str(subclade_taxa) # lists are not hashable
+if subclade_taxa in clades:
+clades[subclade_taxa]+=float(t.weight)/total
+else:
+clades[subclade_taxa]=float(t.weight)/total
+#if subclade_taxa in countclades:
+#    countclades[subclade_taxa]+=t.weight
+#else:
+#    countclades[subclade_taxa]=t.weight
+# weed out clades below threshold
+for (c,p) in clades.items():
+if p<threshold:
+del clades[c]
+# create a tree with a root node
+consensus=Tree(name='consensus_%2.1f' % float(threshold),data=dataclass)
+# each clade needs a node in the new tree, add them as isolated nodes
+for (c,s) in clades.items():
+node=Nodes.Node(data=dataclass())
+node.data.support=s
+node.data.taxon=sets.Set(eval(c))
+consensus.add(node)
+# set root node data
+consensus.node(consensus.root).data.support=None
+consensus.node(consensus.root).data.taxon=alltaxa
+# we sort the nodes by no. of taxa in the clade, so root will be the last
+consensus_ids=consensus.all_ids()
+consensus_ids.sort(lambda x,y:len(consensus.node(x).data.taxon)-len(consensus.node(y).data.taxon))
+# now we just have to hook each node to the next smallest node that includes all taxa of the current
+for i,current in enumerate(consensus_ids[:-1]): # skip the last one which is the root
+#print '----'
+#print 'current: ',consensus.node(current).data.taxon
+# search remaining nodes
+for parent in consensus_ids[i+1:]:
+#print 'parent: ',consensus.node(parent).data.taxon
+if consensus.node(parent).data.taxon.issuperset(consensus.node(current).data.taxon):
+break
+else:
+sys.exit('corrupt tree structure?')
+# internal nodes don't have taxa
+if len(consensus.node(current).data.taxon)==1:
+consensus.node(current).data.taxon=consensus.node(current).data.taxon.pop()
+# reset the support for terminal nodes to maximum
+#consensus.node(current).data.support=max_support
+else:
+consensus.node(current).data.taxon=None
+consensus.link(parent,current)
+# eliminate root taxon name
+consensus.node(consensus_ids[-1]).data.taxon=None
+return consensus

Mercurial > repos > davidmurphy > codonlogo

comparison corebio/seq_io/_nexus/Trees.py @ 0:c55bdc2fb9fa