【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)

参考资料

药企,独角兽,苏州。团队长期招人,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn