多个序列比对对象
本章描述老年人 MultipleSeqAlignment 类和中的解析器 Bio.AlignIO 分析序列比对软件的输出，生成 MultipleSeqAlignment 对象所谓多序列对齐，我们是指已对齐在一起的多个序列的集合-通常插入空位字符，并添加引导或尾部空位-使得所有序列串的长度相同。这样的对齐可以被视为字母矩阵，其中每一行都被保存为 SeqRecord 对象内部。
我们将介绍 MultipleSeqAlignment 保存此类数据的对象，以及 Bio.AlignIO 用于将它们读写为各种文件格式的模块（遵循 Bio.SeqIO 模块（上一章）。注意这两 Bio.SeqIO 和 Bio.AlignIO 可以读写序列对齐文件。适当的选择在很大程度上取决于您想对数据做什么。
本章的最后一部分是关于使用常见的多序列比对工具，例如Python的CustalW和MUSYS，并使用Biopython解析结果。
解析或读取序列对齐
我们有两个功能用于读取序列比对， Bio.AlignIO.read() 和 Bio.AlignIO.parse() 它是在大会之后推出的 Bio.SeqIO 分别用于包含一个或多个对齐的文件。
使用 Bio.AlignIO.parse() 将返回一个 iterator 这给 MultipleSeqAlignment 对象迭代器通常用于for循环中。您将拥有多个不同比对的情况示例包括从PHYLIP工具重新采样的比对 seqboot ，或来自CLARSS工具的多个成对比对 water 或 needle ，或Bill Pearson的FASTA工具。
然而，在许多情况下，您将处理仅包含单个对齐方式的文件。在这种情况下，您应该使用 Bio.AlignIO.read() 返回单个 MultipleSeqAlignment object.
这两个函数都需要两个强制参数：
第一个参数是一个 handle 从其中读取数据，通常是打开的文件（请参阅第节 把手到底是什么？ ），或文件名。
第二个参数是指定对齐格式的大小写字符串。如在
Bio.SeqIO我们不会试图为您猜测文件格式！有关支持格式的完整列表，请参阅http://biopython.org/wiki/AlignIO。
还有一个可选的 seq_count 第一节讨论的论点 Ambient Alignments 下面用于处理可能包含多个对齐方式的模糊文件格式。
单一路线
例如，考虑以下PFAM或Stockholm文件格式中的富含注释的蛋白质排列：
# STOCKHOLM 1.0
#=GS COATB_BPIKE/30-81 AC P03620.1
#=GS COATB_BPIKE/30-81 DR PDB; 1ifl ; 1-52;
#=GS Q9T0Q8_BPIKE/1-52 AC Q9T0Q8.1
#=GS COATB_BPI22/32-83 AC P15416.1
#=GS COATB_BPM13/24-72 AC P69541.1
#=GS COATB_BPM13/24-72 DR PDB; 2cpb ; 1-49;
#=GS COATB_BPM13/24-72 DR PDB; 2cps ; 1-49;
#=GS COATB_BPZJ2/1-49 AC P03618.1
#=GS Q9T0Q9_BPFD/1-49 AC Q9T0Q9.1
#=GS Q9T0Q9_BPFD/1-49 DR PDB; 1nh4 A; 1-49;
#=GS COATB_BPIF1/22-73 AC P03619.2
#=GS COATB_BPIF1/22-73 DR PDB; 1ifk ; 1-50;
COATB_BPIKE/30-81 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA
#=GR COATB_BPIKE/30-81 SS -HHHHHHHHHHHHHH--HHHHHHHH--HHHHHHHHHHHHHHHHHHHHH----
Q9T0Q8_BPIKE/1-52 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA
COATB_BPI22/32-83 DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA
COATB_BPM13/24-72 AEGDDP...AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
#=GR COATB_BPM13/24-72 SS ---S-T...CHCHHHHCCCCTCCCTTCHHHHHHHHHHHHHHHHHHHHCTT--
COATB_BPZJ2/1-49 AEGDDP...AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA
Q9T0Q9_BPFD/1-49 AEGDDP...AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
#=GR Q9T0Q9_BPFD/1-49 SS ------...-HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH--
COATB_BPIF1/22-73 FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA
#=GR COATB_BPIF1/22-73 SS XX-HHHH--HHHHHH--HHHHHHH--HHHHHHHHHHHHHHHHHHHHHHH---
#=GC SS_cons XHHHHHHHHHHHHHHHCHHHHHHHHCHHHHHHHHHHHHHHHHHHHHHHHC--
#=GC seq_cons AEssss...AptAhDSLpspAT-hIu.sWshVsslVsAsluIKLFKKFsSKA
//
这是Phage_Coat_Gp8（PF 05371）PFAM条目的种子对齐，从https://pfam.xfam.org/的现已过时的PFAM版本下载。我们可以按以下方式加载此文件（假设它已保存到磁盘中，作为当前工作目录中的“PF05371_seed.sth”）：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
此代码将打印对齐摘要：
>>> print(alignment)
Alignment with 7 rows and 52 columns
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRL...SKA COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKL...SRA Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRL...SKA COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73
您会注意到在上面的输出中，序列已被截断。相反，我们可以编写自己的代码来通过迭代行来随意格式化它 SeqRecord 对象：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> print("Alignment length %i" % alignment.get_alignment_length())
Alignment length 52
>>> for record in alignment:
... print("%s - %s" % (record.seq, record.id))
...
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA - COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA - Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA - COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA - COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA - COATB_BPIF1/22-73
您还可以调用Python的内置 format 对齐对象上的函数以特定文件格式显示它-请参阅第节 将对齐对象作为格式化字符串 有关详细信息
您是否在上面的原始文件中注意到，其中几个序列包括对TSB和相关已知二级结构的数据库交叉引用？试试这个：
>>> for record in alignment:
... if record.dbxrefs:
... print("%s %s" % (record.id, record.dbxrefs))
...
COATB_BPIKE/30-81 ['PDB; 1ifl ; 1-52;']
COATB_BPM13/24-72 ['PDB; 2cpb ; 1-49;', 'PDB; 2cps ; 1-49;']
Q9T0Q9_BPFD/1-49 ['PDB; 1nh4 A; 1-49;']
COATB_BPIF1/22-73 ['PDB; 1ifk ; 1-50;']
要查看所有序列注释，请尝试以下操作：
>>> for record in alignment:
... print(record)
...
PFAM在http://pfam.xfam.org/family/PF05371上提供了一个很好的网络界面，实际上可以让您下载多种其他格式的对齐方式。这是FASTA文件格式中的文件外观：
>COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA
>Q9T0Q8_BPIKE/1-52
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA
>COATB_BPI22/32-83
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA
>COATB_BPM13/24-72
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
>COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA
>Q9T0Q9_BPFD/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
>COATB_BPIF1/22-73
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA
请注意，网站应该可以选择将间隙显示为句号（点）或虚线，我们在上面显示了虚线。假设您下载并保存为文件“PF05371_seed.faa”，那么您可以使用几乎完全相同的代码加载它：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.faa", "fasta")
>>> print(alignment)
此代码中唯一更改的是文件名和格式字符串。您将得到与以前相同的输出，序列和记录标识符相同。然而，正如您所期望的那样，如果您检查每个 SeqRecord 没有注释或数据库交叉引用，因为这些不包括在FASTA文件格式中。
请注意，您可以使用，而不是使用Sanger网站 Bio.AlignIO 自行将原始斯德哥尔摩格式文件转换为FASTA文件（见下文）。
对于任何受支持的文件格式，您都可以通过更改格式字符串以完全相同的方式加载路线。例如，对于PHYLIP文件使用“phylip”，对于NEXUS文件使用“nexus”，或者对于由MANUSS工具输出的比对使用“nexus”。在wiki页面（http：//biopython.org/wiki/AlignIO）和内置文档中有一个完整的列表， Bio.AlignIO :
>>> from Bio import AlignIO
>>> help(AlignIO)
多重比对
上一节重点讨论了读取包含单个对齐的文件。但一般来说，文件可以包含多个对齐方式，要读取这些文件，我们必须使用 Bio.AlignIO.parse() 功能
假设您有PHYLIP格式的小对齐：
5 6
Alpha AACAAC
Beta AACCCC
Gamma ACCAAC
Delta CCACCA
Epsilon CCAAAC
如果您想使用PHYLIP工具引导系统发生树，步骤之一是使用该工具创建一组多个重新采样的比对 bootseq .这将给出类似于这样的输出，为了简洁起见，该输出被缩写：
5 6
Alpha AAACCA
Beta AAACCC
Gamma ACCCCA
Delta CCCAAC
Epsilon CCCAAA
5 6
Alpha AAACAA
Beta AAACCC
Gamma ACCCAA
Delta CCCACC
Epsilon CCCAAA
5 6
Alpha AAAAAC
Beta AAACCC
Gamma AACAAC
Delta CCCCCA
Epsilon CCCAAC
...
5 6
Alpha AAAACC
Beta ACCCCC
Gamma AAAACC
Delta CCCCAA
Epsilon CAAACC
如果您想在使用中阅读这篇文章 Bio.AlignIO 您可以用途：
>>> from Bio import AlignIO
>>> alignments = AlignIO.parse("resampled.phy", "phylip")
>>> for alignment in alignments:
... print(alignment)
... print()
...
这将给出以下输出，再次缩写以供显示：
Alignment with 5 rows and 6 columns
AAACCA Alpha
AAACCC Beta
ACCCCA Gamma
CCCAAC Delta
CCCAAA Epsilon
Alignment with 5 rows and 6 columns
AAACAA Alpha
AAACCC Beta
ACCCAA Gamma
CCCACC Delta
CCCAAA Epsilon
Alignment with 5 rows and 6 columns
AAAAAC Alpha
AAACCC Beta
AACAAC Gamma
CCCCCA Delta
CCCAAC Epsilon
...
Alignment with 5 rows and 6 columns
AAAACC Alpha
ACCCCC Beta
AAAACC Gamma
CCCCAA Delta
CAAACC Epsilon
与功能一样 Bio.SeqIO.parse() ，使用 Bio.AlignIO.parse() 返回迭代器。如果您想将所有对齐同时保留在内存中，这将允许您以任何顺序访问它们，那么将迭代器变成列表：
>>> from Bio import AlignIO
>>> alignments = list(AlignIO.parse("resampled.phy", "phylip"))
>>> last_align = alignments[-1]
>>> first_align = alignments[0]
Ambient Alignments
许多路线文件格式可以显式存储多个路线，并且每个路线之间的划分是明确的。然而，当使用通用序列文件格式时，就不存在这样的块结构。最常见的情况是路线已以FASTA文件格式保存时。例如，考虑以下内容：
>Alpha
ACTACGACTAGCTCAG--G
>Beta
ACTACCGCTAGCTCAGAAG
>Gamma
ACTACGGCTAGCACAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Beta
ACTACCGCTAGCTCAGAAG
>Gamma
ACTACGGCTAGCACAGAAG
这可能是包含六个序列（具有重复的标识符）的单个比对。或者，从标识符来看，这可能是两个不同的比对，每个比对有三个序列，恰好都具有相同的长度。
下一个例子怎么样？
>Alpha
ACTACGACTAGCTCAG--G
>Beta
ACTACCGCTAGCTCAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Gamma
ACTACGGCTAGCACAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Delta
ACTACGGCTAGCACAGAAG
同样，这可能是具有六个序列的单一比对。然而，这一次，根据标识符，我们可能会猜测这是三个成对的比对，它们碰巧都具有相同的长度。
最后一个例子类似：
>Alpha
ACTACGACTAGCTCAG--G
>XXX
ACTACCGCTAGCTCAGAAG
>Alpha
ACTACGACTAGCTCAGG
>YYY
ACTACGGCAAGCACAGG
>Alpha
--ACTACGAC--TAGCTCAGG
>ZZZ
GGACTACGACAATAGCTCAGG
在第三个示例中，由于长度不同，因此不能将其视为包含所有六个记录的单个对齐。然而，它可能是三对对齐。
显然，尝试在FASTA文件中存储多个对齐并不理想。然而，如果您被迫将这些作为输入文件处理 Bio.AlignIO 可以应对所有路线具有相同记录数量的最常见情况。其中一个例子是成对比对的集合，可以通过CLARSS工具生成 needle 和 water - 尽管在这种情况下， Bio.AlignIO 应该能够理解使用“RST”作为格式字符串的本地输出。
为了将这些FASTA示例解释为几个单独的比对，我们可以使用 Bio.AlignIO.parse() 与可选的 seq_count 参数，指定每次比对中预期有多少个序列（在这些示例中，分别为3、2和2）。例如，使用第三个例子作为输入数据：
>>> for alignment in AlignIO.parse(handle, "fasta", seq_count=2):
... print("Alignment length %i" % alignment.get_alignment_length())
... for record in alignment:
... print("%s - %s" % (record.seq, record.id))
... print()
...
给予：
Alignment length 19
ACTACGACTAGCTCAG--G - Alpha
ACTACCGCTAGCTCAGAAG - XXX
Alignment length 17
ACTACGACTAGCTCAGG - Alpha
ACTACGGCAAGCACAGG - YYY
Alignment length 21
--ACTACGAC--TAGCTCAGG - Alpha
GGACTACGACAATAGCTCAGG - ZZZ
使用 Bio.AlignIO.read() 或 Bio.AlignIO.parse() 未经 seq_count 参数将给出包含前两个示例的所有六个记录的单个对齐。对于第三个示例，会引发异常，因为长度不同，阻止它们变成单个对齐。
如果文件格式本身具有块结构，允许 Bio.AlignIO 要直接确定每次比对中的序列数，则 seq_count 不需要争论。如果提供了它并且与文件内容不一致，则会引发错误。
请注意，这个可选的 seq_count 参数假设文件中的每个对齐都具有相同数量的序列。假设您可能会遇到陌生的情况，例如一个FASTA文件，其中包含多个比对，每个比对具有不同数量的序列-尽管我很想听到一个现实世界的例子。假设您无法以更好的文件格式获取数据，则没有直接的方法可以使用 Bio.AlignIO .在这种情况下，您可以考虑使用读取序列本身 Bio.SeqIO 并将它们叠加在一起以酌情创建路线。
书写对齐
我们已经讨论过使用 Bio.AlignIO.read() 和 Bio.AlignIO.parse() 用于对齐输入（读取文件），现在我们来看看 Bio.AlignIO.write() 用于对齐输出（写入文件）。这是一个包含三个参数的函数：一些 MultipleSeqAlignment 对象（或者为了向后兼容性，过时的 Alignment 对象）、要写入的手柄或文件名以及序列格式。
这是一个例子，我们首先创建一些 MultipleSeqAlignment 用硬方法（手工，而不是从文件加载它们）。请注意，我们创建了一些 SeqRecord 用于构建对齐的对象。
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> from Bio.Align import MultipleSeqAlignment
>>> align1 = MultipleSeqAlignment(
... [
... SeqRecord(Seq("ACTGCTAGCTAG"), id="Alpha"),
... SeqRecord(Seq("ACT-CTAGCTAG"), id="Beta"),
... SeqRecord(Seq("ACTGCTAGDTAG"), id="Gamma"),
... ]
... )
>>> align2 = MultipleSeqAlignment(
... [
... SeqRecord(Seq("GTCAGC-AG"), id="Delta"),
... SeqRecord(Seq("GACAGCTAG"), id="Epsilon"),
... SeqRecord(Seq("GTCAGCTAG"), id="Zeta"),
... ]
... )
>>> align3 = MultipleSeqAlignment(
... [
... SeqRecord(Seq("ACTAGTACAGCTG"), id="Eta"),
... SeqRecord(Seq("ACTAGTACAGCT-"), id="Theta"),
... SeqRecord(Seq("-CTACTACAGGTG"), id="Iota"),
... ]
... )
>>> my_alignments = [align1, align2, align3]
现在我们有了一个列表 Alignment 对象，我们将它们写入PHYLIP格式文件：
>>> from Bio import AlignIO
>>> AlignIO.write(my_alignments, "my_example.phy", "phylip")
如果您在您最喜欢的文本编辑器中打开此文件，它应该是这样的：
3 12
Alpha ACTGCTAGCT AG
Beta ACT-CTAGCT AG
Gamma ACTGCTAGDT AG
3 9
Delta GTCAGC-AG
Epislon GACAGCTAG
Zeta GTCAGCTAG
3 13
Eta ACTAGTACAG CTG
Theta ACTAGTACAG CT-
Iota -CTACTACAG GTG
更常见的情况是，希望加载现有的对齐并保存它，也许在进行一些简单的操作（例如删除某些行或列）之后。
假设您想知道有多少个路线 Bio.AlignIO.write() 函数被写到手柄？如果您的路线位于上面示例中，则可以使用 len(my_alignments) ，但是，当您的记录来自生成器/迭代器时，您就无法做到这一点。因此 Bio.AlignIO.write() 函数返回写入文件的对齐数。
Note - 如果你告诉 Bio.AlignIO.write() 函数写入已存在的文件，旧文件将被覆盖，而没有任何警告。
序列对齐文件格式之间转换
在序列比对文件格式之间转换 Bio.AlignIO 与在序列文件格式之间转换的方式相同 Bio.SeqIO （部分 序列文件格式之间转换 ).我们通常使用以下方式加载路线 Bio.AlignIO.parse() 然后使用 Bio.AlignIO.write() - 或只是使用 Bio.AlignIO.convert() 助手功能。
对于本示例，我们将加载之前使用的PFAM/Stockholm格式文件并将其保存为Custal W格式文件：
>>> from Bio import AlignIO
>>> count = AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.aln", "clustal")
>>> print("Converted %i alignments" % count)
Converted 1 alignments
或者使用 Bio.AlignIO.parse() 和 Bio.AlignIO.write() :
>>> from Bio import AlignIO
>>> alignments = AlignIO.parse("PF05371_seed.sth", "stockholm")
>>> count = AlignIO.write(alignments, "PF05371_seed.aln", "clustal")
>>> print("Converted %i alignments" % count)
Converted 1 alignments
的 Bio.AlignIO.write() 函数期望被赋予多个对齐对象。在上面的例子中，我们为它提供了由 Bio.AlignIO.parse() .
在这种情况下，我们知道文件中只有一个对齐方式，因此我们可以使用 Bio.AlignIO.read() 相反，但请注意，我们必须将此对齐传递给 Bio.AlignIO.write() 作为单个元素列表：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> AlignIO.write([alignment], "PF05371_seed.aln", "clustal")
无论哪种方式，您最终都应该得到相同的新Classal W格式文件“PF05371_seed.aln”，其中包含以下内容：
CLUSTAL X (1.81) multiple sequence alignment
COATB_BPIKE/30-81 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSS
Q9T0Q8_BPIKE/1-52 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVS
COATB_BPI22/32-83 DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSS
COATB_BPM13/24-72 AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTS
COATB_BPZJ2/1-49 AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFAS
Q9T0Q9_BPFD/1-49 AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTS
COATB_BPIF1/22-73 FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVS
COATB_BPIKE/30-81 KA
Q9T0Q8_BPIKE/1-52 RA
COATB_BPI22/32-83 KA
COATB_BPM13/24-72 KA
COATB_BPZJ2/1-49 KA
Q9T0Q9_BPFD/1-49 KA
COATB_BPIF1/22-73 RA
或者，您可以制作一个PHYLIP格式文件，我们将其命名为“PF05371_seed.phy”：
>>> from Bio import AlignIO
>>> AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.phy", "phylip")
这次的输出看起来像这样：
7 52
COATB_BPIK AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIRLFKKFSS
Q9T0Q8_BPI AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIKLFKKFVS
COATB_BPI2 DGTSTATSYA TEAMNSLKTQ ATDLIDQTWP VVTSVAVAGL AIRLFKKFSS
COATB_BPM1 AEGDDP---A KAAFNSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPZJ AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFAS
Q9T0Q9_BPF AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPIF FAADDATSQA KAAFDSLTAQ ATEMSGYAWA LVVLVVGATV GIKLFKKFVS
KA
RA
KA
KA
KA
KA
RA
原始PHYLIP对齐文件格式的一大障碍是序列标识符被严格截断为十个字符。在这个例子中，正如您所看到的那样，结果的名称仍然是唯一的-但它们不是很可读。因此，原始PHYLIP格式的一种更宽松的变体现在被广泛使用：
>>> from Bio import AlignIO
>>> AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.phy", "phylip-relaxed")
这次的输出看起来像这样，使用更长的凹痕来允许完整给出所有标识符：
7 52
COATB_BPIKE/30-81 AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIRLFKKFSS
Q9T0Q8_BPIKE/1-52 AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIKLFKKFVS
COATB_BPI22/32-83 DGTSTATSYA TEAMNSLKTQ ATDLIDQTWP VVTSVAVAGL AIRLFKKFSS
COATB_BPM13/24-72 AEGDDP---A KAAFNSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPZJ2/1-49 AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFAS
Q9T0Q9_BPFD/1-49 AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPIF1/22-73 FAADDATSQA KAAFDSLTAQ ATEMSGYAWA LVVLVVGATV GIKLFKKFVS
KA
RA
KA
KA
KA
KA
RA
如果您必须使用原始严格的PHYLIP格式，那么您可能需要以某种方式压缩标识符-或分配您自己的名称或编号系统。以下这段代码在保存输出之前操作记录标识符：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> name_mapping = {}
>>> for i, record in enumerate(alignment):
... name_mapping[i] = record.id
... record.id = "seq%i" % i
...
>>> print(name_mapping)
{0: 'COATB_BPIKE/30-81', 1: 'Q9T0Q8_BPIKE/1-52', 2: 'COATB_BPI22/32-83', 3: 'COATB_BPM13/24-72', 4: 'COATB_BPZJ2/1-49', 5: 'Q9T0Q9_BPFD/1-49', 6: 'COATB_BPIF1/22-73'}
>>> AlignIO.write([alignment], "PF05371_seed.phy", "phylip")
此代码使用Python字典记录从新序列系统到原始标识符的简单映射：
{
0: "COATB_BPIKE/30-81",
1: "Q9T0Q8_BPIKE/1-52",
2: "COATB_BPI22/32-83",
# ...
}
以下是新的（严格）PHYLIP格式输出：
7 52
seq0 AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIRLFKKFSS
seq1 AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIKLFKKFVS
seq2 DGTSTATSYA TEAMNSLKTQ ATDLIDQTWP VVTSVAVAGL AIRLFKKFSS
seq3 AEGDDP---A KAAFNSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
seq4 AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFAS
seq5 AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
seq6 FAADDATSQA KAAFDSLTAQ ATEMSGYAWA LVVLVVGATV GIKLFKKFVS
KA
RA
KA
KA
KA
KA
RA
一般来说，由于标识符限制，使用 strict PHYLIP文件格式不应成为您的首选。另一方面，使用PFAM/Stockholm格式也允许您记录大量额外的注释。
将对齐对象作为格式化字符串
的 Bio.AlignIO 界面基于手柄，这意味着如果您想将对齐方式转换为特定文件格式的字符串，则需要做更多的工作（见下文）。然而，您可能更喜欢调用Python的内置 format 对对齐对象的函数。这将输出格式规范作为单个参数，即由支持的大小写字符串 Bio.AlignIO 作为输出格式。例如：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> print(format(alignment, "clustal"))
CLUSTAL X (1.81) multiple sequence alignment
COATB_BPIKE/30-81 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSS
Q9T0Q8_BPIKE/1-52 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVS
COATB_BPI22/32-83 DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSS
...
如果没有输出格式规范， format 返回与 str .
节所述 format方法 ， SeqRecord 对象具有类似的方法，使用支持的输出格式 Bio.SeqIO .
境内 format 呼吁 Bio.AlignIO.write() 与 StringIO 把手进行例如，如果您使用的是较旧版本的Biopython，则可以在自己的代码中执行此操作：
>>> from io import StringIO
>>> from Bio import AlignIO
>>> alignments = AlignIO.parse("PF05371_seed.sth", "stockholm")
>>> out_handle = StringIO()
>>> AlignIO.write(alignments, out_handle, "clustal")
1
>>> clustal_data = out_handle.getvalue()
>>> print(clustal_data)
CLUSTAL X (1.81) multiple sequence alignment
COATB_BPIKE/30-81 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSS
Q9T0Q8_BPIKE/1-52 AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVS
COATB_BPI22/32-83 DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSS
COATB_BPM13/24-72 AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTS
...
操纵路线
既然我们已经介绍了加载和保存路线，我们将看看您还可以使用它们做什么。
切片对齐
首先，在某种意义上，对齐对象的行为就像Python list 的 SeqRecord 对象（行）。考虑到这个模型，希望 len() (the行数）和迭代（每行作为 SeqRecord ）有意义：
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> print("Number of rows: %i" % len(alignment))
Number of rows: 7
>>> for record in alignment:
... print("%s - %s" % (record.seq, record.id))
...
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA - COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA - Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA - COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA - COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA - COATB_BPIF1/22-73
您还可以使用类似列表的 append 和 extend 向对齐添加更多行的方法（如 SeqRecord 对象）。请记住列表比喻，对对齐方式的简单切片也应该是有意义的--它选择一些行来返回另一个对齐对象：
>>> print(alignment)
Alignment with 7 rows and 52 columns
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRL...SKA COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKL...SRA Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRL...SKA COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73
>>> print(alignment[3:7])
Alignment with 4 rows and 52 columns
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73
如果你想按列选择呢？那些使用过NumPy矩阵或数组对象的人不会对此感到惊讶-您使用的是双索引。
>>> print(alignment[2, 6])
T
使用两个整指数拉出一个字母（简写）：
>>> print(alignment[2].seq[6])
T
您可以将一列拉出为字符串，如下所示：
>>> print(alignment[:, 6])
TTT---T
您还可以选择一系列列。例如，要选出我们之前提取的相同三行，但只选取它们的前六列：
>>> print(alignment[3:6, :6])
Alignment with 3 rows and 6 columns
AEGDDP COATB_BPM13/24-72
AEGDDP COATB_BPZJ2/1-49
AEGDDP Q9T0Q9_BPFD/1-49
将第一个指数保留为 : 意味着占据所有行：
>>> print(alignment[:, :6])
Alignment with 7 rows and 6 columns
AEPNAA COATB_BPIKE/30-81
AEPNAA Q9T0Q8_BPIKE/1-52
DGTSTA COATB_BPI22/32-83
AEGDDP COATB_BPM13/24-72
AEGDDP COATB_BPZJ2/1-49
AEGDDP Q9T0Q9_BPFD/1-49
FAADDA COATB_BPIF1/22-73
这为我们带来了一种删除部分的简洁方法。请注意第7、8和9列，它们是七个序列中的三个序列中的空白：
>>> print(alignment[:, 6:9])
Alignment with 7 rows and 3 columns
TNY COATB_BPIKE/30-81
TNY Q9T0Q8_BPIKE/1-52
TSY COATB_BPI22/32-83
--- COATB_BPM13/24-72
--- COATB_BPZJ2/1-49
--- Q9T0Q9_BPFD/1-49
TSQ COATB_BPIF1/22-73
同样，您可以切片以获得第九列之后的所有内容：
>>> print(alignment[:, 9:])
Alignment with 7 rows and 43 columns
ATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA COATB_BPIKE/30-81
ATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA Q9T0Q8_BPIKE/1-52
ATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA COATB_BPI22/32-83
AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA COATB_BPM13/24-72
AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA COATB_BPZJ2/1-49
AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA Q9T0Q9_BPFD/1-49
AKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA COATB_BPIF1/22-73
现在，有趣的是，添加对齐对象是按列进行的。这可以让您这样做，作为删除列块的一种方式：
>>> edited = alignment[:, :6] + alignment[:, 9:]
>>> print(edited)
Alignment with 7 rows and 49 columns
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA COATB_BPIKE/30-81
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA Q9T0Q8_BPIKE/1-52
DGTSTAATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA COATB_BPI22/32-83
AEGDDPAKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA COATB_BPM13/24-72
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA COATB_BPZJ2/1-49
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA Q9T0Q9_BPFD/1-49
FAADDAAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA COATB_BPIF1/22-73
比对添加的另一个常见用途是将几个不同基因的比对组合成元比对。但请注意-标识符需要匹配（请参阅第节 添加SeqRecord对象 关于如何添加 SeqRecord 对象作品）。您可能会发现首先按id按字母顺序对对齐行进行排序很有帮助：
>>> edited.sort()
>>> print(edited)
Alignment with 7 rows and 49 columns
DGTSTAATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA COATB_BPI22/32-83
FAADDAAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA COATB_BPIF1/22-73
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA COATB_BPIKE/30-81
AEGDDPAKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA COATB_BPM13/24-72
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA COATB_BPZJ2/1-49
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA Q9T0Q8_BPIKE/1-52
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA Q9T0Q9_BPFD/1-49
请注意，只有当两个对齐的行数相同时，您才能将它们添加在一起。
排列为阵列
根据您正在做的事情，将对齐对象转换为字母数组可能会更有用-您可以使用NumPy做到这一点：
>>> import numpy as np
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> align_array = np.array(alignment)
>>> print("Array shape %i by %i" % align_array.shape)
Array shape 7 by 52
>>> align_array[:, :10]
array([['A', 'E', 'P', 'N', 'A', 'A', 'T', 'N', 'Y', 'A'],
['A', 'E', 'P', 'N', 'A', 'A', 'T', 'N', 'Y', 'A'],
['D', 'G', 'T', 'S', 'T', 'A', 'T', 'S', 'Y', 'A'],
['A', 'E', 'G', 'D', 'D', 'P', '-', '-', '-', 'A'],
['A', 'E', 'G', 'D', 'D', 'P', '-', '-', '-', 'A'],
['A', 'E', 'G', 'D', 'D', 'P', '-', '-', '-', 'A'],
['F', 'A', 'A', 'D', 'D', 'A', 'T', 'S', 'Q', 'A']],...
请注意，这将内存中的原始Biopython对齐对象和NumPy数组作为单独的对象-编辑其中一个不会更新另一个！
计数替换
的 substitutions 对齐的属性报告对齐中的字母相互替换的频率。计算方法是：取对齐中的所有行对，计算两个字母相互对齐的次数，然后对所有行对进行相加。例如，
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> from Bio.Align import MultipleSeqAlignment
>>> msa = MultipleSeqAlignment(
... [
... SeqRecord(Seq("ACTCCTA"), id="seq1"),
... SeqRecord(Seq("AAT-CTA"), id="seq2"),
... SeqRecord(Seq("CCTACT-"), id="seq3"),
... SeqRecord(Seq("TCTCCTC"), id="seq4"),
... ]
... )
>>> print(msa)
Alignment with 4 rows and 7 columns
ACTCCTA seq1
AAT-CTA seq2
CCTACT- seq3
TCTCCTC seq4
>>> substitutions = msa.substitutions
>>> print(substitutions)
A C T
A 2.0 4.5 1.0
C 4.5 10.0 0.5
T 1.0 0.5 12.0
由于对的顺序是任意的，因此计数在对角线上下平分。例如，9个路线 A 到 C 位置存储为4.5 ['A', 'C'] 位置为4.5 ['C', 'A'] .这种安排有助于在根据这些计数计算替换矩阵时简化数学运算，如第 取代矩阵 .
注意 msa.substitutions 仅包含对齐中出现的字母的条目。您可以使用 select 添加缺失字母条目的方法，例如
>>> m = substitutions.select("ATCG")
>>> print(m)
A T C G
A 2.0 1.0 4.5 0.0
T 1.0 12.0 0.5 0.0
C 4.5 0.5 10.0 0.0
G 0.0 0.0 0.0 0.0
这还允许您更改字母表中字母的顺序。
获取新型对齐对象
使用 alignment 房产打造新型 Alignment 对象（请参阅部分 路线对象 ）来自老式的 MultipleSeqAlignment 对象：
>>> type(msa)
<class 'Bio.Align.MultipleSeqAlignment'>
>>> print(msa)
Alignment with 4 rows and 7 columns
ACTCCTA seq1
AAT-CTA seq2
CCTACT- seq3
TCTCCTC seq4
>>> alignment = msa.alignment
>>> type(alignment)
<class 'Bio.Align.Alignment'>
>>> print(alignment)
seq1 0 ACTCCTA 7
seq2 0 AAT-CTA 6
seq3 0 CCTACT- 6
seq4 0 TCTCCTC 7
注意到 alignment 属性创建并返回新的 Alignment 与存储在 MultipleSeqAlignment 当时的对象 Alignment 对象被创建。的任何变更 MultipleSeqAlignment 打电话给 alignment 属性不会传播到 Alignment object.但是，您当然可以调用 alignment 再次创建新的房产 Alignment 对象与更新的一致 MultipleSeqAlignment object.
根据多序列比对计算替代矩阵
您可以根据对齐创建自己的替代矩阵。在这个例子中，我们首先从Custalw文件中读取蛋白质序列比对 protein.aln (also在线提供 here )
>>> from Bio import AlignIO
>>> filename = "protein.aln"
>>> msa = AlignIO.read(filename, "clustal")
部分 ClustalW 包含有关此操作的更多信息。
的 substitutions 对齐的属性存储不同残基相互替代的次数：
>>> observed_frequencies = msa.substitutions
为了使示例更具可读性，我们将仅选择具有带电荷侧链的氨基酸：
>>> observed_frequencies = observed_frequencies.select("DEHKR")
>>> print(observed_frequencies)
D E H K R
D 2360.0 255.5 7.5 0.5 25.0
E 255.5 3305.0 16.5 27.0 2.0
H 7.5 16.5 1235.0 16.0 8.5
K 0.5 27.0 16.0 3218.0 116.5
R 25.0 2.0 8.5 116.5 2079.0
从矩阵中去除其他氨基酸的柱和柱。
接下来，我们对矩阵进行规范化：
>>> import numpy as np
>>> observed_frequencies /= np.sum(observed_frequencies)
行或列的总和给出了每个残基的相对出现频率：
>>> residue_frequencies = np.sum(observed_frequencies, 0)
>>> print(residue_frequencies.format("%.4f"))
D 0.2015
E 0.2743
H 0.0976
K 0.2569
R 0.1697
>>> sum(residue_frequencies) == 1.0
True
那么，残基对的预期频率是
>>> expected_frequencies = np.dot(
... residue_frequencies[:, None], residue_frequencies[None, :]
... )
>>> print(expected_frequencies.format("%.4f"))
D E H K R
D 0.0406 0.0553 0.0197 0.0518 0.0342
E 0.0553 0.0752 0.0268 0.0705 0.0465
H 0.0197 0.0268 0.0095 0.0251 0.0166
K 0.0518 0.0705 0.0251 0.0660 0.0436
R 0.0342 0.0465 0.0166 0.0436 0.0288
在这里， residue_frequencies[:, None] 创建由单个列组成的2D数组，值为 residue_frequencies ，而且 residue_frequencies[None, :] 将这些值作为一行的2D数组。取其点积（内部积）创建预期频率矩阵，其中每个条目由两个组成 residue_frequencies 价值观相互叠加。例如， expected_frequencies['D', 'E'] 等于 residue_frequencies['D'] * residue_frequencies['E'] .
我们现在可以通过将观察到的频率除以预期频率并取对数来计算对比矩阵：
>>> m = np.log2(observed_frequencies / expected_frequencies)
>>> print(m)
D E H K R
D 2.1 -1.5 -5.1 -10.4 -4.2
E -1.5 1.7 -4.4 -5.1 -8.3
H -5.1 -4.4 3.3 -4.4 -4.7
K -10.4 -5.1 -4.4 1.9 -2.3
R -4.2 -8.3 -4.7 -2.3 2.5
当执行比对时，该矩阵可以用作替代矩阵。例如，
>>> from Bio.Align import PairwiseAligner
>>> aligner = PairwiseAligner()
>>> aligner.substitution_matrix = m
>>> aligner.gap_score = -3.0
>>> alignments = aligner.align("DEHEK", "DHHKK")
>>> print(alignments[0])
target 0 DEHEK 5
0 |.|.| 5
query 0 DHHKK 5
>>> print("%.2f" % alignments.score)
-2.18
>>> score = m["D", "D"] + m["E", "H"] + m["H", "H"] + m["E", "K"] + m["K", "K"]
>>> print("%.2f" % score)
-2.18
对准工具
有 lots 有多种用于序列比对的算法，包括成对比对和多序列比对。这些计算相对较慢，您通常不想用Python编写这样的算法。对于成对对齐，您可以使用Biopython的 PairwiseAligner (see章 双序列比对 ），它是用C语言实现的，因此速度很快。或者，您可以通过从Python中调用外部对齐程序来运行该程序。通常您会：
准备未对齐序列的输入文件，通常这是一个FASTA文件，您可以使用它创建
Bio.SeqIO(see章 序列输入/输出 ).通过使用Python的运行其命令来运行对齐程序
subprocessmodule.读取工具的输出，即对齐的序列，通常使用
Bio.AlignIO(see本章前面）。
在这里，我们将展示此工作流程的一些示例。
ClustalW
CustalW是一款用于多序列比对的流行命令行工具（还有一个名为CustalX的图形界面）。在尝试在Python中使用ClassalW之前，您应该首先尝试在命令行中亲自手动运行ClassalW工具，以熟悉其他选项。
对于最基本的用途，您需要的只是拥有一个FASTA输入文件，例如 opuntia.fasta （在线或在Biopython源代码的Doc/示例子目录中提供）。这是一个小型FASTA文件，包含七个花椒DNA序列（来自仙人掌科 Opuntia ).默认情况下，ClassalW将根据输入的FASTA文件生成一个对齐和引导树文件，在本例中，该文件的名称是基于输入的FASTA文件 opuntia.aln 和 opuntia.dnd ，但您可以覆盖此内容或使其显式：
>>> import subprocess
>>> cmd = "clustalw2 -infile=opuntia.fasta"
>>> results = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, text=True)
请注意，我们将可执行文件名称指定为 clustalw2 ，表明我们安装了版本二，该版本的文件名与版本一不同 (clustalw ，默认）。幸运的是，两个版本在命令行支持相同的参数集（事实上，功能上应该相同）。
您可能会发现，即使您安装了ClassalW，上述命令也不起作用-您可能会收到有关“命令未找到”的消息（尤其是在Windows上）。这表明ClassalW可执行文件不在PATH（环境变量、要搜索的目录列表）上。您可以更新PATH设置以包括您的ClassalW工具副本的位置（如何执行此操作将取决于您的操作系统），或者简单地输入工具的完整路径。记住，在Python字符串中 \n 和 \t 默认情况下被解释为新行和制表符-这就是为什么我们在未以这种方式翻译的原始字符串的开头放置字母“r”。在指定Windows风格文件名时，这通常是一个好做法。
>>> import os
>>> clustalw_exe = r"C:\Program Files\new clustal\clustalw2.exe"
>>> assert os.path.isfile(clustalw_exe), "Clustal W executable missing"
>>> cmd = clustalw_exe + " -infile=opuntia.fasta"
>>> results = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, text=True)
现在，了解命令行工具如何“工作”很有帮助。当您在命令行上运行工具时，它通常会将文本输出直接打印到屏幕上。可以通过两个“管道”捕获或重定向此文本，称为标准输出（正常结果）和标准错误（用于错误消息和调试消息）。还有标准输入，即输入到工具中的任何文本。这些名称被缩短为stdin、stdout和stderr。当工具完成时，它会有一个返回代码（一个整数），按照惯例，该代码为零表示成功，而非零返回代码则表示发生了错误。
在上面的ClassalW示例中，当在命令行运行时，所有重要的输出都直接写入输出文件。等待时通常打印到屏幕上的所有内容都被捕获在 results.stdout 和 results.stderr ，而返回代码存储在 results.returncode .
我们关心的是两个输出文件：对齐和引导树。我们没有告诉ClassalW要使用哪些文件名，但它默认会根据输入文件选择名称。在这种情况下，输出应该在文件中 opuntia.aln .您应该能够了解如何使用读取对齐方式 Bio.AlignIO 现在：
>>> from Bio import AlignIO
>>> align = AlignIO.read("opuntia.aln", "clustal")
>>> print(align)
Alignment with 7 rows and 906 columns
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273285|gb|AF191659.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273284|gb|AF191658.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273287|gb|AF191661.1|AF191
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273286|gb|AF191660.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273290|gb|AF191664.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273289|gb|AF191663.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273291|gb|AF191665.1|AF191
如果您感兴趣（这是本章主旨的一个例外）， opuntia.dnd ClassalW创建的文件只是一个标准的Newick树文件，并且 Bio.Phylo 可以解析这些：
>>> from Bio import Phylo
>>> tree = Phylo.read("opuntia.dnd", "newick")
>>> Phylo.draw_ascii(tree)
_______________ gi|6273291|gb|AF191665.1|AF191665
__________________________|
| | ______ gi|6273290|gb|AF191664.1|AF191664
| |__|
| |_____ gi|6273289|gb|AF191663.1|AF191663
|
_|_________________ gi|6273287|gb|AF191661.1|AF191661
|
|__________ gi|6273286|gb|AF191660.1|AF191660
|
| __ gi|6273285|gb|AF191659.1|AF191659
|___|
| gi|6273284|gb|AF191658.1|AF191658
章 系统发生学与Bio.Phylo 更深入地涵盖了Biopython对系统发生树的支持。
MUSCLE
MUSCLE是比CustalW更新的多序列比对工具。与之前一样，我们建议您尝试从命令行使用MUSCLE，然后再尝试从Python运行它。
对于最基本的用途，您需要的只是拥有一个FASTA输入文件，例如 opuntia.fasta （在线或在Biopython源代码的Doc/示例子目录中提供）。然后，您可以告诉MUSCLE读取该FASTA文件，并将对齐写入名为 opuntia.txt :
>>> import subprocess
>>> cmd = "muscle -align opuntia.fasta -output opuntia.txt"
>>> results = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, text=True)
MUSCLE将将比对输出为FASTA文件（使用间隔序列）。的 Bio.AlignIO 模块能够使用读取此对齐 format="fasta" :
>>> from Bio import AlignIO
>>> align = AlignIO.read("opuntia.txt", "fasta")
>>> print(align)
Alignment with 7 rows and 906 columns
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273289|gb|AF191663.1|AF191663
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273291|gb|AF191665.1|AF191665
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273290|gb|AF191664.1|AF191664
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273287|gb|AF191661.1|AF191661
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273286|gb|AF191660.1|AF191660
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273285|gb|AF191659.1|AF191659
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273284|gb|AF191658.1|AF191658
您还可以设置其他可选参数;有关详细信息，请参阅MUSCLE的内置帮助。
注射针头和水
的 EMBOSS 套件包括 water 和 needle Smith-Waterman算法局部对齐和Needleman-Wunsch全局对齐的工具。这些工具共享相同的风格界面，因此在两者之间切换很简单-我们只使用 needle 这里.
假设您想在两个序列之间进行全局成对比对，以FASTA格式准备，如下所示：
>HBA_HUMAN
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR
在文件中 alpha.faa ，其次是在文件中 beta.faa :
>HBB_HUMAN
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPK
VKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFG
KEFTPPVQAAYQKVVAGVANALAHKYH
您可以在 Doc/examples/ 目录.
使用将这两个序列相互对齐的命令 needle 如下：
needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5
为什么不尝试在命令提示符下手动运行此命令呢？您应该看到它进行成对比较并将输出记录在文件中 needle.txt (in默认的CLARSS对齐文件格式）。
即使您安装了CLARSS，运行此命令也可能不起作用-您可能会收到有关“命令未找到”的消息（尤其是在Windows上）。这可能意味着CLARSS工具不在PATH环境变量上。您可以更新PATH设置，也可以简单地使用工具的完整路径，例如：
C:\EMBOSS\needle.exe -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5
接下来我们要使用Python为我们运行这个命令。如上所述，为了完全控制，我们建议您使用Python的内置 subprocess 模块：
>>> import sys
>>> import subprocess
>>> cmd = "needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5"
>>> results = subprocess.run(
... cmd,
... stdout=subprocess.PIPE,
... stderr=subprocess.PIPE,
... text=True,
... shell=(sys, platform != "win32"),
... )
>>> print(results.stdout)
>>> print(results.stderr)
Needleman-Wunsch global alignment of two sequences
接下来我们可以用以下方式加载输出文件 Bio.AlignIO 正如本章前面讨论的那样，作为 emboss 格式：
>>> from Bio import AlignIO
>>> align = AlignIO.read("needle.txt", "emboss")
>>> print(align)
Alignment with 2 rows and 149 columns
MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY...KYR HBA_HUMAN
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRF...KYH HBB_HUMAN
在本例中，我们告诉CLARSS将输出写入文件，但您 can 告诉它将输出写入stdout（如果您不想摆脱临时输出文件，这很有用- use outfile=stdout 论点）：
>>> cmd = "needle -outfile=stdout -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5"
>>> child = subprocess.Popen(
... cmd,
... stdout=subprocess.PIPE,
... stderr=subprocess.PIPE,
... text=True,
... shell=(sys.platform != "win32"),
... )
>>> align = AlignIO.read(child.stdout, "emboss")
>>> print(align)
Alignment with 2 rows and 149 columns
MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY...KYR HBA_HUMAN
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRF...KYH HBB_HUMAN
同样，可以阅读 one 来自stdin的输入（例如 asequence="stdin" ).
这只是触及了您可以做的事情的表面 needle 和 water .一个有用的技巧是第二个文件可以包含多个序列（比如五个），然后CLARSS将进行五次成对比对。