The Phylogenetic Bootstrap

image here
The image of the bootstrap is, or ought to be, a familiar one to every computer user. It is occurs in the phrase "boot the computer," which draws on James Joyce's phrase:

"Ladies who like distinctive underclothing should, and every welltailored man must, trying to make the gap wider between them by innuendo and give more of a genuine filip to acts of impropriety between the two, she unbuttoned his and then he untied her, mind the pin, whereas savages in the cannibal islands, say, at ninety degrees in the shade not caring a continental. However, reverting to the original, there were on the other hand others who had forced their way to the top from the lowest rung by the aid of their bootstraps. Sheer force of natural genius, that. With brains, sir."

The long quote here courtesy of Project Gutenburg (my post here).

Wikipedia has a different provenance (here).

The bootstrap is a relatively new statistical technique. It was applied by Felsentstein to phylogenies in 1985 (Evolution 39(4):783-791); he cites B. Efron (1979 Ann. Statist. 7:1-26). Neither paper is in PubMed, though Efron discusses Felsenstein's approach here. (A disclaimer: I have yet to actually read these, though I certainly plan to.)

In the phylogenetics context, the bootstrap simply samples columns from an alignment, with replacement, as many times as there are columns. A tree is constructed using the sampled data and evaluated for a node or clade of interest. The process is repeated many times (100 or 1000) and the number of times the majority rule clade is recovered is noted.

Until recently, I had wrongly conceived of the bootstrap as a sort of data sufficiency test, a measure of where we might be on a hypothetical saturation curve for data acquisition. But that's not it. What the bootstrap really is about is an assessment of what fraction of the data points support the hypothesis. As Felsenstein says:

"Majority-rule consensus trees can be used to construct a phylogeny showing all of the inferred monophyletic groups that occurred in a majority of the bootstrap samples. If a group shows up 95% of the time or more, the evidence for it is taken to be statistically significant."

I found the value quite surprising. I hadn't conceived of the bootstrap value as being like a p-value. Often I've looked at numbers for nodes/clades in phylogenies and thought, "70%, that's not so bad.."

It doesn't really matter how muddled my thinking is, what you might find useful here is the result of a simulation of the bootstrap in Python, including my testing it using tools PyCogent.

[UPDATE: A little more detail: I study the simplest possible case of a tree (A,B),(C,D). Data supporting the correct tree is in the form: 'AAGG'---that is otu A has 'A', otu B has 'A', otu C has 'G' and otu D has 'G'. Data not supporting the correct tree is in the form: 'TCTC'. We specify the fraction of columns supporting ('pro') and not supporting ('con') to be contained in the alignment. ]

The code to draw the bootstrap samples is trivial.

For the analysis, I ended up using the fact that the number of "connecting edges" for nodes A and B is 3 for the correct tree (themselves, plus their parent). I couldn't get the function sameTopology() to evaluate to True in comparisons with the prototype: ((A,B),(C,D)). A most surprising result of the simulation was how little effect substantial contradictory evidence has, and how steep is the curve. Here, the fraction of times the correct tree is recovered is almost never less than 1.0 until about 40% of the values in the sample are from contradictory alignments. That's a whole lot of contradictory data!

A secondary observation was that "uninformative" positions in the alignment don't seem to alter the value obtained at all.



code listing:

import random, sys from cogent import LoadSeqs, DNA, LoadTree from cogent.phylo import distance, nj from cogent.evolve.models import F81, JC69 import numpy as np import matplotlib.pyplot as plt debug = False # true tree is: ((A,B),(C,D)) # alignment that is consistent: AAGG # alignment that is not: AGAG # uninformative alignment: AAAA tr = LoadTree(treestring='((A,B),(C,D));') def makeAlignment(pro,con,neutral): seqs = list() for otu in 'ABCD': temp = ['A'*neutral] if otu in 'AB': temp.append('A'*pro) else: temp.append('G'*pro) if otu in 'AC': temp.append('T'*con) else: temp.append('C'*con) seqs.append(''.join(temp)) return seqs def sample(seqs): N = len(seqs[0]) iL = [random.choice(range(N)) for j in range(N)] rL = list() for seq in seqs: temp = [seq[i] for i in iL] rL.append(''.join(temp)) return rL def evaluate_tree(aln): d = distance.EstimateDistances(aln, submodel=JC69()) d.run(show_progress=False) njtree = nj.nj(d.getPairwiseDistances()) if debug: print d print njtree.asciiArt() print njtree.sameTopology(tr) for otu in 'BCD': print njtree.getConnectingEdges('A',otu) L = njtree.getConnectingEdges('A','B') return len(L) == 3 R = [35,40,42,44] + range(44,57) + [58,60,65] C = [] for con in R: pro = 100 - con seqs = makeAlignment(pro,con,0) if debug: print '\n'.join(seqs) count = 0 reps = 100 for i in range(reps): #print i s = sample(seqs) #print '\n'.join(s) dna = dict(zip('ABCD',s)) aln = LoadSeqs(data=dna, moltype=DNA) if evaluate_tree(aln): count += 1 C.append(count*1.0/reps) print con, count, reps plt.scatter(R,np.array(C),s=250) ax = plt.axes() ax.set_xlim(30,70) ax.set_ylim(-0.05,1.05) plt.savefig('example.png')

Phylogenetic tree surgery 1

One of the biggest problems in phylogenetic analysis is the very large number of possible trees. I don't think I've posted about that yet, but you probably appreciate that for N species the number of trees is approximately 10^N. (Speaking very roughly). For 25 species we're talking more than a mole of trees. That's a lot!

One consequence is that we can never evaluate all the possible trees for any problem of significant size. Furthermore, if we have a tree that we like (it's the one with the maximum likelihood that we've seen so far), it becomes important to look at its near neighbors in "tree space" and see if they are any better. Felsenstein has a great discussion of this in Chapter 4.

I want to study a bit about the basic methods of tree surgery. My hunch is that representing trees in the way that was introduced previously (here), as lists of connections to each internal node, will make this pretty easy. My goal in this post is simply to present one example of each of the methods that I know about, and show the reorganization of the list that results. The examples are shown in the figures before, which were constructed in Keynote. We'll explore these in future posts, as time permits. The plan is to develop rules to emulate the types of list reorganization observed here, and see if the structure of the resulting trees can be explained by that type of surgery.

In each example, the list entries that change are highlighted.

In NNI we pick a branch, and rearrange the quartet of clades to form a new tree. There are two possible rearrangements. Here is one of them for the 1-3 branch.



