用于计算分子生物学的 Python 工具

Biopython 版本 1.84
© 2024. 保留所有权利。
在 GitHub 上编辑此页面

Bio.Phylo 食谱。

以下是一些使用 Bio.Phylo 执行一些常见任务的示例。其中一些功能可能会在 Biopython 的后续版本中添加，但您可以在 Biopython 1.54 中使用它们。

便捷功能

获取分支的父节点

Bio.Phylo 中的树数据结构不存储每个分支的父节点引用。相反，可以使用 get_path 方法追踪从树根到目标分支的父节点-子节点链接路径。

def get_parent(tree, child_clade):
    node_path = tree.get_path(child_clade)
    return node_path[-2]


# Select a clade
myclade = tree.find_clades("foo").next()
# Test the function
parent = get_parent(tree, myclade)
assert myclade in parent

请注意，get_path 的运行时间与树的大小成线性关系，即为了获得最佳性能，不要在时间关键循环中调用 get_parent 或 get_path。如果可能，在循环外调用 get_path，并在该函数返回的列表中查找父节点。

或者，如果您需要重复查找任意树元素的父节点，请创建一个字典，将所有节点映射到其父节点。

def all_parents(tree):
    parents = {}
    for clade in tree.find_clades(order="level"):
        for child in clade:
            parents[child] = clade
    return parents


# Example
parents = all_parents(tree)
myclade = tree.find_clades("foo").next()
parent_of_myclade = parents[myclade]
assert myclade in parent_of_myclade

按名称索引分支

对于大型树，能够通过名称或其他唯一标识符选择分支，而不是在每次操作期间搜索整个树，这非常有用。

def lookup_by_names(tree):
    names = {}
    for clade in tree.find_clades():
        if clade.name:
            if clade.name in names:
                raise ValueError("Duplicate key: %s" % clade.name)
            names[clade.name] = clade
    return names

现在您可以以常数时间检索一个分支。

tree = Phylo.read("ncbi_taxonomy.xml", "phyloxml")
names = lookup_by_names(tree)
for phylum in ("Apicomplexa", "Euglenozoa", "Fungi"):
    print("Phylum size: %d" % len(names[phylum].get_terminals()))

潜在问题：上面实现的 lookup_by_names 不包含未命名的分支，通常是内部节点。我们可以通过为每个分支添加一个唯一标识符来解决这个问题。在这里，所有分支名称都以一个唯一编号作为前缀（这也对搜索有用）。

def tabulate_names(tree):
    names = {}
    for idx, clade in enumerate(tree.find_clades()):
        if clade.name:
            clade.name = "%d_%s" % (idx, clade.name)
        else:
            clade.name = str(idx)
        names[clade.name] = clade
    return names

计算相邻末端之间的距离

由 Joel Berendzen 提出

import itertools


def terminal_neighbor_dists(self):
    """Return a list of distances between adjacent terminals."""

    def generate_pairs(self):
        pairs = itertools.tee(self)
        next(pairs[1])
        return zip(pairs[0], pairs[1])

    return [self.distance(*i) for i in generate_pairs(self.find_clades(terminal=True))]

测试“半末端”分支

由 Joel Berendzen 提出

现有的树方法 is_preterminal 如果所有直接子节点都是末端节点，则返回 True。这段代码片段将返回 True，如果任何直接子节点是末端节点，但如果给定分支本身是末端节点，则仍然返回 False。

def is_semipreterminal(clade):
    """True if any direct descendent is terminal."""
    for child in clade:
        if child.is_terminal():
            return True
    return False

在 Python 2.5 及更高版本中，这可以通过内置的 any 函数简化。

def is_semipreterminal(clade):
    return any(child.is_terminal() for child in clade)

比较树

待办事项

共识方法

待办事项

生根方法

Tree 类（而不是 TreeMixin）上的基本方法是 root_with_outgroup

tree = Phylo.read("example.nwk", "newick")
print(tree)
# ...
tree.root_with_outgroup({"name": "A"})  # Operates in-place
print(tree)

通常，您希望外群是一个单系群，而不是单个分类单元。这不会自动检查，但您可以使用 is_monophyletic 方法自行检查。

为了节省一些输入，请尝试将查询保存在一个单独的列表中并重复使用它。

outgroup = [{"name": taxon_name} for taxon_name in ("E", "F", "G")]
if tree.is_monophyletic(outgroup):
    tree.root_with_outgroup(*outgroup)
else:
    raise ValueError("outgroup is paraphyletic")

待办事项

图形

待办事项

导出到其他类型

通过 Rpy2 转换为“ape”树

