多个序列比对对象

本章介绍了较旧的 MultipleSeqAlignment 类和 Bio.AlignIO 中的解析器，这些解析器解析序列比对软件的输出，生成 MultipleSeqAlignment 对象。所谓多个序列比对，是指将多个序列比对在一起的集合，通常会插入间隙字符，并添加前导或尾随间隙，使得所有序列字符串的长度相同。这种比对可以被看作一个字母矩阵，其中每一行在内部都作为 SeqRecord 对象存储。

我们将介绍 MultipleSeqAlignment 对象，它包含这种类型的数据，以及 Bio.AlignIO 模块，用于以各种文件格式读取和写入它们（遵循上一章中 Bio.SeqIO 模块的设计）。请注意，Bio.SeqIO 和 Bio.AlignIO 都可以读取和写入序列比对文件。你应该根据你想要对数据执行的操作来选择合适的工具。

本章的最后一部分将介绍如何使用 Python 中常见的多个序列比对工具（如 ClustalW 和 MUSCLE），以及如何使用 Biopython 解析结果。

解析或读取序列比对

我们有两个用于读取序列比对的函数，Bio.AlignIO.read() 和 Bio.AlignIO.parse()，它们遵循 Bio.SeqIO 中介绍的约定，分别用于包含一个或多个比对的文件。

使用 Bio.AlignIO.parse() 将返回一个迭代器，它提供 MultipleSeqAlignment 对象。迭代器通常用于 for 循环中。在你拥有多个不同比对的情况（例如 PHYLIP 工具 seqboot 的重采样比对，或 EMBOSS 工具 water 或 needle 的多个成对比对，或 Bill Pearson 的 FASTA 工具）下，可以使用迭代器。

然而，在许多情况下，你将处理只包含一个比对的文件。在这种情况下，你应该使用 Bio.AlignIO.read() 函数，它返回一个单独的 MultipleSeqAlignment 对象。

这两个函数都期望两个必填参数

1. 第一个参数是用于读取数据的句柄，通常是打开的文件（见第 什么鬼句柄？ 节），或文件名。

2. 第二个参数是一个小写字符串，用于指定比对格式。与 Bio.SeqIO 中一样，我们不会尝试为你猜测文件格式！请查看 https://biopython.pythonlang.cn/wiki/AlignIO 以获取支持格式的完整列表。

还有一个可选的 seq_count 参数，它将在下面的 模糊比对 节中讨论，用于处理可能包含多个比对的模糊文件格式。

单个比对

例如，以下是在 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 (PF05371) PFAM 条目的种子比对，从现在已经过时的 PFAM 版本中下载的 https://pfam.xfam.org/。我们可以通过以下方式加载该文件（假设它已经保存到磁盘，名为当前工作目录中的 “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 函数，在比对对象上以特定的文件格式显示它，有关详细信息，请参阅 将你的比对对象作为格式化的字符串获取 节。

你有没有注意到在上面的原始文件中，几个序列包含指向 PDB 的数据库交叉引用以及相关的已知二级结构？请尝试以下操作

>>> 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 上提供了友好的 Web 界面，它实际上允许你以几种其他格式下载该比对。以下是该文件在 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 文件格式中。

请注意，你可以使用 Bio.AlignIO 将原始 Stockholm 格式文件转换为 FASTA 文件，而不是使用 Sanger 网站（见下文）。

