计算分子生物学的 Python 工具
计算分子生物学的 Python 工具
此模块提供用于处理系统发育树的类、函数和 I/O 支持。
有关更完整的文档，请参阅 Biopython 教程 的系统发育学章节和从源代码生成的 Bio.Phylo API 页面。 Phylo 食谱 页面提供了更多有关如何使用此模块的示例，而 PhyloXML 页面描述了如何将图形提示和附加信息附加到树上。
此模块包含在 Biopython 1.54 及更高版本中。如果您有兴趣在下一个正式版本发布之前测试此代码的更新内容，请参阅 源代码 以获取有关获取开发分支副本的说明。
要绘制树（可选），您还需要这些软件包
to_networkx（以及已弃用的函数 draw_graphviz）to_networkx（以及已弃用的函数 draw_graphviz）I/O 和树操作功能无需它们；它们在调用函数 draw()、draw_graphviz() 和 to_networkx() 时按需导入。
Phylo 模块还在 Jython 2.5.1 上成功测试过，减去基于 Graphviz 和 NetworkX 的函数。但是，解析 phyloXML 文件明显更慢，因为 Jython 使用了不同版本的底层 XML 解析库。
支持的文件格式的包装器可从模块的顶层获得
from Bio import Phylo
与 SeqIO 和 AlignIO 一样，此模块提供四个 I/O 函数：parse()、read()、write() 和 convert()。每个函数都接受文件名或打开的文件句柄，因此数据也可以从压缩文件、StringIO 对象等加载。如果文件名作为字符串传递，则在函数完成时会自动关闭文件；否则，您有责任自己关闭句柄。
每个函数的第二个参数是目标格式。当前，支持以下格式
有关使用树对象的更多示例，请参阅 PhyloXML 页面。
增量解析给定文件或句柄中的每棵树，返回一个 Tree 对象迭代器（即 Bio.Phylo.BaseTree Tree 类的某个子类，具体取决于文件格式）。
>>> trees = Phylo.parse("phyloxml_examples.xml", "phyloxml")
>>> for tree in trees:
... print(tree.name)
...
example from Prof. Joe Felsenstein's book "Inferring Phylogenies"
example from Prof. Joe Felsenstein's book "Inferring Phylogenies"
same example, with support of type "bootstrap"
same example, with species and sequence
same example, with gene duplication information and sequence relationships
similar example, with more detailed sequence data
network, node B is connected to TWO nodes: AB and C
...
如果只有一棵树，则对所得生成器的 next() 方法将返回它。
>>> tree = Phylo.parse('phyloxml_examples.xml', 'phyloxml').next()
>>> tree.name
'example from Prof. Joe Felsenstein\'s book "Inferring Phylogenies"'
请注意，这不会立即显示是否存在任何剩余的树——如果您想验证这一点，请使用 read() 而不是。
从给定文件或句柄中解析并返回恰好一棵树。如果文件包含零棵或多棵树，则会引发 ValueError。如果您知道一个文件只包含一棵树，这很有用，可以用来直接加载该树对象，而不是通过 parse() 和 next() 加载，以及作为安全检查以确保输入文件实际上在顶层恰好包含一棵系统发育树。
tree = Phylo.read("example.dnd", "newick")
print(tree)
如果您已经将树数据作为 Python 字符串加载，则可以使用 StringIO（在 Python 的标准库中）解析它
from cStringIO import StringIO
treedata = "(A, (B, C), (D, E))"
handle = StringIO(treedata)
tree = Phylo.read(handle, "newick")
在一行中
tree = Phylo.read(StringIO("(A, (B, C), (D, E))"), "newick")
其他 I/O 函数也可以与 StringIO 一起使用。
（一般提示：如果您写入 StringIO 对象并想要重新读取内容，则需要调用 seek(0) 方法将句柄移回 StringIO 数据的开头——与打开的文件句柄相同。在 Biopython 源代码中 针对 Phylo 的单元测试 中查看此示例）。
将一系列 Tree 对象写入给定文件或句柄。传递单个 Tree 对象而不是列表或可迭代对象也将起作用（请参阅，Phylo 很友好）。
tree1 = Phylo.read("example1.xml", "phyloxml")
tree2 = Phylo.read("example2.xml", "phyloxml")
Phylo.write([tree1, tree2], "example-both.xml", "phyloxml")
给定两个文件（或句柄）和两种格式，两者都受 Bio.Phylo 支持，将第一个文件从第一种格式转换为第二种格式，将输出写入第二个文件。
Phylo.convert("example.nhx", "newick", "example2.nex", "nexus")
在 Phylo 模块中是特定文件格式的解析器和写入器，符合基本的顶层 API，有时还会添加其他功能。
PhyloXMLIO： 支持 phyloXML 格式。有关详细信息，请参阅 PhyloXML 页面。
NeXMLIO： 支持 NeXML 格式。
NewickIO： Bio.Nexus.Trees 中解析器的端口，通过 Phylo API 支持 Newick（又名新罕布什尔）格式。
NexusIO： 围绕 Bio.Nexus 的包装器，以支持 Nexus 树格式。
CDAOIO： 支持比较数据分析本体 (CDAO)。需要 RDFlib。
Nexus 格式实际上包含几种用于不同类型数据的子格式；为了表示树，Nexus 提供了一个包含一些元数据和一个或多个 Newick 树的块（另一种 Nexus 块可以表示比对；这在 AlignIO 中处理。因此，要解析包含所有块类型的完整 Nexus 文件，请直接使用 Bio.Nexus，要仅提取树，请使用 Bio.Phylo。
基本对象在 Bio.Phylo.BaseTree 中定义。
为了支持特定文件格式中存储的附加信息，Tree 中的子模块提供了从 BaseTree 类继承的附加类。
每个 BaseTree.Tree 或 Node 的子类都有一个类方法来将对象从基本类型提升到特定于格式的类型。这些子类对象通常可以被视为基本类型的实例，无需任何显式转换。
PhyloXML：支持 phyloXML 格式。有关详细信息，请参阅 PhyloXML 页面。
Newick：Newick 模块对 BaseTree 类提供了一些小幅改进，以及一些与现有 Bio.Nexus 模块兼容的垫片。此模块的 API 正在开发中，不应依赖于它，除了 BaseTree 已提供的功能以外。
一些额外的工具位于 Bio.Phylo 下的 Utils 模块中。这些函数也会在导入时加载到 Phylo 模块的顶层，以便于访问。
如果需要第三方软件包，则在调用函数本身时导入该软件包，因此这些依赖项对于安装或使用 Tree 模块的其余部分而言并非必需。
str(tree) 生成整个树的纯文本表示形式。字符串会自动截断，以确保合理的显示。
将此与 print 函数一起使用，以快速了解您的树
>>> tree = Phylo.parse('phyloxml_examples.xml', 'phyloxml').next()
>>> print(tree)
Phylogeny(description='phyloXML allows to use either a "branch_length"
attribute or element to indicate branch lengths.', name='example from
Prof. Joe Felsenstein s book "Inferring Phylogenies"')
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')
...
draw() 使用 Matplotlib 或 Pylab 显示一个有根的系统发育图。Biopython 1.58 中的新增功能。
尝试一下
tree = Phylo.read("apaf.xml", "phyloxml")
tree.ladderize() # Flip branches so deeper clades are displayed at top
Phylo.draw(tree)
draw_graphviz 模仿同名 networkx 函数，并进行了一些调整以改善图形的显示。如果指定了文件名，则图形将直接绘制到该文件，并且可以使用诸如图像格式（默认值为 PDF）之类的选项。
遗憾的是，draw_graphviz 绘制的图具有误导性，因此我们已弃用此方法。分支长度被忽略，并且图中节点之间的距离是任意的，这取决于 graphviz 布局引擎。但它乍一看就像一个合适的径向系统发育图，这可能会导致对数据的错误解释。用户最好使用其他库（如 ETE 或 DendroPy）创建径向图，或者只使用 Phylo.draw 生成的简单矩形图。
先决条件：除了 NetworkX 之外，您还需要在本地安装 Graphviz、Matplotlib 以及 PyGraphviz 或 pydot。
绘制基本树状图很简单
import pylab
tree = Phylo.read("apaf.xml", "phyloxml")
Phylo.draw_graphviz(tree)
pylab.show()
这是同一棵树，没有每个标记节点上的圆圈
Phylo.draw_graphviz(tree, node_size=0)
有关更多绘图功能和选项，请参阅 Phylo 食谱 页面。
draw_ascii 将 ASCII 艺术有根系统发育图打印到标准输出，或在指定的情况下打印到另一个文件句柄。仅显示终端节点标签；这些是 str(clade) 的结果（通常是进化枝名称）。用于绘制的文本字段的宽度默认情况下为 80 个字符，可以使用 column_width 关键字参数进行调整，字符行的高度是树中终端数量的两倍。
具有定义的分支长度的简单树如下所示
>>> tree = Phylo.parse('phyloxml_examples.xml', 'phyloxml').next()
>>> Phylo.draw_ascii(tree)
_____________ A
_______|
_| |_______________________________ B
|
|_______________________________________________________ C
没有分支长度的相同拓扑结构以等长分支绘制
___________________________ A
___________________________|
_| |___________________________ B
|
|___________________________ C
使用默认列宽绘制的较大树（apaf.xml，31 个叶节点）展示了如何处理相对较短的分支
>>> apaf = Phylo.read('apaf.xml', 'phyloxml')
>>> Phylo.draw_ascii(apaf)
_ 22_MOUSE
|
_| Apaf-1_HUMAN
| |
,| | 12_CANFA
||
_||___ 11_CHICK
| |
| |___________ 16_XENLA
_______|
| | , 14_FUGRU
| | __|
| |____| |__ 15_TETNG
_____| |
| | |____ 17_BRARE
| |
| | ______ 1_BRAFL
| | __|
______| || |_________ 18_NEMVE
| | |
| | |____________ 23_STRPU
| |
_| | _________ 26_STRPU
| | |_________|
| | |________ 25_STRPU
| |
| | ___ CED4_CAEEL
| |___________________________________|
____| |_ 31_CAEBR
| |
| | ___ 28_DROPS
| | _____________________|
| | ______| |____ Dark_DROME
| | | |
| |__| |_______________________ 29_AEDAE
| |
| |__________________________ 30_TRICA
|
| _ 34_BRAFL
| _________________________|
_| _____| |_ 35_BRAFL
| | |
| __| |_______________ 8_BRAFL
| | |
| | | ___________________ 20_NEMVE
| ______________| |_______|
| | | |__________________________ 21_NEMVE
| | |
| ___| |______________________________ 9_BRAFL
| | |
| | | _____________ 3_BRAFL
| | | _____|
| | |_________| |_________________ 2_BRAFL
|____| |
| |_______________ 19_NEMVE
|
| _____ 37_BRAFL
| ________________________|
|___________| |____ 36_BRAFL
|
|______________________ 33_BRAFL
虽然任何系统发育树都可以合理地用有向无环图表示，但 Phylo 模块没有尝试提供一个通用的图形库——只有表示系统发育树的最小功能。相反，它提供了一些函数，可以使用第三方库将 Tree 对象导出到标准图形表示形式、邻接列表（字典）和邻接矩阵。
to_networkx 将给定的树作为 NetworkX LabeledDiGraph 或 LabeledGraph 对象返回（取决于树是否已根植）。您可能需要直接导入 NetworkX 以便对图对象进行后续操作。从这一点开始，您也可以尝试使用 NetworkX 的绘图函数之一来显示树，并且对于简单、完全标记的树，它甚至可能有效 - 但您将通过 Phylo 自身的 draw_graphviz 函数获得更好的结果，如上所述。
import networkx, pylab
tree = Phylo.read("example.xml", "phyloxml")
net = Phylo.to_networkx(tree)
networkx.draw(net)
pylab.show()
导出到其他库（包括 **ape**（通过 Rpy2）和 **PyCogent**）的食谱可在 Phylo 食谱 页面上找到。
Yanbo Ye 为 Google Summer of Code 2013 开发了许多构建和处理系统发育树的新功能。它们在开发分支上可用（参见 SourceCode），但尚未包含在正式的 Biopython 版本中。请注意，这些功能的行为和 API 可能在即将发布的正式版本之前发生变化。
除了树构建程序的包装器（通过 Bio.Emboss.Applications 中的 EMBOSS 包装器提供的 PHYLIP 程序）外，Biopython 现在还在纯 Python 中提供了几个树构建算法实现，位于 Bio.Phylo.TreeConstruction 模块中。
所有算法都被设计为基类 TreeConstructor 的工作子类。所有构造函数都具有相同的 build_tree 方法，该方法接受一个 MultipleSeqAlignment 对象来构建树。目前有两种类型的树构造函数：DistanceTreeConstructor 和 ParsimonyTreeConstructor。
DistanceTreeConstructor 有两种算法：UPGMA（带算术平均值的未加权配对组方法）和 NJ（邻接连接）。
这两种算法都是基于距离矩阵构建树的。因此，在使用这些算法之前，让我介绍一下 DistanceCalculator，它可以从 MultipleSeqAlignment 对象生成距离矩阵。以下代码展示了一种常见的做法
>>> from Bio.Phylo.TreeConstruction import DistanceCalculator
>>> from Bio import AlignIO
>>> aln = AlignIO.read('Tests/TreeConstruction/msa.phy', 'phylip')
>>> print(aln)
Alignment with 5 rows and 13 columns
AACGTGGCCACAT Alpha
AAGGTCGCCACAC Beta
GAGATTTCCGCCT Delta
GAGATCTCCGCCC Epsilon
CAGTTCGCCACAA Gamma
>>> calculator = DistanceCalculator('identity')
>>> dm = calculator.get_distance(aln)
>>> dm
DistanceMatrix(names=['Alpha', 'Beta', 'Gamma', 'Delta', 'Epsilon'], matrix=[[0], [0.23076923076923073, 0], [0.3846153846153846, 0.23076923076923073, 0], [0.5384615384615384, 0.5384615384615384, 0.5384615384615384, 0], [0.6153846153846154, 0.3846153846153846, 0.46153846153846156, 0.15384615384615385, 0]])
>>> print(dm)
Alpha 0
Beta 0.230769230769 0
Gamma 0.384615384615 0.230769230769 0
Delta 0.538461538462 0.538461538462 0.538461538462 0
Epsilon 0.615384615385 0.384615384615 0.461538461538 0.153846153846 0
Alpha Beta Gamma Delta Epsilon
如您所见，我们使用字符串 'identity' 创建一个 DistanceCalculator 对象，该字符串是用于计算距离的模型（评分矩阵）的名称。'identity' 模型是默认模型，可用于 DNA 和蛋白质序列。要检查 DNA、蛋白质或所有可用模型，请分别使用计算器 dna_models、protein_models、models 属性。
在使用模型创建计算器后，只需使用 get_distance() 方法即可获得给定比对对象的距离矩阵。然后，您将获得一个 DistanceMatrix 对象，它是 Matrix 的子类（稍后我们将讨论它）。
现在，让我们回到 DistanceTreeConstructor。我们可以将 DistanceCalculator 对象和一个字符串参数（'nj' 或 'upgma'）传递给它进行初始化，然后调用之前提到的 build_tree()。
>>> from TreeConstruction import DistanceTreeConstructor
>>> constructor = DistanceTreeConstructor(calculator, 'nj')
>>> tree = constructor.build_tree(aln)
>>> print(tree)
Tree(rooted=False)
Clade(branch_length=0, name='Inner3')
Clade(branch_length=0.182692307692, name='Alpha')
Clade(branch_length=0.0480769230769, name='Beta')
Clade(branch_length=0.0480769230769, name='Inner2')
Clade(branch_length=0.278846153846, name='Inner1')
Clade(branch_length=0.0512820512821, name='Epsilon')
Clade(branch_length=0.102564102564, name='Delta')
Clade(branch_length=0.144230769231, name='Gamma')
虽然有时您可能希望直接使用自己的 DistanceMatrix 而不是原始比对，但我们提供另一种直接使用这两种算法的方法。
>>> from TreeConstruction import DistanceTreeConstructor
>>> constructor = DistanceTreeConstructor()
>>> tree = constructor.nj(dm)
>>> print(tree)
Tree(rooted=False)
Clade(branch_length=0, name='Inner3')
Clade(branch_length=0.182692307692, name='Alpha')
Clade(branch_length=0.0480769230769, name='Beta')
Clade(branch_length=0.0480769230769, name='Inner2')
Clade(branch_length=0.278846153846, name='Inner1')
Clade(branch_length=0.0512820512821, name='Epsilon')
Clade(branch_length=0.102564102564, name='Delta')
Clade(branch_length=0.144230769231, name='Gamma')
>>> tree = constructor.upgma(dm)
>>> print(tree)
Tree(rooted=True)
Clade(branch_length=0, name='Inner4')
Clade(branch_length=0.1875, name='Inner1')
Clade(branch_length=0.0769230769231, name='Epsilon')
Clade(branch_length=0.0769230769231, name='Delta')
Clade(branch_length=0.110576923077, name='Inner3')
Clade(branch_length=0.0384615384615, name='Inner2')
Clade(branch_length=0.115384615385, name='Gamma')
Clade(branch_length=0.115384615385, name='Beta')
Clade(branch_length=0.153846153846, name='Alpha')
与 DistanceTreeConstructor 不同，ParsimonyTreeConstructor 的具体算法委托给两个不同的工作类：ParsimonyScorer 用于根据给定的比对计算目标树的简约得分，以及 TreeSearcher 用于搜索最优树，从而最小化简约得分。一个典型的使用示例如下
>>> from Bio import AlignIO
>>> from TreeConstruction import *
>>> aln = AlignIO.read(open('Tests/TreeConstruction/msa.phy'), 'phylip')
>>> starting_tree = Phylo.read('Tests/TreeConstruction/nj.tre', 'newick')
>>> scorer = ParsimonyScorer()
>>> searcher = NNITreeSearcher(scorer)
>>> constructor = ParsimonyTreeConstructor(searcher, starting_tree)
>>> pars_tree = constructor.build_tree(aln)
>>> print(pars_tree)
Tree(weight=1.0, rooted=True)
Clade(branch_length=0.0)
Clade(branch_length=0.197335, name='Inner1')
Clade(branch_length=0.13691, name='Delta')
Clade(branch_length=0.08531, name='Epsilon')
Clade(branch_length=0.041935, name='Inner2')
Clade(branch_length=0.01421, name='Inner3')
Clade(branch_length=0.17523, name='Gamma')
Clade(branch_length=0.07477, name='Beta')
Clade(branch_length=0.29231, name='Alpha')
ParsimonyScorer 是 Fitch 算法和 Sankoff 算法的组合。如果没有提供参数，它将默认作为 Fitch 算法运行，如果提供了简约评分矩阵（一个 Matrix 对象），它将作为 Sankoff 算法运行。
然后，将评分器传递给 TreeSearcher 以告知它在搜索期间如何评估不同的树。目前，仅实现了一个搜索器 NNITreeSearcher，即最近邻交换 (NNI) 算法。
通过将搜索器和起始树传递给 ParsimonyTreeConstructor，我们最终获得了它的实例。如果没有提供起始树，将改为创建一个简单的 upgma 树，使用 'identity' 模型。要使用此简约构造函数，只需简单地使用比对调用 build_tree 方法。
与树构建算法相同，纯 Python 中的三个共识树算法（严格、多数规则和 Adam 共识）也在 Bio.Phylo.Consensus 模块中实现。您可以直接调用 strict_consensus、majority_consensus 和 adam_consensus 来使用这些算法，以树列表作为输入。
from Bio import Phylo
from Bio.Phylo.Consensus import *
trees = list(Phylo.parse("Tests/TreeConstruction/trees.tre", "newick"))
strict_tree = strict_consensus(trees)
majority_tree = majority_consensus(trees, 0.5)
adam_tree = adam_consensus(trees)
majority_consensus 方法允许您通过提供额外的 cutoff 参数（0~1，默认值为 0）来设置自己的截止值，而不是使用 50% 作为截止值。这意味着当 cutoff 等于 1 时，它也可以作为严格共识算法运行。与严格共识树和 Adam 共识树的不同之处在于，多数规则共识树的结果具有在计算过程中自动分配的分支支持值。
要获得共识树，我们必须构建引导复制树列表。因此，在 Bio.Phylo.Consensus 模块中，我们还提供了一些有用的引导方法来实现这一点。
from Bio import Phylo
from Bio.Phylo.Consensus import *
msa = AlignIO.read("Tests/TreeConstruction/msa.phy", "phylip")
msas = bootstrap(msa, 100)
如您所见，bootstrap 方法接受一个 MultipleSeqAlignment 对象，并生成其引导复制 100 次。然后，您可以使用它们来构建复制树。同时，我们还提供了一种方便的方法来做到这一点。
calculator = DistanceCalculator("blosum62")
constructor = DistanceTreeConstructor(calculator)
trees = bootstrap_trees(msa, 100, constructor)
这一次，我们将额外的 DistanceTreeConstructor 对象传递给 bootstrap_trees 方法，并最终获得了复制树。请注意，bootstrap 和 bootstrap_trees 都是生成器函数。您需要使用 list() 函数将结果转换为比对或树列表。
另一个有用的函数是 bootstrap_consensus。通过将共识方法作为另一个额外参数传递，我们可以直接获得共识树。
consensus_tree = bootstrap_consensus(msa, 100, constructor, majority_consensus)
要获得特定树的分支支持，我们可以使用 get_support 方法。
from Bio import Phylo
from Bio.Phylo.Consensus import *
trees = list(Phylo.parse("Tests/TreeConstruction/trees.tre", "newick"))
target_tree = trees[0]
support_tree = get_support(target_tree, trees)
在上面的代码中，我们使用第一棵树作为我们要计算其分支支持的目标树。get_support 方法接受目标树和树列表，并返回一棵树，所有内部分支都被分配了分支支持值。
在 TreeConstruction 和 Consensus 模块中，还有一些其他类用于这些算法中。它们可能不会经常使用，但在某些情况下可能对您有用。
_Matrix 类位于 TreeConstruction 模块中，是 _DistanceMatrix 的超类。它们实际上都是通过名称列表和嵌套数字列表（以下三角矩阵格式）构造的。它们之间的唯一区别是，_DistanceMatrix 中的对角线元素始终为 0，无论分配了哪些值。
要创建 _Matrix 对象
>>> from Bio.Phylo.TreeConstruction import _Matrix
>>> names = ['Alpha', 'Beta', 'Gamma', 'Delta']
>>> matrix = [[0], [1, 0], [2, 3, 0], [4, 5, 6, 0]]
>>> m = _Matrix(names, matrix)
>>> m
_Matrix(names=['Alpha', 'Beta', 'Gamma', 'Delta'], matrix=[[0], [1, 0], [2, 3, 0], [4, 5, 6, 0]])
>>> print(m)
Alpha 0
Beta 1 0
Gamma 2 3 0
Delta 4 5 6 0
Alpha Beta Gamma Delta
您可以使用两个索引来获取或分配矩阵中的元素，并且索引是可以交换的。
>>> m[1,2]
3
>>> m[2,1]
3
>>> m['Beta','Gamma']
3
>>> m['Beta','Gamma'] = 4
>>> m['Beta','Gamma']
4
此外，您可以使用一个索引来获取或分配与该索引相关的元素列表
>>> m[0]
[0, 1, 2, 4]
>>> m['Alpha']
[0, 1, 2, 4]
>>> m['Alpha'] = [0, 7, 8, 9]
>>> m[0]
[0, 7, 8, 9]
>>> m[0,1]
7
您还可以通过索引删除或插入元素的列和行
>>> m
_Matrix(names=['Alpha', 'Beta', 'Gamma', 'Delta'], matrix=[[0], [7, 0], [8, 4, 0], [9, 5, 6, 0]])
>>> del m['Alpha']
>>> m
_Matrix(names=['Beta', 'Gamma', 'Delta'], matrix=[[0], [4, 0], [5, 6, 0]])
>>> m.insert('Alpha', [0, 7, 8, 9] , 0)
>>> m
_Matrix(names=['Alpha', 'Beta', 'Gamma', 'Delta'], matrix=[[0], [7, 0], [8, 4, 0], [9, 5, 6, 0]])
_BitString 是 Consensus 模块中算法中经常使用的辅助类。它是 str 对象的子类，只接受两个字符（'0' 和 '1'），并具有用于二进制操作（& | ^ ~）的附加函数。常见的用法如下
>>> from Bio.Phylo.Consensus import _BitString
>>> bitstr1 = _BitString('11111')
>>> bitstr2 = _BitString('11100')
>>> bitstr3 = _BitString('01101')
>>> bitstr1
_BitString('11111')
>>> bitstr2 & bitstr3
_BitString('01100')
>>> bitstr2 | bitstr3
_BitString('11101')
>>> bitstr2 ^ bitstr3
_BitString('10001')
>>> bitstr2.index_one()
[0, 1, 2]
>>> bitstr3.index_one()
[1, 2, 4]
>>> bitstr3.index_zero()
[0, 3]
>>> bitstr1.contains(bitstr2)
True
>>> bitstr2.contains(bitstr3)
False
>>> bitstr2.independent(bitstr3)
False
>>> bitstr2.independent(bitstr4)
True
>>> bitstr1.iscompatible(bitstr2)
True
>>> bitstr2.iscompatible(bitstr3)
False
>>> bitstr2.iscompatible(bitstr4)
True
在共识和分支支持算法中，它用于对多棵树中的分支进行计数和存储。在计数过程中，如果分支的终端（就 name 属性而言）相同，则这些分支将被视为相同。
例如，假设提供了以下两棵树来搜索它们的严格共识树
tree1: (((A, B), C),(D, E))
tree2: ((A, (B, C)),(D, E))
对于这两棵树，_BitString 对象 '11111' 将表示它们的根分支。每个 '1' 代表列表 [A, B, C, D, E] 中的终端分支（顺序可能不同，由提供的第一个树的 get_terminal 方法确定）。对于树 1 中的 ((A, B), C) 分支和树 2 中的 (A, (B, C)) 分支，它们都可以表示为 '11100'。类似地，'11000' 代表树 1 中的 (A, B) 分支，'01100' 代表树 2 中的 (B, C) 分支，'00011' 代表两棵树中的 (D, E) 分支。
因此，使用此模块中的 _count_clades 函数，我们最终可以获得分支计数及其 _BitString 表示（根和终端被省略）
clade _BitString count
ABC '11100' 2
DE '00011' 2
AB '11000' 1
BC '01100' 1
要获得分支的 _BitString 表示，我们可以使用以下代码段
# suppose we are provided with a tree list, the first thing
# to do is to get all the terminal names in the first tree
term_names = [term.name for term in trees[0].get_terminals()]
# for a specific clade in any of the tree, also get its terminal names
clade_term_names = [term.name for term in clade.get_terminals()]
# then create a boolean list
boolvals = [name in clade_term_names for name in term_names]
# create the string version and pass it to _BitString
bitstr = _BitString("".join(map(str, map(int, boolvals))))
要转换回
# get all the terminal clades of the first tree
terms = [term for term in trees[0].get_terminals()]
# get the index of terminal clades in bitstr
index_list = bitstr.index_one()
# get all terminal clades by index
clade_terms = [terms[i] for i in index_list]
# create a new calde and append all the terminal clades
new_clade = BaseTree.Clade()
new_clade.clades.extend(clade_terms)
这就是 _BitString 在共识和分支支持算法中的使用方式。我认为它可以在许多其他情况下使用。
例如，我们甚至可以使用它来检查两棵树的结构是否相同。
# store and return all _BitStrings
def _bitstrs(tree):
bitstrs = set()
term_names = [term.name for term in tree.get_terminals()]
term_names.sort()
for clade in tree.get_nonterminals():
clade_term_names = [term.name for term in clade.get_terminals()]
boolvals = [name in clade_term_names for name in term_names]
bitstr = _BitString("".join(map(str, map(int, boolvals))))
bitstrs.add(bitstr)
return bitstrs
def compare(tree1, tree2):
term_names1 = [term.name for term in tree1.get_terminals()]
term_names2 = [term.name for term in tree2.get_terminals()]
# false if terminals are not the same
if set(term_names1) != set(term_names2):
return False
# true if _BitStrings are the same
if _bitstrs(tree1) == _bitstrs(tree2):
return True
else:
return False
要使用它
>>> tree1 = Phylo.read('Tests/TreeConstruction/upgma.tre', 'newick')
>>> tree2 = Phylo.read('Tests/TreeConstruction/nj.tre', 'newick')
>>> tree3 = Phylo.read('Tests/TreeConstruction/pars1.tre', 'newick')
>>> compare(tree1, tree2)
False
>>> compare(tree1, tree3)
True
>>> compare(tree2, tree3)
False
参见 Biopython 教程 中关于序列比对和 BLAST 的部分，以了解此处显示的前几个步骤的说明。
使用 BLAST 搜索蛋白质序列的同源物。
from Bio.Blast import NBCIStandalone, NCBIXML
query_fname = "AAG35789.fasta"
result_handle, error_handle = NCBIStandalone.blastall(
"/usr/bin/blastall", "blastp", "/db/fasta/swissprot", query_fname
)
blast_record = NCBIXML.read(result_handle) # This takes some time to run
从 BLAST 结果中提取最佳匹配。
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
def get_seqrecs(alignments, threshold):
for aln in alignments:
for hsp in aln.hsps:
if hsp.expect < threshold:
yield SeqRecord(Seq(hsp.sbjct), id=aln.accession)
break
best_seqs = get_seqrecs(blast_record.alignments, 1e-90)
SeqIO.write(best_seqs, "egfr-family.fasta", "fasta")
为了帮助您稍后对树进行注释，在此处选择一个查找键（例如，登录号），并构建一个字典，将该键映射到您从 BLAST 输出中获得的任何其他数据，例如分类学和 GI 号。在此示例中，我们只保留原始序列和登录号。
使用 Muscle 重新比对序列。该程序以 Clustal 格式创建比对文件，“egfr-family.aln”。
from Bio.Align.Applications import MuscleCommandline
cmdline = MuscleCommandline(input="egfr-family.fasta", out="efgr-family.aln", clw=True)
cmdline()
使用 PhyML 推断基因树。首先，将步骤 3 中的对齐结果转换为“松散的 Phylip”格式（Biopython 1.58 新增功能）。
from Bio import AlignIO
AlignIO.convert("egfr-family.aln", "clustal", "egfr-family.phy", "phylip-relaxed")
使用命令行包装器将对齐结果输入到 PhyML 中。
from Bio.Phylo.Applications import PhymlCommandline
cmdline = PhymlCommandline(
input="egfr-family.phy", datatype="aa", model="WAG", alpha="e", bootstrap=100
)
out_log, err_log = cmdline()
使用 Phylo 加载基因树，并快速查看拓扑结构。（PhyML 会将树写入以输入文件加“_phyml_tree.txt”命名的文件。）
from Bio import Phylo
egfr_tree = Phylo.read("egfr-family.phy_phyml_tree.txt", "newick")
Phylo.draw_ascii(egfr_tree)
向树中添加登录号和序列 - 现在我们正在使用 PhyloXML 的额外功能。
from Bio.Phylo import PhyloXML
# Promote the basic tree to PhyloXML
egfr_phy = egfr_tree.as_phyloxml()
# Make a lookup table for sequences
lookup = dict((rec.id, str(rec.seq)) for rec in best_seqs)
for clade in egfr_phy.get_terminals():
key = clade.name
accession = PhyloXML.Accession(key, "NCBI")
mol_seq = PhyloXML.MolSeq(lookup[key], is_aligned=True)
sequence = PhyloXML.Sequence(type="aa", accession=accession, mol_seq=mol_seq)
clade.sequences.append(sequence)
# Save the annotated phyloXML file
Phylo.write(egfr_phy, "egfr-family-annotated.xml", "phyloxml")