In SPR we imagine cutting in the middle of one branch, and then connecting one of the pieces to the middle of a second branch in the other piece. I had to relabel the internal nodes for this one, I did it in a way that makes the changes to the list easiest.



In TBR we cut a branch and dissolve the stubs, then pick one branch in each of the smaller pieces, and connect there.

Neighbor-joining in Python: doing the plot

As discussed last time, I'm using Python lists to hold the information for each node on a phylogenetic tree. This seems as if it will have significant advantages, especially when it comes to tree "surgery." I'll leave that for another post. What I want to do today is actually draw a tree rooted at each of the internal nodes. Here is our example:

example_tree = { '0':['A','B','1'], '1':['C','0','3'], '2':['E','D','3'], '3':['F','1','2'] }

Code to traverse the tree was shown last time (here). I've modified it slightly to use a dictionary. Previously we had this representation for each node and its parent (in the traversal from a particular root):

A:0, C:1, B:0, E:2, D:2, F:3, 1:root, 0:1, 3:1, 2:3

Now these keys and values are in a dict. The results of the traversals from each internal node as root are:

A:0, C:1, B:0, E:2, D:2, F:3, 1:1, 0:1, 3:1, 2:3 A:0, C:1, B:0, E:2, D:2, F:3, 1:0, 0:0, 3:1, 2:3 A:0, C:1, B:0, E:2, D:2, F:3, 1:3, 0:1, 3:3, 2:3 A:0, C:1, B:0, E:2, D:2, F:3, 1:3, 0:1, 3:2, 2:2

The root node is now marked by having itself as its "parent."

The drawing code works by building a repr (representation) for each node from the farthest tips working up to the root. The method was inspired by PyCogent's ancestors() function, which returns a list of ancestors (naturally enough). The output with debug enabled looks like this:

ancestors E ['2', '3', '1', '0'] D ['2', '3', '1', '0'] F ['3', '1', '0'] 2 ['3', '1', '0'] C ['1', '0'] 3 ['1', '0'] A ['0'] B ['0'] 1 ['0'] 0 []

I simply sort on the length (largest first) and work in that order. So E's repr is just 'E'. When we process the first i_node (2), it's repr will be (E,D)2. If we're using branch lengths, it will be: (E:2.25,D:2.75)2. The distance data is in a separate dictionary:

dist = { ('A','0'):1.0, ('B','0'):4.0, ('C','1'):2.0, ('D','2'):2.75, ('E','2'):2.25, ('F','3'):4.75, ('0','1'):1.0, ('1','3'):1.25, ('2','3'):0.75 }

Here's the resulting Newick tree:

(A:1.00,B:4.00,(C:2.00,(F:4.75,(E:2.25,D:2.75)2:0.75)3:1.25)1:1.00)0;

One of the plots is shown at the top, a second one is below. Although I had it working before, be advised there's a bug somewhere that blows up R when I try to plot the tree as 'unrooted.' So now my to_do list includes the challenge of fixing that. [UPDATE: It seems to only happen when calling R from Python using RPy---probably related to the use of a Python keyword type as an argument to the R function. ] Zipped project files are here. And a final note, I've included the upgma code from before. However, we processed the nodes with a different algorithm there, so the drawing code is different, and I haven't tried to convert it to the new method yet.

Tree representation and one kind of traversal

Shown above is the tree we generated using the neighbor-joining method (here and here). I've annotated the plot from ape to label the internal nodes and show all the branch lengths. The internal node labels were assigned in the same order as they were clustered to build the tree. I didn't mention it before, but it's important to remember that this is an unrooted tree, and that's why it's drawn as it is, rather than with all the branches laid out horizontally, as in the UPGMA example (here).

I want to spend some time on trees in general, without worrying about how they were built. One goal is to understand the different modes of tree traversal (preorder and postorder and so on). Another is to explore different internal representations for trees in Python code, and conversions to and back from Newick format (and possibly XML format). In addition, I'd like to understand what it takes (in code) to "root" a tree, and how to implement various methods for tree "surgery" like "tree pruning and re-grafting" (e.g. here), and what the implications are for the Newick structure, which I find a bit hard to grasp.

In today's (very simple) example we're just going to try to visit all the nodes. The first thing is to decide on a way to represent the tree. The basic structure is naturally a dictionary, because we want to be able to "look up" the data for any individual node quickly. In a simplifying (and liberating) decision we're going to use a simple Python list to hold that data. One advantage of this is that we could accommodate star topologies easily. (And we could remember distance information by using a tuple of two lists, the names in one and distances in the second). Although a string would also do here, in general it won't work because we want to allow labels with len(label) > 1.

Only the internal nodes are in the dict; it looks like this:

tree_dict = { '0':['A','B','1'], '1':['C','0','3'], '2':['E','D','3'], '3':['F','1','2'] }

At this point, we have no notion of directionality. That's the liberating part. We will gain more structure as we traverse the tree.

Code to visit all the nodes naturally involves recursion, a function which processes nodes and calls itself for each of the child nodes, but returns if the current node is a tip (external). The first script (below) gives the following output, the node where we start the traversal is listed first:

['1', 'C', '0', 'A', 'B', '3', 'F', '2', 'E', 'D'] ['0', 'A', 'B', '1', 'C', '3', 'F', '2', 'E', 'D'] ['3', 'F', '1', 'C', '0', 'A', 'B', '2', 'E', 'D'] ['2', 'E', 'D', '3', 'F', '1', 'C', '0', 'A', 'B']

There's one simple modification we can make to allow us to actually draw the tree. That is to record the relationship between two nodes as we descend into the tree. The second script is modified to do that. It prints:

A:0, C:1, B:0, E:2, D:2, F:3, 1:root, 0:1, 3:1, 2:3 A:0, C:1, B:0, E:2, D:2, F:3, 1:0, 0:root, 3:1, 2:3 A:0, C:1, B:0, E:2, D:2, F:3, 1:3, 0:1, 3:root, 2:3 A:0, C:1, B:0, E:2, D:2, F:3, 1:3, 0:1, 3:2, 2:root

Here, each value is a pair where the first name is the node, and the second is its parent. The "root" node is identified by using that label in place of a parent. Now we will be able to draw the tree, but actually doing that will involve some complexities which are better dealt with in a separate post.

utils contains flatten from here.

code listings:

import utils tree_dict = { '0':['A','B','1'], '1':['C','0','3'], '2':['E','D','3'], '3':['F','1','2'] } all = list(utils.flatten(tree_dict.values())) external = [e for e in all if not e in tree_dict] def descend(node,seen): seen.append(node) if node in tree_dict: for child in tree_dict[node]: if not child in seen: descend(child,seen) def traverse_tree(root): seen = list() descend(i_node,seen) print seen for i_node in tree_dict: traverse_tree(i_node)