对于任何支持的文件格式，你都可以通过更改格式字符串来以完全相同的方式加载比对。例如，对 PHYLIP 文件使用 “phylip”，对 NEXUS 文件使用 “nexus”，或对 EMBOSS 工具输出的比对使用 “emboss”。维基页面 (https://biopython.pythonlang.cn/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]

模糊比对

许多比对文件格式可以显式地存储多个比对，并且每个比对之间的分隔是明确的。然而，当使用一般的序列文件格式时，没有这样的块结构。最常见的情况是，比对已保存为 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 可以应对最常见的情况，即所有比对都具有相同数量的记录。这种情况的一个例子是成对比对的集合，它们可以通过 EMBOSS 工具 needle 和 water 生成，尽管在这种情况下，Bio.AlignIO 应该能够理解它们的原生输出，使用 “emboss” 作为格式字符串。

要将这些 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() 函数返回写入文件的对齐方式数量。

注意 - 如果你告诉 Bio.AlignIO.write() 函数写入已经存在的文件，旧文件将被覆盖，没有任何警告。

在序列对齐方式文件格式之间转换

使用 Bio.AlignIO 在序列对齐方式文件格式之间转换的工作方式与使用 Bio.SeqIO 在序列文件格式之间转换的方式相同（第 在序列文件格式之间转换节）。我们通常使用 Bio.AlignIO.parse() 加载对齐方式，然后使用 Bio.AlignIO.write() 保存它们——或者只使用 Bio.AlignIO.convert() 辅助函数。

对于此示例，我们将加载之前使用的 PFAM/Stockholm 格式文件并将其保存为 Clustal 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")

无论哪种方式，你最终都应该得到同一个新的 Clustal 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

一般来说，由于标识符限制，使用严格 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 使用 StringIO 句柄调用 Bio.AlignIO.write()。如果你正在使用旧版本的 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()（行数）和迭代（每一行为一个 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

由于对的排序是任意的，因此计数在对角线以上和以下平均分配。例如，A 与 C 的 9 次比对存储在位置 ['A', 'C'] 处为 4.5，存储在位置 ['C', 'A'] 处为 4.5。这种安排有助于在从这些计数计算替换矩阵时使数学更简单，如第 替换矩阵 节所述。

请注意，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 的比对对象，例如使用 Bio.AlignIO.read(...) 读取，如第 多序列比对对象 章所述。要获取一个将计算摘要信息的的对象，我们只需

>>> from Bio.Align import AlignInfo
>>> summary_align = AlignInfo.SummaryInfo(msa)

summary_align 对象非常有用，它将为你完成以下几件事

1. 计算一个快速共识序列 - 请参见第 计算快速共识序列 节

2. 获取比对的特定位置评分矩阵 - 请参见第 特定位置评分矩阵 节

3. 计算比对的信息内容 - 请参见第 信息内容 节

4. 生成有关比对中替换的信息 - 第 替换矩阵 节详细介绍了如何使用它来生成替换矩阵。

计算快速共识序列

第 计算摘要信息 节中描述的 SummaryInfo 对象提供计算比对快速共识的功能。假设我们有一个名为 summary_align 的 SummaryInfo 对象，我们可以通过以下方式计算共识

>>> consensus = summary_align.dumb_consensus()
>>> consensus
Seq('XCTXCTX')

顾名思义，这是一个非常简单的共识计算器，它只会将共识中每个点的所有残基加起来，如果最常见的残基值高于某个阈值，就会将该残基添加到共识中。如果它没有达到阈值，它会将一个模糊字符添加到共识中。返回的共识对象是一个 Seq 对象。

你可以通过传递可选参数来调整 dumb_consensus 的工作方式

阈值

这是指定在某个位置添加残基之前，该残基必须在该位置的频率的阈值。默认值为 \(0.7\)（表示 \(70\%\)）。

模糊字符

这是要使用的模糊字符。默认值为 ’N’。

或者，你可以将多序列比对对象 msa 转换为新的 Alignment 对象（见 比对对象 节），方法是使用 alignment 属性（见 获取新的比对对象 节）

>>> alignment = msa.alignment

然后你可以创建一个 Motif 对象（见 基序对象 节）

>>> from Bio.motifs import Motif
>>> motif = Motif("ACGT", alignment)

并获得一个快速共识序列

>>> motif.consensus
Seq('ACTCCTA')

motif.counts.calculate_consensus 方法（见 获取共识序列 节）允许你详细指定如何计算共识序列。例如，

>>> motif.counts.calculate_consensus(identity=0.7)
'NCTNCTN'

特定位置评分矩阵

特定位置评分矩阵 (PSSM) 以不同于共识的方式总结比对信息，并且可能对不同的任务有用。基本上，PSSM 是一个计数矩阵。对于比对中的每一列，计算每个字母表字母的数量并将其相加。总数相对于左侧轴上的某个代表性序列显示。该序列可能是共识序列，但也可能是比对中的任何序列。

例如，对于上面的比对

>>> print(msa)
Alignment with 4 rows and 7 columns
ACTCCTA seq1
AAT-CTA seq2
CCTACT- seq3
TCTCCTC seq4

我们可以使用以下命令获得一个带有共识序列的 PSSM

>>> my_pssm = summary_align.pos_specific_score_matrix(consensus, chars_to_ignore=["N"])
>>> print(my_pssm)
    A   C   T
X  2.0 1.0 1.0
C  1.0 3.0 0.0
T  0.0 0.0 4.0
X  1.0 2.0 0.0
C  0.0 4.0 0.0
T  0.0 0.0 4.0
X  2.0 1.0 0.0

其中我们在计算 PSSM 时忽略任何 N 模糊残基。

关于这一点需要注意两点

1. 为了保持字母表的严格性，你只能在 PSSM 的顶部包含存在于比对对象字母表中的字符。PSSM 的顶部轴上不包含间隙。

2. 传递到显示在左侧轴上的序列不必是共识。例如，如果你想在该轴上显示比对中的第二条序列，你需要执行以下操作

   >>> second_seq = msa[1]
   >>> my_pssm = summary_align.pos_specific_score_matrix(second_seq, chars_to_ignore=["N"])
   >>> print(my_pssm)
       A   C   T
   A  2.0 1.0 1.0
   A  1.0 3.0 0.0
   T  0.0 0.0 4.0
   -  1.0 2.0 0.0
   C  0.0 4.0 0.0
   T  0.0 0.0 4.0
   A  2.0 1.0 0.0
   

上面的命令返回一个 PSSM 对象。你可以通过下标访问 PSSM 中的任何元素，例如 your_pssm[sequence_number][residue_count_name]。例如，要获取上面 PSSM 中第二个元素的 ’A’ 残基的计数，你可以执行以下操作

>>> print(my_pssm[5]["T"])
4.0

PSSM 类结构希望使访问元素和美化打印矩阵变得容易。

或者，你可以将多序列比对对象 msa 转换为新的 Alignment 对象（见 比对对象 节），方法是使用 alignment 属性（见 获取新的比对对象 节）

>>> alignment = msa.alignment

然后你可以创建一个 Motif 对象（见 基序对象 节）

>>> from Bio.motifs import Motif
>>> motif = Motif("ACGT", alignment)

并获取每个位置上每个核苷酸的计数

>>> counts = motif.counts
>>> print(counts)
        0      1      2      3      4      5      6
A:   2.00   1.00   0.00   1.00   0.00   0.00   2.00
C:   1.00   3.00   0.00   2.00   4.00   0.00   1.00
G:   0.00   0.00   0.00   0.00   0.00   0.00   0.00
T:   1.00   0.00   4.00   0.00   0.00   4.00   0.00

>>> print(counts["T"][5])
4.0

信息内容

一个潜在的有用的进化保守性度量是序列的信息内容。

针对分子生物学家的信息论的实用介绍可以在 http://www.lecb.ncifcrf.gov/~toms/paper/primer/ 找到。出于我们的目的，我们将关注共识序列或共识序列的一部分的信息内容。我们使用以下公式计算多序列比对中特定列的信息内容

\[IC_{j} = \sum_{i=1}^{N_{a}} P_{ij} \mathrm{log}\left(\frac{P_{ij}}{Q_{i}}\right)\]

其中

• \(IC_{j}\) – 比对中第 \(j\) 列的信息内容。

• \(N_{a}\) – 字母表中字母的数量。

• \(P_{ij}\) – 第 \(j\) 列中特定字母 \(i\) 的频率（即，如果 G 在比对列中出现了 6 次中的 3 次，这将是 0.5）

• \(Q_{i}\) – 字母 \(i\) 的预期频率。这是一个可选参数，使用方式由用户自行决定。默认情况下，它会自动分配给蛋白质字母表中的 \(0.05 = 1/20\)，分配给核酸字母表中的 \(0.25 = 1/4\)。这是为了在没有先验分布假设的情况下获取信息内容。在假设先验或使用非标准字母表时，你应该提供 \(Q_{i}\) 的值。

现在我们对 Biopython 中计算的信息内容有了一定的了解，让我们看看如何获取比对特定区域的信息内容。

首先，我们需要使用比对获取比对摘要对象，我们假设它名为 summary_align（见第 计算摘要信息 节），获取该对象的步骤如下。获取此对象后，计算区域的信息内容非常简单

>>> e_freq_table = {"A": 0.3, "G": 0.2, "T": 0.3, "C": 0.2}
>>> info_content = summary_align.information_content(
...     2, 6, e_freq_table=e_freq_table, chars_to_ignore=["N"]
... )
>>> info_content  
6.3910647...

现在，info_content 将包含区域 [2:6] 中相对于预期频率的相对信息内容。

返回的值使用 2 为底的对数在上面的公式中计算。你可以通过传递参数 log_base 作为你想要使用的底来修改它

>>> info_content = summary_align.information_content(
...     2, 6, e_freq_table=e_freq_table, log_base=10, chars_to_ignore=["N"]
... )
>>> info_content  
1.923902...

默认情况下，在计算列的相对信息列时，不会考虑在列中频率为 0 的核苷酸或氨基酸残基。如果这不是想要的结果，你可以使用 pseudo_count 代替。

>>> info_content = summary_align.information_content(
...     2, 6, e_freq_table=e_freq_table, chars_to_ignore=["N", "-"], pseudo_count=1
... )
>>> info_content  
4.299651...

在这种情况下，特定字母 \(i\) 在第 \(j\) 列中的观察频率 \(P_{ij}\) 计算如下

\[P_{ij} = \frac{n_{ij} + k\times Q_{i}}{N_{j} + k}\]

其中

• \(k\) – 你作为参数传递的伪计数。

• \(k\) – 你作为参数传递的伪计数。

• \(Q_{i}\) – 如上所述的字母 \(i\) 的预期频率。

现在，您已准备好计算信息内容。如果您想尝试将它应用到一些现实问题中，最好深入研究信息内容的文献，了解它的使用方法。希望您的探索不会发现对编码此函数造成的任何错误！

获取新式 Alignment 对象

使用 alignment 属性从旧式 MultipleSeqAlignment 对象创建一个新式 Alignment 对象（参见第 Alignment 对象 节）。

>>> 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 对象。但是，您可以再次调用 alignment 属性来创建一个新的 Alignment 对象，该对象与更新后的 MultipleSeqAlignment 对象一致。

从多序列比对计算替换矩阵

您可以从比对创建自己的替换矩阵。在这个例子中，我们首先从 Clustalw 文件 protein.aln（也可在线获取 此处）读取蛋白质序列比对。

>>> 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] 创建一个二维数组，该数组包含一个列，其中包含 residue_frequencies 的值，而 residue_frequencies[None, :] 创建一个二维数组，该数组将这些值作为单行。取它们的点积（内积）将创建一个预期频率矩阵，其中每个条目包含两个 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

