【3.5.2】基于bioptyhon绘制进化树
一、初识tree的结构
创建测试文本 simple.dnd
cat simple.dnd
(((A,B),(C,D)),(E,F,G));
代码：
from Bio import Phylo
tree = Phylo.read("simple.dnd", "newick")
print(tree)
结果：
Tree(rooted=False, weight=1.0)
Clade()
Clade()
Clade()
Clade(name='A')
Clade(name='B')
Clade()
Clade(name='C')
Clade(name='D')
Clade()
Clade(name='E')
Clade(name='F')
Clade(name='G')
Tree 对象包含树的全局信息，如树是有根树还是无根树。它包含一个根进化枝， 和以此往下以列表嵌套的所有进化枝，直至叶子分支。
二、绘制tree
2.1 基本绘制
函数 draw_ascii 创建一个简单的ASCII-art(纯文本)系统发生图。在没有更好 图形工具的情况下，这对于交互研究来说是一个方便的可视化展示方式。
Phylo.draw_ascii(tree)
结果：
________________________ A
________________________|
| |________________________ B
________________________|
| | ________________________ C
| |________________________|
_| |________________________ D
|
| ________________________ E
| |
|________________________|________________________ F
|
|________________________ G
如果你安装有 matplotlib 或者 pylab, 你可以使用 draw 函数一个图像
tree.rooted = True
Phylo.draw(tree)
2.2 给树的分支上颜色
函数 draw 和 draw_graphviz 支持在树中显示不同的颜色和分支宽度。 从Biopython 1.59开始，Clade对象就开始支持 color 和 width 属性， 且使用他们不需要额外支持。这两个属性都表示导向给定的进化枝前面的分支的 属性，并依次往下作用，所以所有的后代分支在显示时也都继承相同的宽度和颜 色。
在早期的Biopython版本中，PhyloXML树有些特殊的特性，使用这些属性需要首先 将这个树转换为一个基本树对象的子类Phylogeny，该类在Bio.Phylo.PhyloXML模 块中。
在Biopython 1.55和之后的版本中，这是一个很方便的树方法：
tree = tree.as_phyloxml()
在Biopython 1.54中, 你能通过导入一个额外的模块实现相同的事情：
from Bio.Phylo.PhyloXML import Phylogeny
tree = Phylogeny.from_tree(tree)
注意Newick和Nexus文件类型并不支持分支颜色和宽度，如果你在Bio.Phylo中使用 这些属性，你只能保存这些值到PhyloXML格式中。（你也可以保存成Newick或Nexus 格式，但是颜色和宽度信息在输出的文件时会被忽略掉。）
现在我们开始指定颜色。首先，我们将设置根进化枝为灰色。我们能通过赋值24位 的颜色值来实现，用三位数的RGB值、HTML格式的十六进制字符串、或者预先设置好的 颜色名称。
tree.root.color = (128, 128, 128)
或者：
tree.root.color = "#808080"
或者：
tree.root.color = "gray"
一个进化枝的颜色会被当作从上而下整个进化枝的颜色，所以我们这里设置根的 的颜色会将整个树的颜色变为灰色。我们能通过在树中下面分支赋值不同的颜色 来重新定义某个分支的颜色。
让我们先定位“E”和“F”最近祖先（MRCA）节点。方法 common_ancestor 返回 原始树中这个进化枝的引用，所以当我们设置该进化枝为“salmon”颜色时，这个颜 色则会在原始的树中显示出来。
mrca = tree.common_ancestor({"name": "E"}, {"name": "F"})
mrca.color = "salmon"
当我们碰巧明确地知道某个进化枝在树中的位置，以嵌套列表的形式，我们就能 通过索引的方式直接跳到那个位置。这里，索引 [0,1] 表示根节点的第一个 子代节点的第二个子代。
tree.clade[0,1].color = "blue"
tree.clade[0,1,1].color = "yellow"
最后，展示一下我们的工作结果
Phylo.draw(tree)
注意进化枝的颜色包括导向它的分支和它的子代的分支。E和F的共同祖先结果刚好 在根分支下面，而通过这样上色，我们能清楚的看出这个树的根在哪里。
我们已经完成了很多！现在让我们休息一下，保存一下我们的工作。使用一个文件 名或句柄（这里我们使用标准输出来查看将会输出什么）和 phyloxml 格式来 调用 write 函数。PhyloXML格式保存了我们设置的颜色，所以你能通过其他树 查看工具，如Archaeopteryx，打开这个phyloXML文件，这些颜色也会显示出来。
result_phyloxml = 'result_phyloxml'
Phylo.write(tree, result_phyloxml, "phyloxml")
result_phyloxml结果如下：
<phyloxml xmlns="http://www.phyloxml.org" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.phyloxml.org http://www.phyloxml.org/1.10/phyloxml.xsd">
<phylogeny rooted="true">
<clade>
<color>
<red>128</red>
<green>128</green>
<blue>128</blue>
</color>
<clade>
<clade>
<clade>
<name>A</name>
</clade>
<clade>
<name>B</name>
</clade>
</clade>
三、I/O 函数
和SeqIO、AlignIO类似, Phylo使用四个函数处理文件的输入输出： parse 、 read 、 write 和 convert ，所有的函数都支持Newick、NEXUS、 phyloXML和NeXML等树文件格式。
read 函数解析并返回给定文件中的单个树。注意，如果文件中包含多个或不包含任何树，它将抛出一个错误。
from Bio import Phylo
tree = Phylo.read("Tests/Nexus/int_node_labels.nwk", "newick")
print tree
处理多个（或者未知个数）的树文件，需要使用 parse 函数迭代给定文件中的每一个树。
trees = Phylo.parse("Tests/PhyloXML/phyloxml_examples.xml", "phyloxml")
for tree in trees:
print tree
使用 write 函数输出一个或多个可迭代的树。
>>> trees = list(Phylo.parse("phyloxml_examples.xml", "phyloxml"))
>>> tree1 = trees[0]
>>> others = trees[1:]
>>> Phylo.write(tree1, "tree1.xml", "phyloxml")
1
>>> Phylo.write(others, "other_trees.xml", "phyloxml")
12
使用 convert 函数转换任何支持的树格式。
>>> Phylo.convert("tree1.dnd", "newick", "tree1.xml", "nexml")
1
>>> Phylo.convert("other_trees.xml", "phyloxml", "other_trees.nex", 'nexus")
12
和SeqIO和AlignIO类似，当使用字符串而不是文件作为输入输出时，需要使用 ‵‵StringIO`` 函数。
>>> from Bio import Phylo
>>> from StringIO import StringIO
>>> handle = StringIO("(((A,B),(C,D)),(E,F,G));")
>>> tree = Phylo.read(handle, "newick")
四、查看和导出树
了解一个 Tree 对象概况的最简单的方法是用 print 函数将它打印出来：
>>> tree = Phylo.read("Tests/PhyloXML/example.xml", "phyloxml")
>>> print tree
Phylogeny(rooted='True', description='phyloXML allows to use either a "branch_length"
attribute...', name='example from Prof. Joe Felsenstein's book "Inferring Phyl...')
Clade()
Clade(branch_length='0.06')
Clade(branch_length='0.102', name='A')
Clade(branch_length='0.23', name='B')
Clade(branch_length='0.4', name='C')
上面实际上是Biopython的树对象层次结构的一个概况。然而更可能的情况是，你希望见到 画出树的形状，这里有三个函数来做这件事情。
如我们在demo中看到的一样， draw_ascii 打印一个树的ascii-art图像（有根进化树） 到标准输出，或者一个打开的文件句柄，若有提供。不是所有关于树的信息被显示出来，但是它提供了一个 不依靠于任何外部依赖的快速查看树的方法。
>>> tree = Phylo.read("example.xml", "phyloxml")
>>> Phylo.draw_ascii(tree)
__________________ A
__________|
_| |___________________________________________ B
|
|___________________________________________________________________________ C
draw 函数则使用matplotlib类库画出一个更加好看的图像。查看API文档以获得关于它所接受的 用来定制输出的参数。
>>> tree = Phylo.read("example.xml", "phyloxml")
>>> Phylo.draw(tree, branch_labels=lambda c: c.branch_length)
draw_graphviz 则画出一个无根的进化分枝图（cladogram），但是它要求你安装有Graphviz、 PyDot或PyGraphviz、Network和matplotlib（或pylab）。使用上面相同的例子，和Graphviz中的 dot 程序，让我们来画一个有根树（见图. 13.3 ）：
>>> tree = Phylo.read("example.xml", "phyloxml")
>>> Phylo.draw_graphviz(tree, prog='dot')
>>> import pylab
>>> pylab.show() # Displays the tree in an interactive viewer
>>> pylab.savefig('phylo-dot.png') # Creates a PNG file of the same graphic
（提示：如果你使用 -pylab 选项执行IPython，调用 draw_graphviz 将导致matplotlib 查看器自动运行，而不需要手动的调用 show() 方法。）
这将输出树对象到一个NetworkX图中，使用Graphviz来布局节点的位置，并使用matplotlib来显示 它。这里有几个关键词参数来修改结果图像，包括大多数被NetworkX函数 networkx.draw 和 networkx.draw_graphviz 所接受的参数。
最终的显示也受所提供的树对象的 rooted 属性的影响。有根树在每个分支（branch）上显示 一个“head”来表明它的方向（见图. 13.3 ）：
>>> tree = Phylo.read("simple.dnd", "newick")
>>> tree.rooted = True
>>> Phylo.draw_graphiz(tree)
五、我的案例
5.1 继续序列，生成进化树，并对某些选中的序列名进行颜色标记
def build_pho(fasta_file,mark_color_fp,phy_png='tree.png',phy_html='tree.html'):
max_seqname_len = 30
mark_colors = {}
if os.path.exists(mark_color_fp):
with open(mark_color_fp) as data1:
for each_line in data1:
if each_line.strip() == '':
continue
mark_colors[each_line.strip()[:max_seqname_len]] = 'red'
else:
mark_colors = {}
aln_file = '.'.join(fasta_file.split('.')[:-1]) + '.aln'
# mutiple sequence alignment using clustalw2
cline = ClustalwCommandline("clustalw2", infile=fasta_file) # you need download clustalw2
stdout, stderr = cline() # run clustalw2, will generate aln_file. clustalw2 max_seqname_len = 30
# read in the alignment
aln = AlignIO.read(aln_file, "clustal")
# generate tree from identity of alignment
calculator = DistanceCalculator('identity') # identity
# calculator = DistanceCalculator('blosum62') # similarity
constructor = DistanceTreeConstructor(calculator, 'upgma')
tree = constructor.build_tree(aln)
for clade in tree.get_nonterminals():
clade.name = '' # get rid of the "inner" names added automatically
print(tree)
# draw phylo plot
n_terms = tree.count_terminals() # number of terminals to show
height = int(n_terms / 3) + 5 # height grows with number of seqs
# pylab.rcParams['figure.figsize'] = (15.0, height) # scalable figure size
plt.rcParams['figure.figsize'] = (20.0, height) # scalable figure size
plt.rcParams["lines.linewidth"] = 1
plt.rcParams["font.size"] = 12 # scalable figure size
Phylo.draw(tree, do_show=False, branch_labels=lambda c: "%.2f%%" % (float(c.branch_length) * 100),
label_colors=mark_colors) # label format
plt.savefig(phy_png)
参考资料