import utils tree_dict = { '0':['A','B','1'], '1':['C','0','3'], '2':['E','D','3'], '3':['F','1','2'] } all = list(utils.flatten(tree_dict.values())) external = [e for e in all if not e in tree_dict] def descend(node,caller,pD): pD[node] = caller if node in tree_dict: for child in tree_dict[node]: if not child in pD: descend(child,node,pD) def traverse_tree(root): pD = dict() descend(root,'root',pD) # might filter for k in tree_dict L = [k + ':' + pD[k] for k in pD] print ', '.join(L) for i_node in tree_dict: traverse_tree(root=i_node)

My version of neighbor-joining in Python

I finished debugging my version of the neighbor-joining method in Python. The zipped project files are here. I also made an nj tree from the same data using PyCogent, ape, and Phylip's neighbor. Everybody agrees!

One last thing I haven't done is to write the code to assemble the tree's representation. What I actually have is a list of internal nodes with distances to their child nodes. There's always something more to do.

PyCogent:

edge.0 A:1.0; edge.0 B:4.0; edge.1 C:2.0; edge.2 F:4.75; root D:2.75; root E:2.25; edge.2 edge.1 1.25 edge.2 edge.0 2.25 edge.1 edge.0 1.0 root edge.2 0.75 root edge.1 2.0 root edge.0 3.0 ((((A:1.0,B:4.0):1.0,C:2.0):1.25,F:4.75):0.75,D:2.75,E:2.25);

Phylip:

Between And Length ------- --- ------ 1 B 4.00000 1 2 1.00000 2 C 2.00000 2 4 1.25000 4 3 0.75000 3 E 2.25000 3 D 2.75000 4 F 4.75000 1 A 1.00000 (B:4.00000,(C:2.00000,((E:2.25000,D:2.75000):0.75000,F:4.75000):1.25000):1.00000,A:1.00000);

ape:

> tr$edge.length [1] 2.25 2.75 0.75 1.25 1.00 1.00 4.00 2.00 4.75 > tr$edge [,1] [,2] [1,] 7 5 [2,] 7 4 [3,] 7 8 [4,] 8 9 [5,] 9 10 [6,] 10 1 [7,] 10 2 [8,] 9 3 [9,] 8 6 > tr$tip.label [1] "A" "B" "C" "D" "E" "F" (E:2.25,D:2.75,(((A:1,B:4):1,C:2):1.25,F:4.75):0.75);

me:

0 : B 4.000 A 1.000 1 : C 2.000 0 1.000 2 : E 2.250 D 2.750 3 : 1 1.250 2 0.750 F 4.75

Checking our test data in NJ using Phylip

I'd like to compare the results of my own neighbor-joining script to standard methods such as Phylip, ape, and PyCogent. In initial comparisons, I noticed some small differences among those methods themselves with the data from before (UPDATE: I just failed to synch the data properly). (here). I also noticed that in Phylip's input screen for neighbor, they give the option of :

Outgroup root? No, use as outgroup species 1 Randomize input order of species? No. Use input order

I thought it would be interesting to see if there are any differences observed in the tree if we use a different input order of species. We could use Python to rewrite the data file to change to a specific order, but as a simpler test we'll just let Phylip randomize it for us. After getting the kinks out, we do 50 reps, calling neighbor as a process (just like here). One wrinkle: the process returns faster than the file write from Phylip, so we have to wait a second before trying to do anything with the files ourselves.

The data for Phylip look like this (note requirement for 10 characters on each line preceeding the distances, and a count for the otus on line 1):

6 A 0 5 4 7 6 8 B 5 0 7 10 9 11 C 4 7 0 7 6 8 D 7 10 7 0 5 8 E 6 9 6 5 0 8 F 8 11 8 8 8 0

As a bonus, I asked PyCogent to compare the topology of all the resulting trees, but no differences were found. And then, we plot the first one with ape.

import os, subprocess, time, sys, random from cogent import LoadTree prefix = '/Users/telliott_admin/Desktop/' controlfn = prefix + 'responses.txt' datafn = prefix + 'nj_data.phylip.txt' resultfn = prefix + 'outfile' treefn = prefix + 'outtree' def one_run(n): n = str(n) FH = open(controlfn, 'w') FH.write(datafn + '\n') FH.write('J\n') seed = random.choice(range(2,200)) if not seed % 2: seed += 1 FH.write(str(seed) + '\n') FH.write('Y\n') FH.close() cmd = 'neighbor < ' + controlfn + ' > ' cmd += 'screenout &' p = subprocess.Popen(cmd, shell=True) pid,ecode = os.waitpid(p.pid, 0) print ecode time.sleep(1) os.rename(resultfn,prefix+'results/outfile' + n + '.txt') os.rename(treefn,prefix+'results/outtree' + n + '.txt') os.remove(prefix + 'screenout') N = 50 for i in range(1,N+1): one_run(i) L = list() for i in range(1,N+1): tr = LoadTree(prefix + 'results/outtree' + str(i) + '.txt') L.append(tr) for i in range(N): for j in range(N): if i == j: continue if not L[i].sameTopology(L[j]): print L[i] print L[j] else: print '.', ''' R code: library(ape) setwd('Desktop') tr = read.tree('results/outtree1.txt') plot(tr,edge.width=3,cex=2,type='unrooted') axisPhylo() '''

UPGMA in Python 3: Sarich data

In Joe Felsenstein's book, Inferring Phylogenies, he uses some data from Vincent Sarich to demonstrate UPGMA. Although I was tempted to "declare victory and go home", I decided to test my version of UPGMA with that data. And that turned up bugs in both the main upgma and the plotting code! So I guess the moral is: test, test, and test again.

The top figure is a scan from the book. Here is what my program plotted:



It looks pretty good to me.

The data are immunological data from eight species in this order: dog, bear, raccoon, weasel, seal, sea lion, cat, monkey. The data are in this distance matrix:

0 32 48 51 50 48 98 148 32 0 26 34 29 33 84 136 48 26 0 42 44 44 92 152 51 34 42 0 44 38 86 142 50 29 44 44 0 24 89 142 48 33 44 38 24 0 90 142 98 84 92 86 89 90 0 148 148 136 152 142 142 142 148 0

If you grab the zipped files (here) and run them, you'll see a lot of diagnostic output when debug == True. As well as all the branch lengths.

UPGMA in Python 2