比对工具

有很多算法可以用于比对序列，包括成对比对和多序列比对。这些计算相对较慢，一般来说您不会希望在 Python 中编写这样的算法。对于成对比对，您可以使用 Biopython 的 PairwiseAligner（参见第 成对序列比对 章），它是在 C 中实现的，因此速度很快。或者，您可以通过从 Python 调用外部比对程序来运行它。通常，您会

1. 准备一个未比对序列的输入文件，通常这将是一个 FASTA 文件，您可以使用 Bio.SeqIO 创建它（参见第 序列输入/输出 章）。

2. 使用 Python 的 subprocess 模块运行比对程序的命令。

3. 读取工具的输出，即您的比对序列，通常使用 Bio.AlignIO（参见本章前面）。

这里，我们将展示此工作流程的几个示例。

ClustalW

ClustalW 是一个流行的命令行多序列比对工具（还有一个名为 ClustalX 的图形界面）。在尝试从 Python 中使用 ClustalW 之前，您应该先尝试在命令行手动运行 ClustalW 工具，以熟悉其他选项。

对于最基本的使用方法，您只需要有一个 FASTA 输入文件，例如 opuntia.fasta（可在网上或 Biopython 源代码的 Doc/examples 子目录中找到）。这是一个包含七个仙人掌 DNA 序列（来自仙人掌科 Opuntia）的小型 FASTA 文件。默认情况下，ClustalW 将生成一个比对和引导树文件，其名称基于输入 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，默认值）具有不同的文件名。幸运的是，这两个版本都支持命令行中相同的参数集（实际上，应该在功能上是相同的）。