R 统计编程环境通过 ape 包和其他一些在 ape 之上构建的包提供对系统发育学的支持。Python 包 rpy2 提供了 R 和 Python 之间的接口，因此可以将 Bio.Phylo 树转换为 ape 树对象。

import tempfile
from rpy2.robjects import r


def to_ape(tree):
    """Convert a tree to the type used by the R package `ape`, via rpy2.

    Requirements:
        - Python package `rpy2`
        - R package `ape`
    """
    with tempfile.NamedTemporaryFile() as tmpf:
        Phylo.write(tree, tmpf, "newick")
        tmpf.flush()
        rtree = r(
            """
            library('ape')
            read.tree('%s')
            """
            % tmpf.name
        )
    return rtree

看看它是否有效。

>>> from StringIO import StringIO
>>> from Bio import Phylo
>>> tree = Phylo.read(StringIO("(A,(B,C),(D,E));"), "newick")
>>> rtree = to_ape(tree)
>>> len(rtree)
3
>>> print(r.summary(rtree))
Phylogenetic tree: structure(list(edge = structure(c(6, 6, 7, 7, 6, 8, 8, 1, 7,  2, 3, 8, 4, 5),
 .Dim = c(7L, 2L)), tip.label = c("A", "B", "C",  "D", "E"), Nnode = 3L),
 .Names = c("edge", "tip.label", "Nnode" ), class="phylo")

  Number of tips: 5
  Number of nodes: 3
  No branch lengths.
  No root edge.
  Tip labels: A
              B
              C
              D
              E
  No node labels.
NULL
>>> r.plot(rtree)

有关更多指南，请参阅 rpy2 文档。

转换为 DendroPy 或 PyCogent 树

Biopython、DendroPy 和 PyCogent 使用的树对象是不同的。尽管如此，所有三个工具包都支持 Newick 文件格式，因此可以通过使用一个库写入临时文件或 StringIO 对象，然后使用另一个库读取相同的字符串，在该级别实现互操作性。

from Bio import Phylo
import cogent

Phylo.write(bptree, "mytree.nwk", "newick")  # Biopython tree
ctree = cogent.LoadTree("mytree.nwk")  # PyCogent tree

import dendropy

# Create or load a tree in DendroPy
dtree = dendropy.Tree.get_from_string("(A, (B, C), (D, E))", "newick")
dtree.write_to_path("tmp.nwk", "newick", suppress_rooting=True)
# Load the same tree in Biopython
bptree = Phylo.read("tmp.nwk", "newick")

转换为 NumPy 数组或矩阵

邻接矩阵：如果存在父节点-子节点关系，则单元格为 1（真），否则为 0（假）。

import numpy as np


def to_adjacency_matrix(tree):
    """Create an adjacency matrix (NumPy array) from clades/branches in tree.

    Also returns a list of all clades in tree ("allclades"), where the position
    of each clade in the list corresponds to a row and column of the numpy
    array: a cell (i,j) in the array is 1 if there is a branch from allclades[i]
    to allclades[j], otherwise 0.

    Returns a tuple of (allclades, adjacency_matrix) where allclades is a list
    of clades and adjacency_matrix is a NumPy 2D array.
    """
    allclades = list(tree.find_clades(order="level"))
    lookup = {}
    for i, elem in enumerate(allclades):
        lookup[elem] = i
    adjmat = np.zeros((len(allclades), len(allclades)))
    for parent in tree.find_clades(terminal=False, order="level"):
        for child in parent.clades:
            adjmat[lookup[parent], lookup[child]] = 1
    if not tree.rooted:
        # Branches can go from "child" to "parent" in unrooted trees
        adjmat = adjmat + adjmat.transpose()
    return (allclades, np.matrix(adjmat))

距离矩阵：如果存在分支，则单元格值为分支长度，否则为无穷大（这与图算法配合良好）。

import numpy as np


def to_distance_matrix(tree):
    """Create a distance matrix (NumPy array) from clades/branches in tree.

    A cell (i,j) in the array is the length of the branch between allclades[i]
    and allclades[j], if a branch exists, otherwise infinity.

    Returns a tuple of (allclades, distance_matrix) where allclades is a list of
    clades and distance_matrix is a NumPy 2D array.
    """
    allclades = list(tree.find_clades(order="level"))
    lookup = {}
    for i, elem in enumerate(allclades):
        lookup[elem] = i
    distmat = np.repeat(np.inf, len(allclades) ** 2)
    distmat.shape = (len(allclades), len(allclades))
    for parent in tree.find_clades(terminal=False, order="level"):
        for child in parent.clades:
            if child.branch_length:
                distmat[lookup[parent], lookup[child]] = child.branch_length
    if not tree.rooted:
        distmat += distmat.transpose
    return (allclades, np.matrix(distmat))

增强

待办事项