Following up on the previous post (here, and an earlier post here), I wrote some code to extract the data for our UPGMA tree, and plot it using the ape package from R. I updated the zipped project files which are on Dropbox (here). The plot is above, it looks good to me.

To do: test on other examples! [Update: found a bug and updated again.]

Neighbor-joining

Neighbor-joining is a standard way of constructing phylogenetic trees that is fast and does a reasonable job of recovering the correct tree in the case of unequal branch lengths. Consider this tree with vastly unequal amounts of change on the branch to B. I found this example (and the figure above) on the internet, but now I've forgotten where (sorry).

Higgs and Attwood have a description of the method that explains a bit about why the calculation is the way it is, and what guarantees it offers.

Here's the distance data:

A B C D E B 5 C 4 7 D 7 10 7 E 6 9 6 5 F 8 11 8 8 8

Step 1 is to calculate the net divergence for each OTU:

r(A) = 30 = 5 + 4 + 7 + 6 + 8 r(B) = 42 etc. r(C) = 32 r(D) = 37 r(E) = 34 r(F) = 43

Step 2 is to modify the distance matrix by subtracting the value of: (ri + rj)/(N-2) where N is the number of OTUs. For example, the new value for entry 'AB' is: 5 - ((30 + 42)/4) = -13, and that is the smallest value in the whole matrix.

We choose the next two otus to cluster based on the smallest value in the modified distance matrix. In fact, we don't use this modified matrix again, so I just tested each value in turn against the low value observed so far, and kept the relevant indexes. We make a "star" tree, and then group the 2 OTUs we chose:

Step 3. Calculate the distances of A and B from U using the data from the original distance matrix and the divergence values, as follows:

S(AU) = d(AB)/2 + [r(A) - r(B)/(2(N-2))] = 5/2 + (-12 / 8) = 1 S(BU) = d(AB) - S(AU) = 4

Step 4. Now calculate the distances from the new node U to the other nodes:

d(UC) = ( d(AC) + d(BC) - d(AB) ) / 2 = 3 d(UD) = ( d(AD) + d(BD) - d(AB) ) / 2 = 6 d(UE) = ( d(AE) + d(BE) - d(AB) ) / 2 = 5 d(UF) = ( d(AF) + d(BF) - d(AB) ) / 2 = 7

The result is a new distance matrix:

U C D E C 3 D 6 7 E 5 6 5 F 7 8 8 8

Notice that the result from this first step already looks a little more like the real tree. We repeat each of the steps until we've collapsed all the nodes. The result is an unrooted tree.



R code:

m = c(0,5,4,7,6,8, 5,0,7,10,9,11, 4,7,0,7,6,8, 7,10,7,0,5,9, 6,9,6,5,0,8, 8,11,8,8,8,0) dim(m)=c(6,6) colnames(m)=c( 'A','B','C','D','E','F') rownames(m)=colnames(m) d = as.dist(m) library(ape) tr=nj(d) plot(tr,cex=2, tip.color=rainbow(6)) axisPhylo()

R chose a root to draw the tree that is different than in the other graphics. I forgot to fix this (with ape.root).



[UPDATE: I'm writing a Python version of neighbor-joining, and it turned up some errors in my arithmetic. So I updated the table here.

UPGMA in Python

I spent a whole day working on a script to do UPGMA. (It took a lot longer than I thought it should). This first version analyzes the data from the same tree as we constructed in an earlier post (here), because it's simple.

At the end of the run, we have the correct tree, as shown by the first line in the last section of the output:

(((('A', 'B'), 'C'), ('D', 'E')), 'F') A {'up': 1.0, 'parent': 'AB'} C {'up': 2.0, 'parent': 'ABC'} B {'up': 1.0, 'parent': 'AB'} E {'up': 2.0, 'parent': 'DE'} D {'up': 2.0, 'parent': 'DE'} F {'up': 4.0, 'parent': 'ABCDEF'} DE {'to_tips': 2.0, 'right': 'E', 'up': 1.0, 'parent': 'ABCDE', 'left': 'D'} ABCDE {'to_tips': 3.0, 'right': 'DE', 'up': 1.0, 'parent': 'ABCDEF', 'left': 'ABC'} ABC {'to_tips': 2.0, 'right': 'C', 'up': 1.0, 'parent': 'ABCDE', 'left': 'AB'} AB {'to_tips': 1.0, 'right': 'B', 'up': 1.0, 'parent': 'ABC', 'left': 'A'} ABCDEF {'to_tips': 4.0, 'right': 'F', 'left': 'ABCDE'}

I struggled with trying to use the distance information to modify the string representation of the tuple on that first line, but failed. I know I can reconstruct the tree from the data for all the nodes that is in the dictionaries shown. That's for another time.

Once again, I have a better understanding through having coded an (almost) working example in Python. The files are on Dropbox (here).  [Edit:  the files are here].

Likelihood of an evolutionary tree 2

As usual, to improve my understanding of the likelihood calculation, I wrote some code. The first module is called likelihood.py, and it defines three functions for processing nodes, depending on the nature of their child nodes (both external, both internal, or mixed).

It's highly commented and should be self-explanatory after the previous post (here). To keep things simple, I have not implemented variable branch lengths yet.

The problem I'm having is testing. The version that I'm actually developing has lots of statements like: if debug: print .., where I've been rigid about putting the print on an indented line below the if. (I've stripped all of those out in the listing below). I went through the printout trying to verify the results of each calculation manually, and it's pretty exhausting. What would be nice is to have another ML program to test the output against, but without implementing the branch lengths properly, I'm not sure I can do that.

This module itself is sort of boring. I wrote another module that exercises it a lot more, but I'm not quite ready to show that yet.

The code to strip out debug statements is short and sweet:

data = tu.load_data(fn) data = data.strip().split('\n') def lwspace(line): return len(line) - len(line.lstrip()) flag = False target = 'debug:' indent = 0 rL = list() for line in data: if 'debug:' in line: assert 'if' in line flag = True indent = lwspace(line) print '*' + line continue if lwspace(line) <= indent: flag = False if flag: print '*' + line else: rL.append(line) print ' ' + line

likelihood.py