您可能会发现，即使您安装了 ClustalW，上面的命令也不起作用 - 您可能会收到关于“命令未找到”的提示（尤其是在 Windows 上）。这表明 ClustalW 可执行文件不在您的 PATH（环境变量，要搜索的目录列表）中。您可以更新您的 PATH 设置以包含 ClustalW 工具副本的位置（具体操作方法取决于您的操作系统），或者直接输入工具的完整路径。请记住，在 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。当工具完成时，它有一个返回码（一个整数），该码按照惯例为成功时为零，而非零返回码表示发生了错误。

在上面的 ClustalW 示例中，当在命令行运行时，所有重要的输出都直接写入输出文件。您等待时通常打印到屏幕的所有内容都捕获在 results.stdout 和 results.stderr 中，而返回码存储在 results.returncode 中。

我们关心的是两个输出文件，比对和引导树。我们没有告诉 ClustalW 使用什么文件名，但它默认为基于输入文件选择名称。在本例中，输出应该位于文件 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

如果您有兴趣（这是一个与本章主要内容无关的补充），ClustalW 创建的 opuntia.dnd 文件只是一个标准的 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 是一个比 ClustalW 更新的多序列比对工具。如前所述，建议您先尝试从命令行使用 MUSCLE，然后再尝试从 Python 运行它。