import math # list printing def printL(L,name=None): if name: print name pL = [str(round(n,3)).rjust(8) for n in L] N = 6 while pL: line = pL[:N] pL = pL[N:] print ''.join(line) # to start with, we'll assume equal branch lengths # construct the transition-probability matrix nt = 'ACGT' pur = 'AG' pyr = 'CT' f_same = math.log(0.95) f_transition = math.log(0.03) f_transversion = math.log(0.01) pi = math.log(0.25) def get_f(n): m,n = list(n) if m == n: return f_same if m in pur and n in pur: return f_transition if m in pyr and n in pyr: return f_transition return f_transversion k = [m+n for m in nt for n in nt] v = [get_f(n) for n in k] P = dict(zip(k,v)) #--------------------------------------- # three different functions # depending on types of child nodes def ext_ext(u,v): # u,v are nucleotide as char # for each possible n in their parent node # returns an array of log(p) for nu and nv rL = [P[n+u] + P[n+v] for n in nt] return rL # L = list of likelihoods in order, for internal child # v = external child nucleotide as char def int_ext(L,v): rL = list() # for new node in {ACGT} for i,m in enumerate(nt): ep = P[m+v] # log(p) of m -> v sL = list() # each state for internal child # n is a float for j,n in enumerate(L): u = nt[j] # would multiply probs, so add logs p = P[m+u] # log(p) for next branch p += L[j] # log likelihood if that child = u p += ep # log(p) for external child sL.append(p) # will need the actual probs to add them # does this need to be done better? how? sL = [(math.e)**p for p in sL] logS = math.log(sum(sL)) rL.append(logS) return rL # both children are internal def int_int(L1,L2,root=False): rL = list() # for new node = {ACGT} for i,n in enumerate(nt): sL = list() # each state for left child # v1 is a float for j,f1 in enumerate(L1): u = nt[j] # each state for right child # v2 is a float for k,f2 in enumerate(L2): v = nt[k] p = P[n+u] # log(p) for left branch p += f1 # log likelihood if child = u p += P[n+v] # log(p) for right branch p += f2 # log likelihood if child = v sL.append(p) # will need the actual probs to add them sL = [(math.e)**p for p in sL] logS = math.log(sum(sL)) rL.append(math.log(sum(sL))) if root: # do pi calculation rL = [e + pi for e in rL] return rL if __name__ == '__main__': #for k in sorted(P.keys()): print k, P[k] #for u in nt: #for v in nt: #ext_ext(u,v) #print '-'*40 L = [-0.1] * 4 int_ext(L,'C') print '-'*40 int_int(L,L)

Likelihood of an evolutionary tree

Calculating the likelihood of a particular phylogenetic tree sounds complicated and it's easy to get lost in the details, but is really just a matter of bookkeeping. That's why we have computers! This example is taken from my all-time favorite Bioinformatics textbook (Higgs & Attwood).

Let's work with the tree shown above. We're considering one particular site in a sequence alignment from a set of OTUs (operational taxonomic units), where the observed nucleotides are as shown. The identities of the internal nodes are unknown and each could be any one of the four nucleotides. (I've shown the same tree several times at different places in the post for reference).

Consider first the node labeled X. In general, the branches from X to A and from X to G may be different in length. Suppose the first one is described by t₁ and the second by t₂.

We have a particular model of sequence evolution for this site, encapsulated in a transition-probability matrix (as discussed here), which allows us to calculate the probabilities for all 4 possibilities with no change:

PAA, PCC..

and all 12 possible changes:

PAC, PAG, PAT..

The identity of X isn't known, it might be any of {A,C,G,T}. Starting with A, then we can calculate the probabilities for the two branches as:

PAA(t1) PAG(t2)

The probability of the observed data for both of these two terminal nodes with the model that X = A, is the likelihood that X = A, given the data. Namely:

L(X=A) = PAA(t1) PAG(t2)

In the same way:

L(X=C) = PCA(t1) PCG(t2) L(X=G) = PGA(t1) PGG(t2) L(X=T) = PTA(t1) PTG(t2)

Of course, the total probability for X = {A,C,G,T} is the sum:

Σ PNA(t1) PNG(t2) N = {A,C,G,T}

which has to equal 1, since X must be one of the four nucleotides.

Repeated multiplication leads to various difficulties which are avoided by using logarithms. So typically we would talk about the log likelihood, for example the log of L(X=C). All of the calculations shown below as multiplications would be done in practice by adding the logarithm of each term.

Working up the tree

Now consider node Y and its branch to X, with a branch length described by t₃.

We can calculate the probabilty for any choice of nucleotides (say, Y = A, X = C) in the same way as before. For this choice it is:

PAC(t3) PCA(t1) PCG(t2)

It can also be written using the expression for likelihood L(X=C) as:

PAC(t3) L(X=C)

Stick with probabilities for a moment.

There are 16 possible paths here. Suppose Y = M and X = N, where both M and N are in the set {A,C,G,T}. Then there are 16 terms like:

PMN(t3) PNA(t1) PNG(t2)

This is only for the left-hand clade. We usually talk about the tree as if it were rotated with Z at the top and Y on the left, but it's easier to draw as shown. The probability for the right-hand branch from Y to G described by t₄, is calculated for any particular choice of M = {A,C,G,T} as described in the first part:

PMG(t4)

So each of our 16 terms picks up an additional factor

PMG(t4) PMN(t3) PNA(t1) PNG(t2)

and the sum of the four terms involving Y=M (for a particular choice of M = A,C,G or T) is the likelihood L(Y=M).

A node with two internal child nodes

Now consider node W at the top of the tree and the branch t₅ connecting it to Y. We let W assume any of the four states {A,C,G,T} and calculate probabilities as follows:

P(W=A,Y=A) = PAA(t5) P(W=A,Y=C) = PAC(t5) ..

For the t₅ branch, we have three equivalent expressions:

P(W=A,Y=A) L(Y=A) P(W=A,Y=A) P(Y=A,X=A) L(X=A) P(W=A,Y=A) P(Y=A,X=A) P(X=A,left=A) P(X=A,right=G)

There are 16 terms like:

P(W=A,Y=A) L(Y=A)

for each of the 16 possible values for the tuple (W,Y). Once again, this is only for the left-hand clade. Each of the 16 terms picks up another factor from the right-hand branch:

P(W=A,Y=A) L(Y=A) P(W=A,Z=A) L(Z=A)

So the total likelihood L(W=A,Y={A,C,G,T},Z={A,C,G,T}) is the sum of:

P(W=A,Y=A) L(Y=A) P(W=A,Z=A) L(Z=A) P(W=A,Y=C) L(Y=C) P(W=A,Z=A) L(Z=A) P(W=A,Y=G) L(Y=G) P(W=A,Z=A) L(Z=A) P(W=A,Y=T) L(Y=T) P(W=A,Z=A) L(Z=A) P(W=A,Y=A) L(Y=A) P(W=A,Z=C) L(Z=C) P(W=A,Y=C) L(Y=C) P(W=A,Z=C) L(Z=C) P(W=A,Y=G) L(Y=G) P(W=A,Z=C) L(Z=C) P(W=A,Y=T) L(Y=T) P(W=A,Z=C) L(Z=C) P(W=A,Y=A) L(Y=A) P(W=A,Z=G) L(Z=G) P(W=A,Y=C) L(Y=C) P(W=A,Z=G) L(Z=G) P(W=A,Y=G) L(Y=G) P(W=A,Z=G) L(Z=G) P(W=A,Y=T) L(Y=T) P(W=A,Z=G) L(Z=G) P(W=A,Y=A) L(Y=A) P(W=A,Z=T) L(Z=T) P(W=A,Y=C) L(Y=C) P(W=A,Z=T) L(Z=T) P(W=A,Y=G) L(Y=G) P(W=A,Z=T) L(Z=T) P(W=A,Y=T) L(Y=T) P(W=A,Z=T) L(Z=T)

An equivalent description is to call this a sum of sums, or double sum:

Σ Σ P(W=A,Y=M) L(Y=M) P(W=A,Z=N) L(Z=N) M = {A,C,G,T}, N = {A,C,G,T}

The expression is symmetrical and you could use either M or N in the outside sum. But the double sum makes my head spin a bit. So I'll stick with the longer (and perhaps less confusing) formula.

We can abbreviate this as L(W=A).

Top of the tree

At the top of the tree, there is a final factor which comes in. It is π, the equilibrium frequency of each nucleotide according to our evolutionary model. The total likelihood for the tree is the sum over the four nucleotides:

Σ πN L(W=N) N={A,C,G,T}

An alignment of many sites

As Higgs and Attwood say:

The likelihood is calculated for each site in this way, and it is assumed that sites evolve independently. Therefore the likelihood for the complete sequence set, L_tot, is the product of the likelihoods of the sites. It therefore follows that the log of the likelihood is the sum of the logs of the site likelihoods.

log Ltot = Σ {sites} log Lsite

Like I said, it's just a matter of bookkeeping.

RPy and ape working together

I couldn't resist doing a little more with RPy and the R ape package. I used our standard tree file. The code listing is below. I am quite pleased with the results. With this in hand, it will be nice to look more at the ape book, which I already have.

In the first section we do our imports, and get the R object 'T' (True), because I couldn't get True or 'True' to work as a direct argument to the plot function. We also read the tree file. Next, I demonstrate the form of the ape function drop_tip, which is very useful, but I don't show results or a plot here.

In the third part, I define a function to return a list of colors based on the names of the external nodes, and then later convert that to an appropriate R object. This is what had me so frustrated when I first tried ape.

In the last part, we do the plot to a pdf file. The x_lim argument is important, that's how we stretch out the plot along the x-axis to make enough room for the names. Two separate calls after plot show the internal node labels, and the x-axis values.

I agree this is really powerful. The best part is that for someone like me, when you have trouble doing something in R, you can drop into Python without a problem. It's easy when you know how!

One more thing, Laurent Gautier is big in RPy, and one of his interests is Bioconductor. That project is the reason I got into R in the first place. It'll be fun to explore that in the future.

from rpy2 import robjects from rpy2.robjects.packages import importr grdevices = importr('grDevices') T = robjects.r['T'] ape = importr('ape') fn = '/Users/telliott_admin/Desktop/tree.txt' tr = ape.read_tree(fn) print tr #-------------------------------------------- # the ape function drop_tip tip_list = robjects.IntVector([1,6]) #tr = ape.drop_tip(tr,tip_list) #-------------------------------------------- # define some colors: # this works but isn't very flexible # col_list = robjects.r("rep(c('red','blue'),3)") def get_color(n): D = { 'Steno':'maroon','Kingella':'maroon', 'Haemo':'DodgerBlue', 'coli':'magenta','typhi':'magenta' } for k in D: if k in n: return D[k] return 'darkgreen' tip_labels = tr[2] color_list = [get_color(n) for n in tip_labels] color_list = robjects.StrVector(color_list) #-------------------------------------------- ofn = '/Users/telliott_admin/Desktop/plot.pdf' grdevices.pdf(ofn) ape.plot_phylo(tr,tip_color=color_list, cex=1.5,x_lim=0.23) ape.nodelabels(cex=2) ape.axisPhylo() grdevices.dev_off()

Phylo_plotter

I revised this project a bit and replaced the zipped files on Dropbox (here). I got the fonts looking pretty good. And I modified it to use a config file that specifies values for various attributes. The attr dictionary can be modified by including another file that just shows the changed values. That should make it easy enough for a non-programmer to use (I think). It looks like this:

e_node_label_size:10 e_node_dots_visible:True e_node_bar_color:k e_node_dot_size:25

I'm still thinking about how to use forms to run this thing. But the first answer I got from Stack Overflow is that you have to use a server. Maybe I should just write an app that acts like a browser but forwards to a local script?

Notes on some issues to solve with RPy

I upgraded to R version 2.12.0 the other day. The R GUI shows on my 32-bit iMac at work (5 years old):

R version 2.12.0 (2010-10-15) [R.app GUI 1.35 (5632) i386-apple-darwin9.8.0]

$ python Python 2.6.1 (r261:67515, Jun 24 2010, 21:47:49) [GCC 4.2.1 (Apple Inc. build 5646)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> from rpy2 import robjects Bus error

Oops..
Luckily, CRAN still stockpiles old binary installers. Grabbed 2.11.1 here, looks like it installed correctly. The R GUI shows:

R version 2.11.1 (2010-05-31)

However the rpy2 install is messed up. Removed the egg from site-packages, reinstall rpy2 using easy_install. And it works. Need to re-install the ape package in this version of R.

We'll use the tree from the other day (here). Typical R code for this plot:

library(ape) setwd('Desktop') tree = read.tree('tree.txt') plot(tree) axisPhylo()

And it works! (You'll just have to trust me).
Now, how to do it using RPy?

>>> read_tree = robjects.r['read_tree'] Traceback (most recent call last): File "", line 1, in File "build/bdist.macosx-10.6-universal/egg/rpy2/robjects/__init__.py", line 241, in __getitem__ LookupError: 'read_tree' not found

read.dna is not found either..

I'm going to put this up for reference and see if I can get some help. I'll be back later.

And then while re-testing some of these examples I got this:

$ python Python 2.6.1 (r261:67515, Jun 24 2010, 21:47:49) [GCC 4.2.1 (Apple Inc. build 5646)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> from rpy2 import robjects Bus error

Hunhh? How could it stop working like that?

New plotter for phylogenetic trees: version 0.1

This is the last post on this project. The previous posts should be listed in the archive in the sidebar.

I have a phylogenetics project for which I constructed a boutique database. The definition is in the file db.groups.txt, one section of which looks like this:

aeromonas
X60415.1 Aeromonas_trota_ATCC49657
X74677.1 Aeromonas_hydrophila_ATCC7966T
EU770300.1 Aeromonas_enteropelogenes_MS12

set off from surrounding entries by double newlines. The "group" aeromonas has three sequences. There are also super-groups that correspond to bacterial Phyla or sub-Phyla, e.g., the gamma-Proteobacteria. These are defined in the file big.groups.txt like this:

gamma
aeromonas
cardio
pseudo
moraxella
entero
haemo1
haemo2
steno
xantho

Normally the sequences would be fetched from Genbank by a slightly complicated script, but I'd like you to be able to follow along if you want, so we'll do it almost manually. The file fetchSeqs.py contains a function that parses these two files and grabs the info we need for the "gamma" super-group. It also does a manual fetch from Genbank. We save the data to disk like so:

python fetchSeqs.py > seqs.txt

The beginning of one of the entries looks like this:

>gi|312967|emb|X60415.1| Aeromonas trota 16S rRNA gene, strain ATCC 49657 GAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAGCGGCAGCGGGAA

As you can see, the title line of this FASTA-formatted text is really long. So in the second script, analyzeSeqs.py, we use the same info from before to replace the long title with something shorter. This script also uses MUSCLE to align the sequences. You need to have it installed for this to work. Or you could use Clustal. I used the ape library in R to make a neighbor-joining tree and write this to disk as 'tree.txt'. I show a plot of what R gave me (actually the real process for this one was more complicated, resulting the pretty colors). But at least you can see what the tree is supposed to look like.

The last module is test_plotter.py It shows how to make a couple different kinds of plots with our plotter. All of them have colored node labels as defined in the script. The first plot is the standard one, the graphic is at the top of the post. In the next, we plot with the internal node labels showing, so that we can identify the name of a node to re-root the tree, if wanted.

In the third example, I got a balanced tree from PyCogent and plot that.

You might notice that the font used for the external node labels is now italic as it should be. It doesn't look that hot, I think that matplotlib is not using the OS X system fonts. Not sure why that is, yet.

I put a zip of all the files up on Dropbox (here).

New plotter for phylogenetic trees: customization

Continuing with a tree-plotter in Python. Previous posts here, here, and here. These examples depend on PyCogent and matplotlib.

I've got a somewhat bigger tree here, using sequences from a project we're working on. I used clustal and R to make a neighbor-joining tree. As you can see in the code, there is a lot of flexibility in terms of customized colors, etc.

code listing:

import tree_utils as TU import tree_plotter as TP fn = 'veillonella.tree.txt' D = TU.make_tree_dict(fn=fn) tr = D['meta']['tree'] attr = TP.get_default_attr() attr['e_node_dot_size'] = 35 attr['figure_type'] = 'png' attr['e_node_label_size'] = 10 attr['e_node_bar_color'] = 'k' attr['using_alternate_names'] = True e_node_names = D['meta']['e_node_names'] #L = [e.replace('_',' ') for e in e_node_names] L = [e.split('_')[-1] for e in e_node_names] # lose the version numbers on Genbank seqs L = [e.split('.')[0] for e in e_node_names] alt_name_dict = dict(zip(e_node_names,L)) attr['alternate_name_dict'] = alt_name_dict attr['using_e_node_specific_colors'] = True def f(name): if 'DB' in name: return 'b' return 'r' L = [f(e) for e in e_node_names] node_label_color_dict = dict(zip(e_node_names,L)) attr['node_label_color_dict'] = node_label_color_dict def f(name): if 'DB' in name: return 'b' return 'r' L = [f(e) for e in e_node_names] dot_color_dict = dict(zip(e_node_names,L)) attr['dot_color_dict'] = dot_color_dict TP.plot(D,attr=attr) ''' Rcode: library(ape) setwd('Desktop') dna = read.dna('veillonella.fasta',format='fasta') tree = nj(dist.dna(dna)) plot(tree,cex=1) axisPhylo() write.tree(tree,'veillonella.tree.txt') '''

New plotter for phylogenetic trees: re-rooting

Continuing with a tree-plotter in Python. Previous posts here and here. These examples depend on PyCogent and matplotlib.

One of the great things about building on PyCogent is we can use its re-rooting routine. In the output for this example we print the same tree as before, rooted at all possible internal nodes.

I've made a few tweaks to the code. make_tree_dict will take either a tree_string, a filename or a PyCogent tree. For the graphic above, we grab the default plot attributes, make the internal nodes visible, and then call plot. The code for this example is at the bottom. (I'm not going to post modifications to the other code until it's finalized. If you want it, shoot me an email).

To do list:
• learn matplotlib typography (e.g. italics)
• adjust for label sizes in a smart way
• provide for printing bootstrap values or internal node names

output:

$ python test_plotter.py original root /-Stenotrophomonas_maltophilia /edge.0--| | \-Kingella_oralis | -root----|--Pseudomonas_aeruginosa | | /-Salmonella_typhi | /edge.1--| \edge.2--| \-Escherichia_coli | \-Haemophilus_parainfluenzae ------------------------------------------------------------ edge.0 /-Stenotrophomonas_maltophilia | |--Kingella_oralis -root----| | /-Pseudomonas_aeruginosa | | \edge.0--| /-Salmonella_typhi | /edge.1--| \edge.2--| \-Escherichia_coli | \-Haemophilus_parainfluenzae ------------------------------------------------------------ edge.2 /-Salmonella_typhi /edge.1--| | \-Escherichia_coli | -root----|--Haemophilus_parainfluenzae | | /-Stenotrophomonas_maltophilia | /edge.0--| \edge.2--| \-Kingella_oralis | \-Pseudomonas_aeruginosa ------------------------------------------------------------ edge.1 /-Salmonella_typhi | |--Escherichia_coli -root----| | /-Haemophilus_parainfluenzae | | \edge.1--| /-Stenotrophomonas_maltophilia | /edge.0--| \edge.2--| \-Kingella_oralis | \-Pseudomonas_aeruginosa ------------------------------------------------------------ balanced tree /-Stenotrophomonas_maltophilia /edge.0--| | \-Kingella_oralis | -root----|--Pseudomonas_aeruginosa | | /-Salmonella_typhi | /edge.1--| \edge.2--| \-Escherichia_coli | \-Haemophilus_parainfluenzae

code listing:

import tree_utils as TU import tree_plotter as TP fn = 'tree.txt' D = TU.make_tree_dict(fn=fn) tr = D['meta']['tree'] print 'original root' print tr.asciiArt() print '-'*60 L = D['meta']['i_node_names'] for name in L[1:]: print name print tr.rootedAt(name).asciiArt() D = TU.make_tree_dict(tr=tr) print '-'*60 print 'balanced tree' print tr.balanced().asciiArt() tr = tr.rootedAt('edge.2') D = TU.make_tree_dict(tr=tr) attr = TP.get_default_attr() attr['i_node_dots_visible'] = True TP.plot(D,attr=attr,ofn='balanced')

New plotter for phylogenetic trees: plotting

Here is a first pass at a plotter, following up on the last post.

It's ugly code yet, but I thought I would show you what I have. What I need to do is really think harder about how this module should be organized. But you can compare the result with what we got here.

[UPDATE: Modified the code to be easily customizable. Next up: more extensive testing.]
UPDATE2: Modified the code yet again.]

output:

alternate name dict None default e node dot color k dot color dict None e node bar color r e node bar color list None e node dot size 75 e node dots visible True e node label default color k e node label size 14 e node labels visible True figure type png horizontal axis visible True i node bar color k i node dot color magenta i node dots visible False i node vertical bar color k label width factor 1.1 line width 2 node label color dict None r node dot color orange using alternate names False using e node specific colors False vertical axis visible False

code listing:

import matplotlib.pyplot as plt import tree_utils as tu def plot(D,attr=None,ofn=None): if not attr: attr = get_default_attr() SZ = attr['e_node_dot_size'] lw = attr['line_width'] maxx, maxy = D['meta']['max_xy'] all_node_names = D['meta']['all_node_names'] e_node_names = D['meta']['e_node_names'] i_node_names = D['meta']['i_node_names'] # actual node dictionaries e_node_dicts = [D[n] for n in e_node_names] i_node_dicts = [D[n] for n in i_node_names] # extract the values we need e_node_positions = [(nD['x'],nD['y']) for nD in e_node_dicts] e_node_lengths = [(nD['dist_to_parent']) for nD in e_node_dicts] i_node_positions = [(nD['x'],nD['y']) for nD in i_node_dicts] i_node_lengths = [(nD['dist_to_parent']) for nD in i_node_dicts] i_node_verticals = [(nD['y_bott'], nD['y_top']) for nD in i_node_dicts] # external nodes for i in range(len(e_node_names)): x,y = e_node_positions[i] d = e_node_lengths[i] xleft = x - d # bars L = attr['e_node_bar_color_list'] if L: c = L[i] else: c = attr['e_node_bar_color'] plt.plot([xleft,x],[y,y],color=c,lw=lw,zorder=1) # dots name = e_node_names[i] if attr['e_node_dots_visible']: if attr['using_e_node_specific_colors']: c = attr['dot_color_dict'][name] else: c = attr['default_e_node_dot_color'] plt.scatter(x,y,color=c,s=SZ,zorder=2) # internal nodes, root is i = 0 for i in range(len(i_node_names))[1:]: x,y = i_node_positions[i] d = i_node_lengths[i] y_bott, y_top = i_node_verticals[i] x0 = x - d # bars c = attr['i_node_bar_color'] plt.plot([x0,x],[y,y],color=c,lw=lw,zorder=1) # verticals c = attr['i_node_vertical_bar_color'] plt.plot([x,x],[y_bott, y_top],color=c,lw=lw,zorder=1) # dots if attr['i_node_dots_visible']: c = attr['i_node_dot_color'] plt.scatter(x,y,color=c,s=SZ,zorder=2) # root i = 0 x,y = i_node_positions[i] y_bott, y_top = i_node_verticals[i] # verticals c = attr['i_node_vertical_bar_color'] plt.plot([x,x],[y_bott, y_top],color=c,lw=lw,zorder=1) # dots if attr['i_node_dots_visible']: c = attr['r_node_dot_color'] plt.scatter(x,y,color=c,s=SZ,zorder=2) L = [len(name) for name in e_node_names] max_label_width = maxx/100.0 * max(L) max_label_width *= attr['label_width_factor'] # external node labels if attr['e_node_labels_visible']: for i in range(len(e_node_names)): x,y = e_node_positions[i] dx = maxx/100.0 * 4 name = e_node_names[i] if attr['using_alternate_names']: s = attr['alternate_name_dict'][name] else: s = name if attr['using_e_node_specific_colors']: c = attr['node_label_color_dict'][name] else: c = attr['e_node_label_default_color'] plt.text(x + dx, y, s, fontname = 'Helvetica', fontsize = attr['e_node_label_size'], color = c, ha = 'left',va = 'center') ax = plt.axes() if not attr['vertical_axis_visible']: ax.yaxis.set_visible(False) if not attr['horizontal_axis_visible']: ax.xaxis.set_visible(False) ax.set_xlim(-maxx/10.0,maxx*1.1 + max_label_width) ax.set_ylim(-maxy/10.0,maxy*1.1) if ofn: plt.savefig(ofn + '.' + attr['figure_type']) else: plt.savefig('example.' + attr['figure_type']) def get_default_attr(): attr = dict() attr['e_node_dot_size'] = 75 attr['line_width'] = 2 attr['figure_type'] = 'png' attr['horizontal_axis_visible'] = True attr['vertical_axis_visible'] = False attr['e_node_bar_color'] = 'r' attr['e_node_bar_color_list'] = None attr['e_node_dots_visible'] = True attr['e_node_labels_visible'] = True attr['default_e_node_dot_color'] = 'k' attr['e_node_label_size'] = 14 attr['e_node_label_default_color'] = 'k' attr['label_width_factor'] = 1.1 attr['i_node_bar_color'] = 'k' attr['i_node_vertical_bar_color'] = 'k' attr['i_node_dots_visible'] = False attr['i_node_dot_color'] = 'magenta' attr['r_node_dot_color'] = 'orange' attr['using_alternate_names'] = False attr['alternate_name_dict'] = None attr['using_e_node_specific_colors'] = False attr['node_label_color_dict'] = None attr['dot_color_dict'] = None return attr def print_defaults(): attr = get_default_attr() N = max([len(k) for k in attr.keys()]) for k in sorted(attr.keys()): v = attr[k] k = k.replace('_',' ') print k.rjust(N), ' ', v if __name__ == '__main__': fn = 'tree.txt' file_data = tu.load_data(fn) D = tu.make_tree_dict(ts=file_data) plot(D) print_defaults()