对于最基本的使用方法，您只需要有一个 FASTA 输入文件，例如 opuntia.fasta（可在网上或 Biopython 源代码的 Doc/examples 子目录中找到）。然后，您可以告诉 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 needle 和 water

EMBOSS 套件包括 water 和 needle 工具，用于 Smith-Waterman 算法局部比对和 Needleman-Wunsch 全局比对。这些工具共享相同的样式界面，因此在两者之间切换非常容易 - 我们在这里只使用 needle。

假设您要对两个序列进行全局成对比对，这些序列以 FASTA 格式准备如下。

>HBA_HUMAN
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR

在一个文件 alpha.faa 中，另一个在文件 beta.faa 中。

>HBB_HUMAN
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPK
VKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFG
KEFTPPVQAAYQKVVAGVANALAHKYH

您可以在 Biopython 源代码的 Doc/examples/ 目录下找到这些示例文件的副本。

使用 needle 对这两个序列进行比对的命令如下。

needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5

为什么不在命令提示符下尝试手动运行它？您应该看到它执行了成对比较并将输出记录在文件 needle.txt 中（以默认的 EMBOSS 比对文件格式）。

即使您安装了 EMBOSS，运行此命令也可能不起作用 - 您可能会收到关于“命令未找到”的提示（尤其是在 Windows 上）。这很可能是因为 EMBOSS 工具不在您的 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

在这个例子中，我们告诉 EMBOSS 将输出写入文件，但您可以告诉它将输出写入标准输出（如果您不想删除临时输出文件，这很有用 - 使用 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

同样，也可以从标准输入读取一个输入（例如 asequence="stdin"）。

这仅仅触及了您使用 needle 和 water 可以做的事情的表面。一个有用的技巧是，第二个文件可以包含多个序列（例如五个），然后 EMBOSS 将执行五个成对比对。