序列比对

序列比对是两个或多个序列的集合，它们已经彼此比对 - 通常包括插入间隙，以及添加前导或尾随间隙 - 使所有序列字符串具有相同的长度。

比对可以扩展到每个序列的全部长度，也可以限制为每个序列的子部分。在 Biopython 中，所有序列比对都由一个 Alignment 对象表示，在第比对对象节中进行了描述。 Alignment 对象可以通过解析比对软件（如 Clustal 或 BLAT）的输出获得（在第读取和写入比对节中进行了描述），或者通过使用 Biopython 的成对序列比对器获得，该比对器可以将两个序列相互比对（在第成对序列比对章中进行了描述）。

有关旧的 MultipleSeqAlignment 类的描述，以及 Bio.AlignIO 中解析序列比对软件输出的解析器，这些解析器会生成 MultipleSeqAlignment 对象，请参阅第多序列比对对象章。

比对对象

Alignment 类在 Bio.Align 中定义。通常，你会通过解析比对程序的输出（第读取和写入比对节）或通过运行 Biopython 的成对比对器（第成对序列比对章）来获得 Alignment 对象。为了便于理解本节内容，我们将从头开始创建一个 Alignment 对象。

从序列和坐标创建比对对象

假设你有三个序列

>>> seqA = "CCGGTTTTT"
>>> seqB = "AGTTTAA"
>>> seqC = "AGGTTT"
>>> sequences = [seqA, seqB, seqC]

为了创建一个 Alignment 对象，我们还需要定义序列如何彼此比对的坐标。为此，我们使用 NumPy 数组

>>> import numpy as np
>>> coordinates = np.array([[1, 3, 4, 7, 9], [0, 2, 2, 5, 5], [0, 2, 3, 6, 6]])

这些坐标定义了以下序列段的比对

SeqA[1:3]、SeqB[0:2] 和 SeqC[0:2] 彼此比对；
SeqA[3:4] 和 SeqC[2:3] 彼此比对，seqB 中有一个核苷酸的间隙；
SeqA[4:7]、SeqB[2:5] 和 SeqC[3:6] 彼此比对；
SeqA[7:9] 不与 seqB 或 seqC 比对。

请注意，比对不包括 seqA 的第一个核苷酸和 seqB 的最后两个核苷酸。

现在我们可以创建 Alignment 对象

>>> from Bio.Align import Alignment
>>> alignment = Alignment(sequences, coordinates)
>>> alignment  
<Alignment object (3 rows x 8 columns) at ...>

比对对象有一个属性 sequences，指向包含在此比对中的序列

>>> alignment.sequences
['CCGGTTTTT', 'AGTTTAA', 'AGGTTT']

以及一个具有比对坐标的属性 coordinates

>>> alignment.coordinates
array([[1, 3, 4, 7, 9],
       [0, 2, 2, 5, 5],
       [0, 2, 3, 6, 6]])

打印 Alignment 对象以明确显示比对

>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6

每个序列的起始和结束坐标分别显示在比对的左侧和右侧。

从已比对的序列创建比对对象

如果你从已比对的序列开始，用连字符表示间隙，那么可以使用 parse_printed_alignment 类方法来计算坐标。此方法主要用于 Biopython 的比对解析器（见第读取和写入比对节），但它可能对其他目的也有用。例如，你可以从已比对的序列构建 Alignment 对象，如下所示

>>> lines = ["CGGTTTTT", "AG-TTT--", "AGGTTT--"]
>>> for line in lines:
...     print(line)
...
CGGTTTTT
AG-TTT--
AGGTTT--
>>> lines = [line.encode() for line in lines]  # convert to bytes
>>> lines
[b'CGGTTTTT', b'AG-TTT--', b'AGGTTT--']
>>> sequences, coordinates = Alignment.parse_printed_alignment(lines)
>>> sequences
[b'CGGTTTTT', b'AGTTT', b'AGGTTT']
>>> sequences = [sequence.decode() for sequence in sequences]
>>> sequences
['CGGTTTTT', 'AGTTT', 'AGGTTT']
>>> print(coordinates)
[[0 2 3 6 8]
 [0 2 2 5 5]
 [0 2 3 6 6]]

seqA 的初始 G 核苷酸和 seqB 的最终 CC 核苷酸没有包含在比对中，因此在这里缺失。但这个问题很容易解决

>>> from Bio.Seq import Seq
>>> sequences[0] = "C" + sequences[0]
>>> sequences[1] = sequences[1] + "AA"
>>> sequences
['CCGGTTTTT', 'AGTTTAA', 'AGGTTT']
>>> coordinates[0, :] += 1
>>> print(coordinates)
[[1 3 4 7 9]
 [0 2 2 5 5]
 [0 2 3 6 6]]

现在我们可以创建 Alignment 对象

>>> alignment = Alignment(sequences, coordinates)
>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6

它与我们在第从序列和坐标创建比对对象节中创建的 Alignment 对象相同。

默认情况下，Alignment 初始化程序的 coordinates 参数为 None，这意味着假设比对中没有间隙。无间隙比对中的所有序列必须具有相同的长度。如果 coordinates 参数为 None，则初始化程序将为你填充 Alignment 对象的 coordinates 属性

>>> ungapped_alignment = Alignment(["ACGTACGT", "AAGTACGT", "ACGTACCT"])
>>> ungapped_alignment  
<Alignment object (3 rows x 8 columns) at ...>
>>> print(ungapped_alignment.coordinates)
[[0 8]
 [0 8]
 [0 8]]
>>> print(ungapped_alignment)
                  0 ACGTACGT 8
                  0 AAGTACGT 8
                  0 ACGTACCT 8

常见的比对属性

以下属性通常存在于 Alignment 对象中

sequences：这是一个包含相互比对的序列的列表。根据比对的创建方式，序列可以具有以下类型
- 普通 Python 字符串；
- Seq;
- MutableSeq;
- SeqRecord;
- 字节;
- 字节数组;
- 数据类型为 numpy.int32 的 NumPy 数组；
- 任何其他具有连续缓冲区格式为 "c"、"B"、"i" 或 "I" 的对象；
- 在创建比对的 PairwiseAligner 对象的 alphabet 属性中定义的对象的列表或元组（见第广义成对比对节）。
对于成对比对（意味着两个序列的比对），属性 target 和 query 分别是 sequences[0] 和 sequences[1] 的别名。
coordinates：一个包含整数的 NumPy 数组，存储定义序列如何彼此比对的序列索引；
score：比对得分，由解析器在比对文件中找到，或者由 PairwiseAligner 计算（见第基本用法节）；
annotations：一个字典，存储与比对相关的绝大多数其他注释；
column_annotations: 用于存储沿对齐方向扩展且长度与对齐长度相同的注释的字典，例如一致序列（有关示例，请参见第 ClustalW 节）。

由 Bio.Align 中的解析器创建的 Alignment 对象可能具有其他属性，具体取决于读取对齐的比对文件格式。

对齐的切片和索引

形式为 alignment[k, i:j] 的切片（其中 k 是一个整数，而 i 和 j 是整数或不存在）返回一个字符串，该字符串显示目标（如果 k=0）或查询（如果 k=1）的对齐序列（包括间隙），其中仅包含打印对齐中第 i 列到第 j 列。

为了说明这一点，在以下示例中，打印的对齐有 8 列

>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6

>>> alignment.length
8

要分别获取对齐的序列字符串，请使用

>>> alignment[0]
'CGGTTTTT'
>>> alignment[1]
'AG-TTT--'
>>> alignment[2]
'AGGTTT--'
>>> alignment[0, :]
'CGGTTTTT'
>>> alignment[1, :]
'AG-TTT--'
>>> alignment[0, 1:-1]
'GGTTTT'
>>> alignment[1, 1:-1]
'G-TTT-'

也可以使用对整数的迭代来选择要包含的列

>>> alignment[0, (1, 2, 4)]
'GGT'
>>> alignment[1, range(0, 5, 2)]
'A-T'

要获取打印对齐中位置 [i, j] 处的字母，请使用 alignment[i, j]；如果在该位置找到间隙，则这将返回 "-"

>>> alignment[0, 2]
'G'
>>> alignment[2, 6]
'-'

要获取对齐中的特定列，请使用

>>> alignment[:, 0]
'CAA'
>>> alignment[:, 1]
'GGG'
>>> alignment[:, 2]
'G-G'

形式为 alignment[i:j:k] 的切片返回一个新的 Alignment 对象，其中仅包含对齐的序列 [i:j:k]

>>> alignment[1:]  
<Alignment object (2 rows x 6 columns) at ...>
>>> print(alignment[1:])
target            0 AG-TTT 5
                  0 ||-||| 6
query             0 AGGTTT 6

形式为 alignment[:, i:j] 的切片（其中 i 和 j 是整数或不存在）返回一个新的 Alignment 对象，其中仅包含打印对齐中第 i 列到第 j 列。

为上述示例对齐提取前 4 列，得到

>>> alignment[:, :4]  
<Alignment object (3 rows x 4 columns) at ...>
>>> print(alignment[:, :4])
                  1 CGGT 5
                  0 AG-T 3
                  0 AGGT 4

类似地，提取最后 6 列得到

>>> alignment[:, -6:]  
<Alignment object (3 rows x 6 columns) at ...>
>>> print(alignment[:, -6:])
                  3 GTTTTT 9
                  2 -TTT-- 5
                  2 GTTT-- 6

列索引也可以是整数的迭代

>>> print(alignment[:, (1, 3, 0)])
                  0 GTC 3
                  0 GTA 3
                  0 GTA 3

调用 alignment[:, :] 返回对齐的副本。

获取有关对齐的信息

对齐形状

对齐序列的数量由 len(alignment) 返回

>>> len(alignment)
3

对齐长度定义为打印对齐中列的数量。这等于匹配次数、不匹配次数以及每个序列中间隙总长度的总和

>>> alignment.length
8

shape 属性返回一个元组，该元组包含对齐的长度和打印对齐中列的数量

>>> alignment.shape
(3, 8)

比较对齐

如果 alignment1.sequences 和 alignment2.sequences 中的每个序列都彼此相等，并且 alignment1.coordinates 和 alignment2.coordinates 包含相同的坐标，则两个对齐彼此相等（意味着 alignment1 == alignment2 评估为 True）。如果这两个条件中的任何一个不满足，则 alignment1 == alignment2 评估为 False。两个对齐的不等（例如 alignment1 < alignment2）通过首先比较 alignment1.sequences 和 alignment2.sequences 来确定，如果它们相等，则通过比较 alignment1.coordinates 与 alignment2.coordinates 来确定。

查找对齐序列的索引

对于成对对齐，对齐的 aligned 属性返回目标和查询序列中彼此对齐的子序列的起始和结束索引。如果目标 (t) 和查询 (q) 之间的对齐包含 \(N\) 个块，您将获得两个长度为 \(N\) 的元组

(((t_start1, t_end1), (t_start2, t_end2), ..., (t_startN, t_endN)),
 ((q_start1, q_end1), (q_start2, q_end2), ..., (q_startN, q_endN)))

例如，

>>> pairwise_alignment = alignment[:2, :]
>>> print(pairwise_alignment)
target            1 CGGTTTTT 9
                  0 .|-|||-- 8
query             0 AG-TTT-- 5

>>> print(pairwise_alignment.aligned)
[[[1 3]
  [4 7]]

 [[0 2]
  [2 5]]]

请注意，不同的对齐可能具有彼此对齐的相同子序列。特别是，如果对齐仅在间隙放置方面彼此不同，则可能会发生这种情况

>>> pairwise_alignment1 = Alignment(["AAACAAA", "AAAGAAA"],
...                                 np.array([[0, 3, 4, 4, 7], [0, 3, 3, 4, 7]]))  # fmt: skip
...
>>> pairwise_alignment2 = Alignment(["AAACAAA", "AAAGAAA"],
...                                 np.array([[0, 3, 3, 4, 7], [0, 3, 4, 4, 7]]))  # fmt: skip
...
>>> print(pairwise_alignment1)
target            0 AAAC-AAA 7
                  0 |||--||| 8
query             0 AAA-GAAA 7

>>> print(pairwise_alignment2)
target            0 AAA-CAAA 7
                  0 |||--||| 8
query             0 AAAG-AAA 7

>>> pairwise_alignment1.aligned
array([[[0, 3],
        [4, 7]],

       [[0, 3],
        [4, 7]]])
>>> pairwise_alignment2.aligned
array([[[0, 3],
        [4, 7]],

       [[0, 3],
        [4, 7]]])

indices 属性返回一个二维 NumPy 数组，其中包含对齐中每个字母的序列索引，间隙由 -1 表示

>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6

>>> alignment.indices
array([[ 1,  2,  3,  4,  5,  6,  7,  8],
       [ 0,  1, -1,  2,  3,  4, -1, -1],
       [ 0,  1,  2,  3,  4,  5, -1, -1]])

inverse_indices 属性返回一个一维 NumPy 数组列表，每个对齐序列一个，其中包含序列中每个字母的对齐中的列索引。对齐中未包含的字母由 -1 表示

>>> alignment.sequences
['CCGGTTTTT', 'AGTTTAA', 'AGGTTT']
>>> alignment.inverse_indices  
[array([-1,  0,  1,  2,  3,  4,  5,  6,  7]),
 array([ 0,  1,  3,  4,  5, -1, -1]),
 array([0, 1, 2, 3, 4, 5])]

计算同一性、不匹配和间隙

counts 方法计算成对对齐的同一性、不匹配和间隙的数量。对于两个以上序列的对齐，计算所有序列对中的同一性、不匹配和间隙的数量，并将它们求和。这三个数字作为 AlignmentCounts 对象返回，它是一个具有字段 gaps、identities 和 mismatches 的 namedtuple。此方法目前不接受任何参数，但在将来可能会修改为接受可选参数，以便可以自定义其行为。

>>> print(pairwise_alignment)
target            1 CGGTTTTT 9
                  0 .|-|||-- 8
query             0 AG-TTT-- 5

>>> pairwise_alignment.counts()
AlignmentCounts(gaps=3, identities=4, mismatches=1)
>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6

>>> alignment.counts()
AlignmentCounts(gaps=8, identities=14, mismatches=2)

字母频率

frequencies 方法计算每个字母在对齐的每一列中出现的频率

>>> alignment.frequencies  
{'C': array([1., 0., 0., 0., 0., 0., 0., 0.]),
 'G': array([0., 3., 2., 0., 0., 0., 0., 0.]),
 'T': array([0., 0., 0., 3., 3., 3., 1., 1.]),
 'A': array([2., 0., 0., 0., 0., 0., 0., 0.]),
 '-': array([0., 0., 1., 0., 0., 0., 2., 2.])}

替换

使用 substitutions 方法查找每对核苷酸之间的替换次数

>>> m = alignment.substitutions
>>> print(m)
    A   C   G   T
A 1.0 0.0 0.0 0.0
C 2.0 0.0 0.0 0.0
G 0.0 0.0 4.0 0.0
T 0.0 0.0 0.0 9.0

请注意，矩阵不是对称的：行字母 R 和列字母 C 的计数是在序列中字母 R 被出现在其下方的序列中的字母 C 替换的次数。例如，对齐到后面序列中的 A 的 C 的数量是

>>> m["C", "A"]
2.0

而对齐到后面序列中的 C 的 A 的数量是

>>> m["A", "C"]
0.0

要获取对称矩阵，请使用

>>> m += m.transpose()
>>> m /= 2.0
>>> print(m)
    A   C   G   T
A 1.0 1.0 0.0 0.0
C 1.0 0.0 0.0 0.0
G 0.0 0.0 4.0 0.0
T 0.0 0.0 0.0 9.0

>>> m["A", "C"]
1.0
>>> m["C", "A"]
1.0

对齐中 A 和 T 之间的总替换次数为 1.0 + 1.0 = 2。

对齐作为数组

使用 NumPy，您可以将 alignment 对象转换为字母数组。特别是，这可能对对齐内容的快速计算很有用。

>>> align_array = np.array(alignment)
>>> align_array.shape
(3, 8)
>>> align_array  
array([[b'C', b'G', b'G', b'T', b'T', b'T', b'T', b'T'],
       [b'A', b'G', b'-', b'T', b'T', b'T', b'-', b'-'],
       [b'A', b'G', b'G', b'T', b'T', b'T', b'-', b'-']], dtype='|S1')

默认情况下，这将为您提供一个 bytes 字符数组（数据类型为 dtype='|S1'）。您可以通过使用 dtype='U' 来创建一个 Unicode（Python 字符串）字符数组

>>> align_array = np.array(alignment, dtype="U")

>>> align_array  
array([['C', 'G', 'G', 'T', 'T', 'T', 'T', 'T'],
       ['A', 'G', '-', 'T', 'T', 'T', '-', '-'],
       ['A', 'G', 'G', 'T', 'T', 'T', '-', '-']], dtype='<U1')

（打印的 dtype 将是 ‘<U1’ 或 ‘>U1’，具体取决于您的系统是 little-endian 还是 big-endian）。请注意，alignment 对象和 NumPy 数组 align_array 是内存中独立的对象 - 编辑其中一个不会更新另一个！

对齐操作

对齐排序

sort 方法对对齐序列进行排序。默认情况下，排序是基于每个序列的 id 属性（如果可用）或序列内容进行的。

>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6

>>> alignment.sort()
>>> print(alignment)
                  0 AGGTTT-- 6
                  0 AG-TTT-- 5
                  1 CGGTTTTT 9

或者，您可以提供一个 key 函数来确定排序顺序。例如，您可以按 GC 含量递增对序列进行排序

>>> from Bio.SeqUtils import gc_fraction
>>> alignment.sort(key=gc_fraction)
>>> print(alignment)
                  0 AG-TTT-- 5
                  0 AGGTTT-- 6
                  1 CGGTTTTT 9

请注意，key 函数应用于完整序列（包括 seqB 的初始 A 和最终 GG 核苷酸），而不仅仅应用于对齐部分。

reverse 参数允许您反转排序顺序以获得按 GC 含量递减排列的序列

>>> alignment.sort(key=gc_fraction, reverse=True)
>>> print(alignment)
                  1 CGGTTTTT 9
                  0 AGGTTT-- 6
                  0 AG-TTT-- 5

反向互补对齐

反向互补对齐将对每个序列执行反向互补，并重新计算坐标

>>> alignment.sequences
['CCGGTTTTT', 'AGGTTT', 'AGTTTAA']
>>> rc_alignment = alignment.reverse_complement()
>>> print(rc_alignment.sequences)
['AAAAACCGG', 'AAACCT', 'TTAAACT']
>>> print(rc_alignment)
                  0 AAAAACCG 8
                  0 --AAACCT 6
                  2 --AAA-CT 7

>>> alignment[:, :4].sequences
['CCGGTTTTT', 'AGGTTT', 'AGTTTAA']
>>> print(alignment[:, :4])
                  1 CGGT 5
                  0 AGGT 4
                  0 AG-T 3

>>> rc_alignment = alignment[:, :4].reverse_complement()
>>> rc_alignment[:, :4].sequences
['AAAAACCGG', 'AAACCT', 'TTAAACT']
>>> print(rc_alignment[:, :4])
                  4 ACCG 8
                  2 ACCT 6
                  4 A-CT 7

反向互补对齐保留其列注释（反向顺序），但丢弃所有其他注释。

添加对齐

如果对齐具有相同数量的行，则可以将它们加在一起形成扩展的对齐。例如，让我们先创建两个对齐

>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> a1 = SeqRecord(Seq("AAAAC"), id="Alpha")
>>> b1 = SeqRecord(Seq("AAAC"), id="Beta")
>>> c1 = SeqRecord(Seq("AAAAG"), id="Gamma")
>>> a2 = SeqRecord(Seq("GTT"), id="Alpha")
>>> b2 = SeqRecord(Seq("TT"), id="Beta")
>>> c2 = SeqRecord(Seq("GT"), id="Gamma")
>>> left = Alignment(
...     [a1, b1, c1], coordinates=np.array([[0, 3, 4, 5], [0, 3, 3, 4], [0, 3, 4, 5]])
... )
>>> left.annotations = {"tool": "demo", "name": "start"}
>>> left.column_annotations = {"stats": "CCCXC"}
>>> right = Alignment(
...     [a2, b2, c2], coordinates=np.array([[0, 1, 2, 3], [0, 0, 1, 2], [0, 1, 1, 2]])
... )
>>> right.annotations = {"tool": "demo", "name": "end"}
>>> right.column_annotations = {"stats": "CXC"}

现在，让我们看一下这两个对齐

>>> print(left)
Alpha             0 AAAAC 5
Beta              0 AAA-C 4
Gamma             0 AAAAG 5

>>> print(right)
Alpha             0 GTT 3
Beta              0 -TT 2
Gamma             0 G-T 2

添加两个对齐将按行将两个对齐组合起来

>>> combined = left + right
>>> print(combined)
Alpha             0 AAAACGTT 8
Beta              0 AAA-C-TT 6
Gamma             0 AAAAGG-T 7

为了使此方法有效，两种比对必须具有相同数量的序列（这里它们都具有 3 行）。

>>> len(left)
3
>>> len(right)
3
>>> len(combined)
3

这些序列是 SeqRecord 对象，可以将它们加在一起。有关如何处理注释的详细信息，请参阅第序列注释对象章。此示例是一个特殊情况，因为两个原始比对都共享相同的名称，这意味着在添加行时，它们也会获得相同的名称。

任何公共注释都会被保留，但不同的注释将会丢失。这与 SeqRecord 注释中使用的行为相同，旨在防止意外传播不适当的值。

>>> combined.annotations
{'tool': 'demo'}

同样，任何公共的每列注释也会被合并。

>>> combined.column_annotations
{'stats': 'CCCXCCXC'}

映射成对的序列比对 

假设您有一个转录本与染色体的成对比对。

>>> chromosome = "AAAAAAAACCCCCCCAAAAAAAAAAAGGGGGGAAAAAAAA"
>>> transcript = "CCCCCCCGGGGGG"
>>> sequences1 = [chromosome, transcript]
>>> coordinates1 = np.array([[8, 15, 26, 32], [0, 7, 7, 13]])
>>> alignment1 = Alignment(sequences1, coordinates1)
>>> print(alignment1)
target            8 CCCCCCCAAAAAAAAAAAGGGGGG 32
                  0 |||||||-----------|||||| 24
query             0 CCCCCCC-----------GGGGGG 13

以及转录本与序列之间的成对比对（例如，通过 RNA 测序获得）。

>>> rnaseq = "CCCCGGGG"
>>> sequences2 = [transcript, rnaseq]
>>> coordinates2 = np.array([[3, 11], [0, 8]])
>>> alignment2 = Alignment(sequences2, coordinates2)
>>> print(alignment2)
target            3 CCCCGGGG 11
                  0 ||||||||  8
query             0 CCCCGGGG  8

在 alignment1 上使用 map 方法，并以 alignment2 作为参数，来查找 RNA 序列与基因组的比对。

>>> alignment3 = alignment1.map(alignment2)
>>> print(alignment3)
target           11 CCCCAAAAAAAAAAAGGGG 30
                  0 ||||-----------|||| 19
query             0 CCCC-----------GGGG  8

>>> print(alignment3.coordinates)
[[11 15 26 30]
 [ 0  4  4  8]]
>>> format(alignment3, "psl")
'8\t0\t0\t0\t0\t0\t1\t11\t+\tquery\t8\t0\t8\ttarget\t40\t11\t30\t2\t4,4,\t0,4,\t11,26,\n'

为了能够打印序列，在本示例中，我们使用具有定义的序列内容的序列构建了 alignment1 和 alignment2。但是，映射比对并不依赖于序列内容；仅使用 alignment1 和 alignment2 的坐标来构建 alignment3 的坐标。

map 方法还可以用于将不同基因组组装之间的比对提升。在这种情况下，self 是两个基因组组装之间的 DNA 比对，参数是转录本与其中一个基因组组装的比对。

>>> from Bio import Align
>>> chain = Align.read("Blat/panTro5ToPanTro6.over.chain", "chain")
>>> chain.sequences[0].id
'chr1'
>>> len(chain.sequences[0].seq)
228573443
>>> chain.sequences[1].id
'chr1'
>>> len(chain.sequences[1].seq)
224244399
>>> import numpy as np
>>> np.set_printoptions(threshold=5)  # print 5 array elements per row
>>> print(chain.coordinates)  
[[122250000 122250400 122250400 ... 122909818 122909819 122909835]
 [111776384 111776784 111776785 ... 112019962 112019962 112019978]]

表明黑猩猩基因组组装 panTro5 上的 chr1 的 122250000:122909835 范围对应于黑猩猩基因组组装 panTro6 上的 chr1 的 111776384:112019978 范围。有关链文件格式的更多信息，请参阅部分 UCSC 链文件格式。

>>> transcript = Align.read("Blat/est.panTro5.psl", "psl")
>>> transcript.sequences[0].id
'chr1'
>>> len(transcript.sequences[0].seq)
228573443
>>> transcript.sequences[1].id
'DC525629'
>>> len(transcript.sequences[1].seq)
407
>>> print(transcript.coordinates)
[[122835789 122835847 122840993 122841145 122907212 122907314]
 [       32        90        90       242       242       344]]

这表明表达序列标签 DC525629 的 32:344 核苷酸范围对应于黑猩猩基因组组装 panTro5 上的 chr1 的 122835789:122907314 范围。请注意，目标序列 chain.sequences[0].seq 和目标序列 transcript.sequences[0] 具有相同的长度。

>>> len(chain.sequences[0].seq) == len(transcript.sequences[0].seq)
True

我们交换链的目标和查询，以便 chain 的查询对应于 transcript 的目标。

>>> chain = chain[::-1]
>>> chain.sequences[0].id
'chr1'
>>> len(chain.sequences[0].seq)
224244399
>>> chain.sequences[1].id
'chr1'
>>> len(chain.sequences[1].seq)
228573443
>>> print(chain.coordinates)  
[[111776384 111776784 111776785 ... 112019962 112019962 112019978]
 [122250000 122250400 122250400 ... 122909818 122909819 122909835]]
>>> np.set_printoptions(threshold=1000)  # reset the print options

现在，我们可以通过调用 chain.map（以 transcript 作为参数）来获得 DC525629 相对于黑猩猩基因组组装 panTro6 的坐标。

>>> lifted_transcript = chain.map(transcript)
>>> lifted_transcript.sequences[0].id
'chr1'
>>> len(lifted_transcript.sequences[0].seq)
224244399
>>> lifted_transcript.sequences[1].id
'DC525629'
>>> len(lifted_transcript.sequences[1].seq)
407
>>> print(lifted_transcript.coordinates)
[[111982717 111982775 111987921 111988073 112009200 112009302]
 [       32        90        90       242       242       344]]

这表明表达序列标签 DC525629 的 32:344 核苷酸范围对应于黑猩猩基因组组装 panTro6 上的 chr1 的 111982717:112009302 范围。请注意，DC525629 在黑猩猩基因组组装 panTro5 上的基因组跨度为 122907314 - 122835789 = 71525 bp，而在 panTro6 上的基因组跨度为 112009302 - 111982717 = 26585 bp。

映射多序列比对 

考虑黑猩猩、人类、猕猴、狨猴、小鼠和大鼠的基因组序列的多序列比对。

>>> from Bio import Align
>>> path = "Blat/panTro5.maf"
>>> genome_alignment = Align.read(path, "maf")
>>> for record in genome_alignment.sequences:
...     print(record.id, len(record.seq))
...
panTro5.chr1 228573443
hg19.chr1 249250621
rheMac8.chr1 225584828
calJac3.chr18 47448759
mm10.chr3 160039680
rn6.chr2 266435125
>>> print(genome_alignment.coordinates)
[[133922962 133922962 133922970 133922970 133922972 133922972 133922995
  133922998 133923010]
 [155784573 155784573 155784581 155784581 155784583 155784583 155784606
  155784609 155784621]
 [130383910 130383910 130383918 130383918 130383920 130383920 130383943
  130383946 130383958]
 [  9790455   9790455   9790463   9790463   9790465   9790465   9790488
    9790491   9790503]
 [ 88858039  88858036  88858028  88858026  88858024  88858020  88857997
   88857997  88857985]
 [188162970 188162967 188162959 188162959 188162957 188162953 188162930
  188162930 188162918]]
>>> print(genome_alignment)
panTro5.c 133922962 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAACTTAAAggcc
hg19.chr1 155784573 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAACTTAAAggcc
rheMac8.c 130383910 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAACTTAAAggcc
calJac3.c   9790455 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAGCTTAAAggct
mm10.chr3  88858039 TATAATAATTGTATATGTCACAGAAAAAAATGAATTTTCAAT---GACTTAATAGCC
rn6.chr2  188162970 TACAATAATTG--TATGTCATAGAAAAAAATGAATTTTCAAT---AACTTAATAGCC

panTro5.c 133923010
hg19.chr1 155784621
rheMac8.c 130383958
calJac3.c   9790503
mm10.chr3  88857985
rn6.chr2  188162918

假设我们想用较新版本的基因组组装（panTro5、hg19、rheMac8、calJac3、mm10 和 rn6）替换旧版本（panTro6、hg38、rheMac10、calJac4、mm39 和 rn7）。为此，我们需要每个物种的旧组装版本与新组装版本之间的成对比对。这些由 UCSC 作为链文件提供，通常用于 UCSC 的 liftOver 工具。Biopython 源代码分发中的 Tests/Align 子目录中的 .chain 文件是从 UCSC 的 .chain 文件中提取的，仅包含相关的基因组区域。例如，要将 panTro5 提升到 panTro6，我们使用 panTro5ToPanTro6.chain 文件，其内容如下：

chain 1198066 chr1 228573443 + 133919957 133932620 chr1 224244399 + 130607995 130620657 1
4990    0   2
1362    3   0
6308

为了提升每个物种的基因组组装，我们读取相应的 .chain 文件。

>>> paths = [
...     "Blat/panTro5ToPanTro6.chain",
...     "Blat/hg19ToHg38.chain",
...     "Blat/rheMac8ToRheMac10.chain",
...     "Blat/calJac3ToCalJac4.chain",
...     "Blat/mm10ToMm39.chain",
...     "Blat/rn6ToRn7.chain",
... ]
>>> liftover_alignments = [Align.read(path, "chain") for path in paths]
>>> for liftover_alignment in liftover_alignments:
...     print(liftover_alignment.target.id, liftover_alignment.coordinates[0, :])
...
chr1 [133919957 133924947 133924947 133926309 133926312 133932620]
chr1 [155184381 156354347 156354348 157128497 157128497 157137496]
chr1 [130382477 130383872 130383872 130384222 130384222 130388520]
chr18 [9786631 9787941 9788508 9788508 9795062 9795065 9795737]
chr3 [66807541 74196805 74196831 94707528 94707528 94708176 94708178 94708718]
chr2 [188111581 188158351 188158351 188171225 188171225 188228261 188228261
 188236997]

请注意，物种的顺序在 liftover_alignments 和 genome_alignment.sequences 中相同。现在，我们可以将多序列比对提升到新的基因组组装版本。

>>> genome_alignment = genome_alignment.mapall(liftover_alignments)
>>> for record in genome_alignment.sequences:
...     print(record.id, len(record.seq))
...
chr1 224244399
chr1 248956422
chr1 223616942
chr18 47031477
chr3 159745316
chr2 249053267
>>> print(genome_alignment.coordinates)
[[130611000 130611000 130611008 130611008 130611010 130611010 130611033
  130611036 130611048]
 [155814782 155814782 155814790 155814790 155814792 155814792 155814815
  155814818 155814830]
 [ 95186253  95186253  95186245  95186245  95186243  95186243  95186220
   95186217  95186205]
 [  9758318   9758318   9758326   9758326   9758328   9758328   9758351
    9758354   9758366]
 [ 88765346  88765343  88765335  88765333  88765331  88765327  88765304
   88765304  88765292]
 [174256702 174256699 174256691 174256691 174256689 174256685 174256662
  174256662 174256650]]

由于 .chain 文件不包含序列内容，因此我们无法直接打印序列比对。相反，我们单独读取基因组序列（作为 .2bit 文件，因为它允许延迟加载；请参阅部分序列文件作为字典），用于每个物种。

>>> from Bio import SeqIO
>>> names = ("panTro6", "hg38", "rheMac10", "calJac4", "mm39", "rn7")
>>> for i, name in enumerate(names):
...     filename = f"{name}.2bit"
...     genome = SeqIO.parse(filename, "twobit")
...     chromosome = genome_alignment.sequences[i].id
...     assert len(genome_alignment.sequences[i]) == len(genome[chromosome])
...     genome_alignment.sequences[i] = genome[chromosome]
...     genome_alignment.sequences[i].id = f"{name}.{chromosome}"
...
>>> print(genome_alignment)
panTro6.c 130611000 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAACTTAAAggcc
hg38.chr1 155814782 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAACTTAAAggcc
rheMac10.  95186253 ---ACTAGTTA--CA----GTAACAGAAAATAAAATTTAAATAGAAACTTAAAggcc
calJac4.c   9758318 ---ACTAGTTA--CA----GTAACAGAaaataaaatttaaatagaagcttaaaggct
mm39.chr3  88765346 TATAATAATTGTATATGTCACAGAAAAAAATGAATTTTCAAT---GACTTAATAGCC
rn7.chr2  174256702 TACAATAATTG--TATGTCATAGAAAAAAATGAATTTTCAAT---AACTTAATAGCC

panTro6.c 130611048
hg38.chr1 155814830
rheMac10.  95186205
calJac4.c   9758366
mm39.chr3  88765292
rn7.chr2  174256650

mapall 方法还可以用于从相应氨基酸序列的多序列比对创建密码子序列的多序列比对（有关详细信息，请参阅部分生成密码子序列的多序列比对）。

Alignments 类 

Alignments（复数）类继承自 AlignmentsAbstractBaseClass 和 list，可以用作列表来存储 Alignment 对象。Alignments 对象的行为与 list 对象的行为在两个重要方面有所不同。

Alignments 对象是它自己的迭代器，与 Bio.Align.parse 返回的迭代器一致（请参阅部分读取比对）或成对比对程序返回的迭代器（请参阅部分成对序列比对）。在迭代器上调用 iter 将始终返回 Alignments 对象本身。相反，在 list 对象上调用 iter 会每次创建一个新的迭代器，允许您为给定的列表创建多个独立的迭代器。

在本示例中，alignment_iterator1 和 alignment_iterator2 从列表中获取，并彼此独立地运行。

>>> alignment_list = [alignment1, alignment2, alignment3]
>>> alignment_iterator1 = iter(alignment_list)
>>> alignment_iterator2 = iter(alignment_list)
>>> next(alignment_iterator1)  
<Alignment object (2 rows x 24 columns) at ...>
>>> next(alignment_iterator2)  
<Alignment object (2 rows x 24 columns) at ...>
>>> next(alignment_iterator1)  
<Alignment object (2 rows x 8 columns) at ...>
>>> next(alignment_iterator1)  
<Alignment object (2 rows x 19 columns) at ...>
>>> next(alignment_iterator2)  
<Alignment object (2 rows x 8 columns) at ...>
>>> next(alignment_iterator2)  
<Alignment object (2 rows x 19 columns) at ...>

相反，通过在 Alignments 对象上调用 iter 获取的 alignment_iterator1 和 alignment_iterator2 彼此相同。

>>> from Bio.Align import Alignments
>>> alignments = Alignments([alignment1, alignment2, alignment3])
>>> alignment_iterator1 = iter(alignments)
>>> alignment_iterator2 = iter(alignments)
>>> alignment_iterator1 is alignment_iterator2
True
>>> next(alignment_iterator1)  
<Alignment object (2 rows x 24 columns) at ...>
>>> next(alignment_iterator2)  
<Alignment object (2 rows x 8 columns) at ...>
>>> next(alignment_iterator1)  
<Alignment object (2 rows x 19 columns) at ...>
>>> next(alignment_iterator2)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration

在 Alignments 对象上调用 iter 会将迭代器重置为其第一个项目，因此您可以再次循环遍历它。您还可以使用 for 循环多次迭代比对，该循环隐式地对迭代器调用 iter。

>>> for item in alignments:
...     print(repr(item))  
...
<Alignment object (2 rows x 24 columns) at ...>
<Alignment object (2 rows x 8 columns) at ...>
<Alignment object (2 rows x 19 columns) at ...>

>>> for item in alignments:
...     print(repr(item))  
...
<Alignment object (2 rows x 24 columns) at ...>
<Alignment object (2 rows x 8 columns) at ...>
<Alignment object (2 rows x 19 columns) at ...>

此行为与常规 Python 列表一致，以及与 Bio.Align.parse 返回的迭代器一致（请参阅部分读取比对）或成对比对程序返回的迭代器一致（请参阅部分成对序列比对）。

元数据可以作为属性存储在 Alignments 对象上，而普通 list 不接受属性。

>>> alignment_list.score = 100  
Traceback (most recent call last):
 ...
AttributeError: 'list' object has no attribute 'score'...
>>> alignments.score = 100
>>> alignments.score
100

读取和写入比对 

Bio.Align.read 和 Bio.Align.parse 函数可以将 Clustal 等序列比对软件的输出解析为 Alignment 对象。它们的使用类似于 Bio.SeqIO 中的 read 和 parse 函数（请参阅部分解析或读取序列）：read 函数用于读取包含单个比对的输出文件，并返回 Alignment 对象，而 parse 函数返回一个迭代器来迭代存储在包含一个或多个比对的输出文件中的比对。部分比对文件格式描述了可以在 Bio.Align 中解析的比对格式。Bio.Align 还提供了一个 write 函数，可以以这些格式中的大多数格式写入比对。

读取比对 

使用 Bio.Align.parse 解析序列比对文件。例如，文件 ucsc_mm9_chr10.maf 包含 48 个以 MAF（多序列比对格式）格式存储的多序列比对（请参阅部分多序列比对格式 (MAF)）。

>>> from Bio import Align
>>> alignments = Align.parse("MAF/ucsc_mm9_chr10.maf", "maf")
>>> alignments  
<Bio.Align.maf.AlignmentIterator object at 0x...>

其中 "maf" 是文件格式。由 Bio.Align.parse 返回的对齐对象可能包含存储在文件中找到的元数据的属性，例如用于创建对齐的软件的版本号。每个文件格式存储的特定属性在第对齐文件格式节中描述。对于 MAF 文件，我们可以获得文件格式版本和使用的评分方案

>>> alignments.metadata
{'MAF Version': '1', 'Scoring': 'autoMZ.v1'}

由于对齐文件可能非常大，Align.parse 返回对齐的迭代器，因此您不必同时在内存中存储所有对齐。您可以遍历这些对齐并打印出，例如，每个对齐中对齐序列的数量

>>> for a in alignments:
...     print(len(a.sequences))  
...
2
4
5
6
...
15
14
7
6

您还可以对对齐调用 len 以获取对齐的数量。

>>> len(alignments)
48

根据文件格式，对齐的数量可能明确存储在文件中（例如在 bigBed、bigPsl 和 bigMaf 文件的情况下），否则对齐的数量通过循环遍历它们一次（并将迭代器返回到其原始位置）来计算。如果文件很大，因此 len 返回可能需要相当长的时间。但是，由于对齐的数量已缓存，因此后续对 len 的调用将快速返回。

如果对齐的数量不是过大并且可以放入内存，您可以将对齐迭代器转换为对齐列表。为此，您可以对 alignments 调用 list

>>> alignment_list = list(alignments)
>>> len(alignment_list)
48
>>> alignment_list[27]  
<Alignment object (3 rows x 91 columns) at 0x...>
>>> print(alignment_list[27])
mm9.chr10   3019377 CCCCAGCATTCTGGCAGACACAGTG-AAAAGAGACAGATGGTCACTAATAAAATCTGT-A
felCat3.s     46845 CCCAAGTGTTCTGATAGCTAATGTGAAAAAGAAGCATGTGCCCACCAGTAAGCTTTGTGG
canFam2.c  47545247 CCCAAGTGTTCTGATTGCCTCTGTGAAAAAGAAACATGGGCCCGCTAATAagatttgcaa

mm9.chr10   3019435 TAAATTAG-ATCTCAGAGGATGGATGGACCA  3019465
felCat3.s     46785 TGAACTAGAATCTCAGAGGATG---GGACTC    46757
canFam2.c  47545187 tgacctagaatctcagaggatg---ggactc 47545159

但这将丢失元数据信息

>>> alignment_list.metadata  
Traceback (most recent call last):
  ...
AttributeError: 'list' object has no attribute 'metadata'

相反，您可以请求对齐的完整切片

>>> type(alignments)
<class 'Bio.Align.maf.AlignmentIterator'>
>>> alignments = alignments[:]
>>> type(alignments)
<class 'Bio.Align.Alignments'>

这将返回一个 Bio.Align.Alignments 对象，该对象可以作为列表使用，同时保留元数据信息

>>> len(alignments)
48
>>> print(alignments[11])
mm9.chr10   3014742 AAGTTCCCTCCATAATTCCTTCCTCCCACCCCCACA 3014778
calJac1.C      6283 AAATGTA-----TGATCTCCCCATCCTGCCCTG---    6311
otoGar1.s    175262 AGATTTC-----TGATGCCCTCACCCCCTCCGTGCA  175231
loxAfr1.s      9317 AGGCTTA-----TG----CCACCCCCCACCCCCACA    9290

>>> alignments.metadata
{'MAF Version': '1', 'Scoring': 'autoMZ.v1'}

写入对齐

要将对齐写入文件，请使用

>>> from Bio import Align
>>> target = "myfile.txt"
>>> Align.write(alignments, target, "clustal")

其中 alignments 是单个对齐或对齐列表，target 是文件名或打开的文件类对象，而 "clustal" 是要使用的文件格式。由于某些文件格式允许或要求将元数据与对齐一起存储，因此您可能希望使用 Alignments（复数）类而不是简单的对齐列表（参见第 Alignments 类节），允许您将元数据字典存储为 alignments 对象的属性

>>> from Bio import Align
>>> alignments = Align.Alignments(alignments)
>>> metadata = {"Program": "Biopython", "Version": "1.81"}
>>> alignments.metadata = metadata
>>> target = "myfile.txt"
>>> Align.write(alignments, target, "clustal")

打印对齐

对于文本（非二进制）格式，您可以对对齐调用 Python 的内置 format 函数以获取显示以请求格式显示的对齐的字符串，或者在格式化 (f-) 字符串中使用 Alignment 对象。如果在没有参数的情况下调用，则 format 函数将返回对齐的字符串表示形式

>>> str(alignment)
'                  1 CGGTTTTT 9\n                  0 AGGTTT-- 6\n                  0 AG-TTT-- 5\n'
>>> format(alignment)
'                  1 CGGTTTTT 9\n                  0 AGGTTT-- 6\n                  0 AG-TTT-- 5\n'
>>> print(format(alignment))
                  1 CGGTTTTT 9
                  0 AGGTTT-- 6
                  0 AG-TTT-- 5

通过指定第对齐文件格式节中所示的格式之一，format 将创建一个字符串，显示以请求格式显示的对齐

>>> format(alignment, "clustal")
'sequence_0                          CGGTTTTT\nsequence_1                          AGGTTT--\nsequence_2                          AG-TTT--\n\n\n'
>>> print(format(alignment, "clustal"))
sequence_0                          CGGTTTTT
sequence_1                          AGGTTT--
sequence_2                          AG-TTT--



>>> print(f"*** this is the alignment in Clustal format: ***\n{alignment:clustal}\n***")
*** this is the alignment in Clustal format: ***
sequence_0                          CGGTTTTT
sequence_1                          AGGTTT--
sequence_2                          AG-TTT--



***
>>> format(alignment, "maf")
'a\ns sequence_0 1 8 + 9 CGGTTTTT\ns sequence_1 0 6 + 6 AGGTTT--\ns sequence_2 0 5 + 7 AG-TTT--\n\n'
>>> print(format(alignment, "maf"))
a
s sequence_0 1 8 + 9 CGGTTTTT
s sequence_1 0 6 + 6 AGGTTT--
s sequence_2 0 5 + 7 AG-TTT--

由于可选关键字参数不能与 Python 的内置 format 函数或格式化字符串一起使用，因此 Alignment 类具有一个 format 方法，该方法带有可选参数来自定义对齐格式，如以下小节所述。例如，我们可以使用特定列数以 BED 格式（参见第浏览器扩展数据 (BED) 节）打印对齐

>>> print(pairwise_alignment)
target            1 CGGTTTTT 9
                  0 .|-|||-- 8
query             0 AG-TTT-- 5

>>> print(format(pairwise_alignment, "bed"))  
target  1   7   query   0   +   1   7   0   2   2,3,    0,3,

>>> print(pairwise_alignment.format("bed"))  
target  1   7   query   0   +   1   7   0   2   2,3,    0,3,

>>> print(pairwise_alignment.format("bed", bedN=3))  
target  1   7

>>> print(pairwise_alignment.format("bed", bedN=6))  
target  1   7   query   0   +

对齐文件格式

下表显示了可以在 Bio.Align 中解析的对齐格式。在 Bio.Align 函数中用于指定文件格式的格式参数 fmt 不区分大小写。大多数这些文件格式也可以由 Bio.Align 写入，如表所示。

文件格式 `fmt`	描述	文本 / 二进制	受 `write` 支持	小节
`a2m`	A2M	文本	是	1.7.11
`bed`	浏览器可扩展数据 (BED)	文本	是	1.7.14
`bigbed`	bigBed	二进制	是	1.7.15
`bigmaf`	bigMaf	二进制	是	1.7.19
`bigpsl`	bigPsl	二进制	是	1.7.17
`chain`	UCSC chain 文件	文本	是	1.7.20
`clustal`	ClustalW	文本	是	1.7.2
`emboss`	EMBOSS	文本	否	1.7.5
`` exonerate``	Exonerate	文本	是	1 .7.7
`fasta`	已比对的 FASTA	文本	是	1.7.1
`hhr`	HH-suite 输出文件	文本	否	1.7.10
`maf`	多序列比对格式 (MAF)	文本	是	1.7.18
`mauve`	Mauve 扩展多 FastA (xmfa) 格式	文本	是	1.7.12
`msf`	GCG 多序列格式 (MSF)	文本	否	1.7.6
`nexus`	NEXUS	文本	是	1.7.8
`phylip`	PHYLIP 输出文件	文本	是	1.7.4
`psl`	模式空间布局 (PSL)	文本	是	1.7.16
`sam`	序列对齐/映射 (SAM)	文本	是	1.7.13
`` stockholm``	Stockholm	文本	是	1 .7.3
`表格`	BLAST 或 FASTA 的表格输出	文本	否	1.7.9

对齐 FASTA

对齐 FASTA 格式的文件存储一个（成对或多个）序列对齐，其中对齐中的间隙由连字符 (-) 表示。使用 fmt="fasta" 读取或写入对齐 FASTA 格式的文件。请注意，这不同于 William Pearson 的 FASTA 对齐程序生成的输出（解析此类输出在第 BLAST 或 FASTA 的表格输出节中描述）。

Biopython 测试套件中的 probcons.fa 文件以对齐 FASTA 格式存储一个多重对齐。该文件的内容如下

>plas_horvu
D-VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG-VD-VSKISQEEYLTAPGETFSVTLTV---PGTYGFYCEPHAGAGMVGKVTV
>plas_chlre
--VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG-VN-ADAISRDDYLNAPGETYSVKLTA---AGEYGYYCEPHQGAGMVGKIIV
>plas_anava
--VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKSADLAKSLSHKQLLMSPGQSTSTTFPADAPAGEYTFYCEPHRGAGMVGKITV
>plas_proho
VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG-ES-APALSNTKLRIAPGSFYSVTLGT---PGTYSFYCTPHRGAGMVGTITV
>azup_achcy
VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG-AE-A-------FKSKINENYKVTFTA---PGVYGVKCTPHYGMGMVGVVEV

要读取此文件，请使用

>>> from Bio import Align
>>> alignment = Align.read("probcons.fa", "fasta")
>>> alignment  
<Alignment object (5 rows x 101 columns) at ...>

我们可以打印对齐以查看其默认表示形式

>>> print(alignment)
plas_horv         0 D-VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG-VD-VSKISQE
plas_chlr         0 --VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG-VN-ADAISRD
plas_anav         0 --VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKSADLAKSLSHK
plas_proh         0 VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG-ES-APALSNT
azup_achc         0 VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG-AE-A------

plas_horv        57 EYLTAPGETFSVTLTV---PGTYGFYCEPHAGAGMVGKVTV 95
plas_chlr        56 DYLNAPGETYSVKLTA---AGEYGYYCEPHQGAGMVGKIIV 94
plas_anav        58 QLLMSPGQSTSTTFPADAPAGEYTFYCEPHRGAGMVGKITV 99
plas_proh        56 KLRIAPGSFYSVTLGT---PGTYSFYCTPHRGAGMVGTITV 94
azup_achc        51 -FKSKINENYKVTFTA---PGVYGVKCTPHYGMGMVGVVEV 88

或者我们可以将其打印为对齐的 FASTA 格式

>>> print(format(alignment, "fasta"))
>plas_horvu
D-VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG-VD-VSKISQEEYLTAPGETFSVTLTV---PGTYGFYCEPHAGAGMVGKVTV
>plas_chlre
--VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG-VN-ADAISRDDYLNAPGETYSVKLTA---AGEYGYYCEPHQGAGMVGKIIV
>plas_anava
--VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKSADLAKSLSHKQLLMSPGQSTSTTFPADAPAGEYTFYCEPHRGAGMVGKITV
>plas_proho
VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG-ES-APALSNTKLRIAPGSFYSVTLGT---PGTYSFYCTPHRGAGMVGTITV
>azup_achcy
VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG-AE-A-------FKSKINENYKVTFTA---PGVYGVKCTPHYGMGMVGVVEV

或任何其他可用格式，例如 Clustal（参见第 ClustalW 节）

>>> print(format(alignment, "clustal"))
plas_horvu                          D-VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG-
plas_chlre                          --VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG-
plas_anava                          --VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKS
plas_proho                          VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG-
azup_achcy                          VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG-

plas_horvu                          VD-VSKISQEEYLTAPGETFSVTLTV---PGTYGFYCEPHAGAGMVGKVT
plas_chlre                          VN-ADAISRDDYLNAPGETYSVKLTA---AGEYGYYCEPHQGAGMVGKII
plas_anava                          ADLAKSLSHKQLLMSPGQSTSTTFPADAPAGEYTFYCEPHRGAGMVGKIT
plas_proho                          ES-APALSNTKLRIAPGSFYSVTLGT---PGTYSFYCTPHRGAGMVGTIT
azup_achcy                          AE-A-------FKSKINENYKVTFTA---PGVYGVKCTPHYGMGMVGVVE

plas_horvu                          V
plas_chlre                          V
plas_anava                          V
plas_proho                          V
azup_achcy                          V

与对齐关联的序列是 SeqRecord 对象

>>> alignment.sequences
[SeqRecord(seq=Seq('DVLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSGVDVSKI...VTV'), id='plas_horvu', name='<unknown name>', description='', dbxrefs=[]), SeqRecord(seq=Seq('VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSGVNADAIS...IIV'), id='plas_chlre', name='<unknown name>', description='', dbxrefs=[]), SeqRecord(seq=Seq('VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKSADLAKS...ITV'), id='plas_anava', name='<unknown name>', description='', dbxrefs=[]), SeqRecord(seq=Seq('VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDKVPAGESAPALS...ITV'), id='plas_proho', name='<unknown name>', description='', dbxrefs=[]), SeqRecord(seq=Seq('VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDKGHNVETIKGMIPDGAEAFKS...VEV'), id='azup_achcy', name='<unknown name>', description='', dbxrefs=[])]

请注意，这些序列不包含间隙（“-” 字符），因为对齐信息存储在 coordinates 属性中

>>> print(alignment.coordinates)
[[ 0  1  1 33 34 42 44 48 48 50 50 51 58 73 73 95]
 [ 0  0  0 32 33 41 43 47 47 49 49 50 57 72 72 94]
 [ 0  0  0 32 33 41 43 47 48 50 51 52 59 74 77 99]
 [ 0  1  2 34 35 43 43 47 47 49 49 50 57 72 72 94]
 [ 0  1  2 34 34 42 44 48 48 50 50 51 51 66 66 88]]

使用 Align.write 将此对齐写入文件（在这里，我们将使用 StringIO 对象而不是文件）

>>> from io import StringIO
>>> stream = StringIO()
>>> Align.write(alignment, stream, "FASTA")
1
>>> print(stream.getvalue())
>plas_horvu
D-VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG-VD-VSKISQEEYLTAPGETFSVTLTV---PGTYGFYCEPHAGAGMVGKVTV
>plas_chlre
--VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG-VN-ADAISRDDYLNAPGETYSVKLTA---AGEYGYYCEPHQGAGMVGKIIV
>plas_anava
--VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKSADLAKSLSHKQLLMSPGQSTSTTFPADAPAGEYTFYCEPHRGAGMVGKITV
>plas_proho
VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG-ES-APALSNTKLRIAPGSFYSVTLGT---PGTYSFYCTPHRGAGMVGTITV
>azup_achcy
VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG-AE-A-------FKSKINENYKVTFTA---PGVYGVKCTPHYGMGMVGVVEV

请注意，Align.write 返回写入的对齐数量（在本例中为 1）。

ClustalW

Clustal 是一套多序列比对程序，既可作为独立程序，也可作为网络服务器使用。文件 opuntia.aln（在线或在 Biopython 源代码的 Doc/examples 子目录中可用）是 Clustal 生成的输出文件。其前几行如下

CLUSTAL 2.1 multiple sequence alignment


gi|6273285|gb|AF191659.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273284|gb|AF191658.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273287|gb|AF191661.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273286|gb|AF191660.1|AF191      TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273290|gb|AF191664.1|AF191      TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273289|gb|AF191663.1|AF191      TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273291|gb|AF191665.1|AF191      TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
                                    ******* **** *************************************

...

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("opuntia.aln", "clustal")

alignments 上的 metadata 属性存储文件标题中显示的信息

>>> alignments.metadata
{'Program': 'CLUSTAL', 'Version': '2.1'}

您可以对 alignments 调用 next 以提取第一个（也是唯一）对齐

>>> alignment = next(alignments)
>>> print(alignment)  
gi|627328         0 TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT
gi|627328         0 TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT
gi|627328         0 TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT
gi|627328         0 TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT
gi|627329         0 TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT
gi|627328         0 TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT
gi|627329         0 TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT

gi|627328        60 CTAAATGATATACGATTCCACTATGTAAGGTCTTTGAATCATATCATAAAAGACAATGTA
gi|627328        60 CTAAATGATATACGATTCCACTATGTAAGGTCTTTGAATCATATCATAAAAGACAATGTA
gi|627328        60 CTAAATGATATACGATTCCACTATGTAAGGTCTTTGAATCATATCATAAAAGACAATGTA
gi|627328        60 CTAAATGATATACGATTCCACTA...

如果您对元数据不感兴趣，那么使用 Align.read 函数会更方便，因为无论如何每个 Clustal 文件都只包含一个对齐

>>> from Bio import Align
>>> alignment = Align.read("opuntia.aln", "clustal")

Clustal 输出文件中的每个对齐块下面的共识行在序列在每个位置保守的情况下包含一个星号。此信息存储在 alignment 的 column_annotations 属性中

>>> alignment.column_annotations  
{'clustal_consensus': '******* **** **********************************...

以 clustal 格式打印 alignment 将显示序列对齐，但不包括元数据

>>> print(format(alignment, "clustal"))  
gi|6273285|gb|AF191659.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273284|gb|AF191658.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273287|gb|AF191661.1|AF191      TATACATT...

以 clustal 格式写入 alignments 将包括元数据和序列对齐

>>> from io import StringIO
>>> stream = StringIO()
>>> Align.write(alignments, stream, "clustal")
1
>>> print(stream.getvalue())  
CLUSTAL 2.1 multiple sequence alignment


gi|6273285|gb|AF191659.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273284|gb|AF191658.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273287|gb|AF191661.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273286|gb|AF191660.1|AF191      TATACATAAAAGAAG...

如果您要手动创建对齐，并且想要在输出中包含元数据信息，请使用 Alignments（复数）对象（参见第 Alignments 类节）。

Stockholm

这是 Stockholm 文件格式中蛋白质序列对齐的示例，该格式由 PFAM 使用

# STOCKHOLM 1.0
#=GF ID   7kD_DNA_binding
#=GF AC   PF02294.20
#=GF DE   7kD DNA-binding domain
#=GF AU   Mian N;0000-0003-4284-4749
#=GF AU   Bateman A;0000-0002-6982-4660
#=GF SE   Pfam-B_8148 (release 5.2)
#=GF GA   25.00 25.00;
#=GF TC   26.60 46.20;
#=GF NC   23.20 19.20;
#=GF BM   hmmbuild HMM.ann SEED.ann
#=GF SM   hmmsearch -Z 57096847 -E 1000 --cpu 4 HMM pfamseq
#=GF TP   Domain
#=GF CL   CL0049
#=GF RN   [1]
#=GF RM   3130377
#=GF RT   Microsequence analysis of DNA-binding proteins 7a, 7b, and 7e
#=GF RT   from the archaebacterium Sulfolobus acidocaldarius.
#=GF RA   Choli T, Wittmann-Liebold B, Reinhardt R;
#=GF RL   J Biol Chem 1988;263:7087-7093.
#=GF DR   INTERPRO; IPR003212;
#=GF DR   SCOP; 1sso; fa;
#=GF DR   SO; 0000417; polypeptide_domain;
#=GF CC   This family contains members of the hyper-thermophilic
#=GF CC   archaebacterium  7kD DNA-binding/endoribonuclease P2 family.
#=GF CC   There are five 7kD DNA-binding proteins, 7a-7e, found as
#=GF CC   monomers in the cell. Protein 7e shows the  tightest DNA-binding
#=GF CC   ability.
#=GF SQ   3
#=GS DN7_METS5/4-61   AC A4YEA2.1
#=GS DN7A_SACS2/3-61  AC P61991.2
#=GS DN7A_SACS2/3-61  DR PDB; 1SSO A; 2-60;
#=GS DN7A_SACS2/3-61  DR PDB; 1JIC A; 2-60;
#=GS DN7A_SACS2/3-61  DR PDB; 2CVR A; 2-60;
#=GS DN7A_SACS2/3-61  DR PDB; 1B4O A; 2-60;
#=GS DN7E_SULAC/3-60  AC P13125.2
DN7_METS5/4-61              KIKFKYKGQDLEVDISKVKKVWKVGKMVSFTYDD.NGKTGRGAVSEKDAPKELLNMIGK
DN7A_SACS2/3-61             TVKFKYKGEEKQVDISKIKKVWRVGKMISFTYDEGGGKTGRGAVSEKDAPKELLQMLEK
#=GR DN7A_SACS2/3-61  SS    EEEEESSSSEEEEETTTEEEEEESSSSEEEEEE-SSSSEEEEEEETTTS-CHHHHHHTT
DN7E_SULAC/3-60             KVRFKYKGEEKEVDTSKIKKVWRVGKMVSFTYDD.NGKTGRGAVSEKDAPKELMDMLAR
#=GC SS_cons                EEEEESSSSEEEEETTTEEEEEESSSSEEEEEE-SSSSEEEEEEETTTS-CHHHHHHTT
#=GC seq_cons               KVKFKYKGEEKEVDISKIKKVWRVGKMVSFTYDD.NGKTGRGAVSEKDAPKELLsMLuK
//

这是 7kD_DNA_binding (PF02294.20) PFAM 条目的种子对齐，从 InterPro 网站 (https://www.ebi.ac.uk/interpro/) 下载。此版本的 PFAM 条目也作为 Biopython 源代码分发中的 pfam2.seed.txt 文件在 Tests/Stockholm/ 子目录中可用。我们可以按如下方式加载此文件

>>> from Bio import Align
>>> alignment = Align.read("pfam2.seed.txt", "stockholm")
>>> alignment  
<Alignment object (3 rows x 59 columns) at ...>

我们可以打印出对齐的摘要

>>> print(alignment)
DN7_METS5         0 KIKFKYKGQDLEVDISKVKKVWKVGKMVSFTYDD-NGKTGRGAVSEKDAPKELLNMIGK
DN7A_SACS         0 TVKFKYKGEEKQVDISKIKKVWRVGKMISFTYDEGGGKTGRGAVSEKDAPKELLQMLEK
DN7E_SULA         0 KVRFKYKGEEKEVDTSKIKKVWRVGKMVSFTYDD-NGKTGRGAVSEKDAPKELMDMLAR

DN7_METS5        58
DN7A_SACS        59
DN7E_SULA        58

您也可以对对齐对象调用 Python 的内置 format 函数以以特定文件格式显示它（有关详细信息，请参见第打印对齐节），例如在 Stockholm 格式中以重新生成文件

>>> print(format(alignment, "stockholm"))
# STOCKHOLM 1.0
#=GF ID   7kD_DNA_binding
#=GF AC   PF02294.20
#=GF DE   7kD DNA-binding domain
#=GF AU   Mian N;0000-0003-4284-4749
#=GF AU   Bateman A;0000-0002-6982-4660
#=GF SE   Pfam-B_8148 (release 5.2)
#=GF GA   25.00 25.00;
#=GF TC   26.60 46.20;
#=GF NC   23.20 19.20;
#=GF BM   hmmbuild HMM.ann SEED.ann
#=GF SM   hmmsearch -Z 57096847 -E 1000 --cpu 4 HMM pfamseq
#=GF TP   Domain
#=GF CL   CL0049
#=GF RN   [1]
#=GF RM   3130377
#=GF RT   Microsequence analysis of DNA-binding proteins 7a, 7b, and 7e from
#=GF RT   the archaebacterium Sulfolobus acidocaldarius.
#=GF RA   Choli T, Wittmann-Liebold B, Reinhardt R;
#=GF RL   J Biol Chem 1988;263:7087-7093.
#=GF DR   INTERPRO; IPR003212;
#=GF DR   SCOP; 1sso; fa;
#=GF DR   SO; 0000417; polypeptide_domain;
#=GF CC   This family contains members of the hyper-thermophilic
#=GF CC   archaebacterium  7kD DNA-binding/endoribonuclease P2 family. There
#=GF CC   are five 7kD DNA-binding proteins, 7a-7e, found as monomers in the
#=GF CC   cell. Protein 7e shows the  tightest DNA-binding ability.
#=GF SQ   3
#=GS DN7_METS5/4-61   AC A4YEA2.1
#=GS DN7A_SACS2/3-61  AC P61991.2
#=GS DN7A_SACS2/3-61  DR PDB; 1SSO A; 2-60;
#=GS DN7A_SACS2/3-61  DR PDB; 1JIC A; 2-60;
#=GS DN7A_SACS2/3-61  DR PDB; 2CVR A; 2-60;
#=GS DN7A_SACS2/3-61  DR PDB; 1B4O A; 2-60;
#=GS DN7E_SULAC/3-60  AC P13125.2
DN7_METS5/4-61                  KIKFKYKGQDLEVDISKVKKVWKVGKMVSFTYDD.NGKTGRGAVSEKDAPKELLNMIGK
DN7A_SACS2/3-61                 TVKFKYKGEEKQVDISKIKKVWRVGKMISFTYDEGGGKTGRGAVSEKDAPKELLQMLEK
#=GR DN7A_SACS2/3-61  SS        EEEEESSSSEEEEETTTEEEEEESSSSEEEEEE-SSSSEEEEEEETTTS-CHHHHHHTT
DN7E_SULAC/3-60                 KVRFKYKGEEKEVDTSKIKKVWRVGKMVSFTYDD.NGKTGRGAVSEKDAPKELMDMLAR
#=GC SS_cons                    EEEEESSSSEEEEETTTEEEEEESSSSEEEEEE-SSSSEEEEEEETTTS-CHHHHHHTT
#=GC seq_cons                   KVKFKYKGEEKEVDISKIKKVWRVGKMVSFTYDD.NGKTGRGAVSEKDAPKELLsMLuK
//

或者作为对齐的 FASTA（参见第对齐 FASTA 节）

>>> print(format(alignment, "fasta"))
>DN7_METS5/4-61
KIKFKYKGQDLEVDISKVKKVWKVGKMVSFTYDD-NGKTGRGAVSEKDAPKELLNMIGK
>DN7A_SACS2/3-61
TVKFKYKGEEKQVDISKIKKVWRVGKMISFTYDEGGGKTGRGAVSEKDAPKELLQMLEK
>DN7E_SULAC/3-60
KVRFKYKGEEKEVDTSKIKKVWRVGKMVSFTYDD-NGKTGRGAVSEKDAPKELMDMLAR

或在 PHYLIP 格式中（参见第 PHYLIP 输出文件节）

>>> print(format(alignment, "phylip"))
3 59
DN7_METS5/KIKFKYKGQDLEVDISKVKKVWKVGKMVSFTYDD-NGKTGRGAVSEKDAPKELLNMIGK
DN7A_SACS2TVKFKYKGEEKQVDISKIKKVWRVGKMISFTYDEGGGKTGRGAVSEKDAPKELLQMLEK
DN7E_SULACKVRFKYKGEEKEVDTSKIKKVWRVGKMVSFTYDD-NGKTGRGAVSEKDAPKELMDMLAR

对齐的一般信息存储在 Alignment 对象的 annotations 属性下，例如

>>> alignment.annotations["identifier"]
'7kD_DNA_binding'
>>> alignment.annotations["clan"]
'CL0049'
>>> alignment.annotations["database references"]
[{'reference': 'INTERPRO; IPR003212;'}, {'reference': 'SCOP; 1sso; fa;'}, {'reference': 'SO; 0000417; polypeptide_domain;'}]

此对齐中的各个序列存储在 alignment.sequences 中，作为 SeqRecord，包括与每个序列记录关联的任何注释

>>> for record in alignment.sequences:
...     print("%s %s %s" % (record.id, record.annotations["accession"], record.dbxrefs))
...
DN7_METS5/4-61 A4YEA2.1 []
DN7A_SACS2/3-61 P61991.2 ['PDB; 1SSO A; 2-60;', 'PDB; 1JIC A; 2-60;', 'PDB; 2CVR A; 2-60;', 'PDB; 1B4O A; 2-60;']
DN7E_SULAC/3-60 P13125.2 []

第二个序列 (DN7A_SACS2/3-61) 的二级结构存储在 SeqRecord 的 letter_annotations 属性中

>>> alignment.sequences[0].letter_annotations
{}
>>> alignment.sequences[1].letter_annotations
{'secondary structure': 'EEEEESSSSEEEEETTTEEEEEESSSSEEEEEE-SSSSEEEEEEETTTS-CHHHHHHTT'}
>>> alignment.sequences[2].letter_annotations
{}

共识序列和二级结构与整个序列对齐相关联，因此存储在 Alignment 对象的 column_annotations 属性中

>>> alignment.column_annotations  
{'consensus secondary structure': 'EEEEESSSSEEEEETTTEEEEEESSSSEEEEEE-SSSSEEEEEEETTTS-CHHHHHHTT',
 'consensus sequence': 'KVKFKYKGEEKEVDISKIKKVWRVGKMVSFTYDD.NGKTGRGAVSEKDAPKELLsMLuK'}

PHYLIP 输出文件

序列比对的 PHYLIP 格式源自 Joe Felsenstein 的 PHYLogeny Interference Package。 PHYLIP 格式的文件以两个数字开头，分别表示打印比对中的行数和列数。序列比对本身可以是顺序格式或交错格式。前者的一个例子是 sequential.phy 文件（在 Biopython 源代码分发包的 Tests/Phylip/ 中提供）。

 3 384
CYS1_DICDI   -----MKVIL LFVLAVFTVF VSS------- --------RG IPPEEQ---- --------SQ
             FLEFQDKFNK KY-SHEEYLE RFEIFKSNLG KIEELNLIAI NHKADTKFGV NKFADLSSDE
             FKNYYLNNKE AIFTDDLPVA DYLDDEFINS IPTAFDWRTR G-AVTPVKNQ GQCGSCWSFS
             TTGNVEGQHF ISQNKLVSLS EQNLVDCDHE CMEYEGEEAC DEGCNGGLQP NAYNYIIKNG
             GIQTESSYPY TAETGTQCNF NSANIGAKIS NFTMIP-KNE TVMAGYIVST GPLAIAADAV
             E-WQFYIGGV F-DIPCN--P NSLDHGILIV GYSAKNTIFR KNMPYWIVKN SWGADWGEQG
             YIYLRRGKNT CGVSNFVSTS II--
ALEU_HORVU   MAHARVLLLA LAVLATAAVA VASSSSFADS NPIRPVTDRA ASTLESAVLG ALGRTRHALR
             FARFAVRYGK SYESAAEVRR RFRIFSESLE EVRSTN---- RKGLPYRLGI NRFSDMSWEE
             FQATRL-GAA QTCSATLAGN HLMRDA--AA LPETKDWRED G-IVSPVKNQ AHCGSCWTFS
             TTGALEAAYT QATGKNISLS EQQLVDCAGG FNNF------ --GCNGGLPS QAFEYIKYNG
             GIDTEESYPY KGVNGV-CHY KAENAAVQVL DSVNITLNAE DELKNAVGLV RPVSVAFQVI
             DGFRQYKSGV YTSDHCGTTP DDVNHAVLAV GYGVENGV-- ---PYWLIKN SWGADWGDNG
             YFKMEMGKNM CAIATCASYP VVAA
CATH_HUMAN   ------MWAT LPLLCAGAWL LGV------- -PVCGAAELS VNSLEK---- --------FH
             FKSWMSKHRK TY-STEEYHH RLQTFASNWR KINAHN---- NGNHTFKMAL NQFSDMSFAE
             IKHKYLWSEP QNCSAT--KS NYLRGT--GP YPPSVDWRKK GNFVSPVKNQ GACGSCWTFS
             TTGALESAIA IATGKMLSLA EQQLVDCAQD FNNY------ --GCQGGLPS QAFEYILYNK
             GIMGEDTYPY QGKDGY-CKF QPGKAIGFVK DVANITIYDE EAMVEAVALY NPVSFAFEVT
             QDFMMYRTGI YSSTSCHKTP DKVNHAVLAV GYGEKNGI-- ---PYWIVKN SWGPQWGMNG
             YFLIERGKNM CGLAACASYP IPLV

在顺序格式中，一个序列的完整比对将在继续下一个序列之前显示。在交错格式中，不同序列的比对将彼此相邻，例如在文件 interlaced.phy 中（在 Biopython 源代码分发包的 Tests/Phylip/ 中提供）。

 3 384
CYS1_DICDI   -----MKVIL LFVLAVFTVF VSS------- --------RG IPPEEQ---- --------SQ
ALEU_HORVU   MAHARVLLLA LAVLATAAVA VASSSSFADS NPIRPVTDRA ASTLESAVLG ALGRTRHALR
CATH_HUMAN   ------MWAT LPLLCAGAWL LGV------- -PVCGAAELS VNSLEK---- --------FH

             FLEFQDKFNK KY-SHEEYLE RFEIFKSNLG KIEELNLIAI NHKADTKFGV NKFADLSSDE
             FARFAVRYGK SYESAAEVRR RFRIFSESLE EVRSTN---- RKGLPYRLGI NRFSDMSWEE
             FKSWMSKHRK TY-STEEYHH RLQTFASNWR KINAHN---- NGNHTFKMAL NQFSDMSFAE

             FKNYYLNNKE AIFTDDLPVA DYLDDEFINS IPTAFDWRTR G-AVTPVKNQ GQCGSCWSFS
             FQATRL-GAA QTCSATLAGN HLMRDA--AA LPETKDWRED G-IVSPVKNQ AHCGSCWTFS
             IKHKYLWSEP QNCSAT--KS NYLRGT--GP YPPSVDWRKK GNFVSPVKNQ GACGSCWTFS

             TTGNVEGQHF ISQNKLVSLS EQNLVDCDHE CMEYEGEEAC DEGCNGGLQP NAYNYIIKNG
             TTGALEAAYT QATGKNISLS EQQLVDCAGG FNNF------ --GCNGGLPS QAFEYIKYNG
             TTGALESAIA IATGKMLSLA EQQLVDCAQD FNNY------ --GCQGGLPS QAFEYILYNK

             GIQTESSYPY TAETGTQCNF NSANIGAKIS NFTMIP-KNE TVMAGYIVST GPLAIAADAV
             GIDTEESYPY KGVNGV-CHY KAENAAVQVL DSVNITLNAE DELKNAVGLV RPVSVAFQVI
             GIMGEDTYPY QGKDGY-CKF QPGKAIGFVK DVANITIYDE EAMVEAVALY NPVSFAFEVT

             E-WQFYIGGV F-DIPCN--P NSLDHGILIV GYSAKNTIFR KNMPYWIVKN SWGADWGEQG
             DGFRQYKSGV YTSDHCGTTP DDVNHAVLAV GYGVENGV-- ---PYWLIKN SWGADWGDNG
             QDFMMYRTGI YSSTSCHKTP DKVNHAVLAV GYGEKNGI-- ---PYWIVKN SWGPQWGMNG

             YIYLRRGKNT CGVSNFVSTS II--
             YFKMEMGKNM CAIATCASYP VVAA
             YFLIERGKNM CGLAACASYP IPLV

Bio.Align 中的解析器会根据文件内容检测它是顺序格式还是交错格式，然后进行适当解析。

>>> from Bio import Align
>>> alignment = Align.read("sequential.phy", "phylip")
>>> alignment  
<Alignment object (3 rows x 384 columns) at ...>
>>> alignment2 = Align.read("interlaced.phy", "phylip")
>>> alignment2  
<Alignment object (3 rows x 384 columns) at ...>
>>> alignment == alignment2
True

这里，两个比对被认为是相等的，如果它们具有相同的序列内容和相同的比对坐标。

>>> alignment.shape
(3, 384)
>>> print(alignment)
CYS1_DICD         0 -----MKVILLFVLAVFTVFVSS---------------RGIPPEEQ------------SQ
ALEU_HORV         0 MAHARVLLLALAVLATAAVAVASSSSFADSNPIRPVTDRAASTLESAVLGALGRTRHALR
CATH_HUMA         0 ------MWATLPLLCAGAWLLGV--------PVCGAAELSVNSLEK------------FH

CYS1_DICD        28 FLEFQDKFNKKY-SHEEYLERFEIFKSNLGKIEELNLIAINHKADTKFGVNKFADLSSDE
ALEU_HORV        60 FARFAVRYGKSYESAAEVRRRFRIFSESLEEVRSTN----RKGLPYRLGINRFSDMSWEE
CATH_HUMA        34 FKSWMSKHRKTY-STEEYHHRLQTFASNWRKINAHN----NGNHTFKMALNQFSDMSFAE

CYS1_DICD        87 FKNYYLNNKEAIFTDDLPVADYLDDEFINSIPTAFDWRTRG-AVTPVKNQGQCGSCWSFS
ALEU_HORV       116 FQATRL-GAAQTCSATLAGNHLMRDA--AALPETKDWREDG-IVSPVKNQAHCGSCWTFS
CATH_HUMA        89 IKHKYLWSEPQNCSAT--KSNYLRGT--GPYPPSVDWRKKGNFVSPVKNQGACGSCWTFS

CYS1_DICD       146 TTGNVEGQHFISQNKLVSLSEQNLVDCDHECMEYEGEEACDEGCNGGLQPNAYNYIIKNG
ALEU_HORV       172 TTGALEAAYTQATGKNISLSEQQLVDCAGGFNNF--------GCNGGLPSQAFEYIKYNG
CATH_HUMA       145 TTGALESAIAIATGKMLSLAEQQLVDCAQDFNNY--------GCQGGLPSQAFEYILYNK

CYS1_DICD       206 GIQTESSYPYTAETGTQCNFNSANIGAKISNFTMIP-KNETVMAGYIVSTGPLAIAADAV
ALEU_HORV       224 GIDTEESYPYKGVNGV-CHYKAENAAVQVLDSVNITLNAEDELKNAVGLVRPVSVAFQVI
CATH_HUMA       197 GIMGEDTYPYQGKDGY-CKFQPGKAIGFVKDVANITIYDEEAMVEAVALYNPVSFAFEVT

CYS1_DICD       265 E-WQFYIGGVF-DIPCN--PNSLDHGILIVGYSAKNTIFRKNMPYWIVKNSWGADWGEQG
ALEU_HORV       283 DGFRQYKSGVYTSDHCGTTPDDVNHAVLAVGYGVENGV-----PYWLIKNSWGADWGDNG
CATH_HUMA       256 QDFMMYRTGIYSSTSCHKTPDKVNHAVLAVGYGEKNGI-----PYWIVKNSWGPQWGMNG

CYS1_DICD       321 YIYLRRGKNTCGVSNFVSTSII-- 343
ALEU_HORV       338 YFKMEMGKNMCAIATCASYPVVAA 362
CATH_HUMA       311 YFLIERGKNMCGLAACASYPIPLV 335

在以 PHYLIP 格式输出比对时，Bio.Align 将每个比对序列写入一行。

>>> print(format(alignment, "phylip"))
3 384
CYS1_DICDI-----MKVILLFVLAVFTVFVSS---------------RGIPPEEQ------------SQFLEFQDKFNKKY-SHEEYLERFEIFKSNLGKIEELNLIAINHKADTKFGVNKFADLSSDEFKNYYLNNKEAIFTDDLPVADYLDDEFINSIPTAFDWRTRG-AVTPVKNQGQCGSCWSFSTTGNVEGQHFISQNKLVSLSEQNLVDCDHECMEYEGEEACDEGCNGGLQPNAYNYIIKNGGIQTESSYPYTAETGTQCNFNSANIGAKISNFTMIP-KNETVMAGYIVSTGPLAIAADAVE-WQFYIGGVF-DIPCN--PNSLDHGILIVGYSAKNTIFRKNMPYWIVKNSWGADWGEQGYIYLRRGKNTCGVSNFVSTSII--
ALEU_HORVUMAHARVLLLALAVLATAAVAVASSSSFADSNPIRPVTDRAASTLESAVLGALGRTRHALRFARFAVRYGKSYESAAEVRRRFRIFSESLEEVRSTN----RKGLPYRLGINRFSDMSWEEFQATRL-GAAQTCSATLAGNHLMRDA--AALPETKDWREDG-IVSPVKNQAHCGSCWTFSTTGALEAAYTQATGKNISLSEQQLVDCAGGFNNF--------GCNGGLPSQAFEYIKYNGGIDTEESYPYKGVNGV-CHYKAENAAVQVLDSVNITLNAEDELKNAVGLVRPVSVAFQVIDGFRQYKSGVYTSDHCGTTPDDVNHAVLAVGYGVENGV-----PYWLIKNSWGADWGDNGYFKMEMGKNMCAIATCASYPVVAA
CATH_HUMAN------MWATLPLLCAGAWLLGV--------PVCGAAELSVNSLEK------------FHFKSWMSKHRKTY-STEEYHHRLQTFASNWRKINAHN----NGNHTFKMALNQFSDMSFAEIKHKYLWSEPQNCSAT--KSNYLRGT--GPYPPSVDWRKKGNFVSPVKNQGACGSCWTFSTTGALESAIAIATGKMLSLAEQQLVDCAQDFNNY--------GCQGGLPSQAFEYILYNKGIMGEDTYPYQGKDGY-CKFQPGKAIGFVKDVANITIYDEEAMVEAVALYNPVSFAFEVTQDFMMYRTGIYSSTSCHKTPDKVNHAVLAVGYGEKNGI-----PYWIVKNSWGPQWGMNGYFLIERGKNMCGLAACASYPIPLV

我们可以以 PHYLIP 格式写入比对，解析结果，并确认它与原始比对对象相同。

>>> from io import StringIO
>>> stream = StringIO()
>>> Align.write(alignment, stream, "phylip")
1
>>> stream.seek(0)
0
>>> alignment3 = Align.read(stream, "phylip")
>>> alignment == alignment3
True
>>> [record.id for record in alignment.sequences]
['CYS1_DICDI', 'ALEU_HORVU', 'CATH_HUMAN']
>>> [record.id for record in alignment3.sequences]
['CYS1_DICDI', 'ALEU_HORVU', 'CATH_HUMAN']

EMBOSS

EMBOSS（欧洲分子生物学开放软件套件）是一组用于分子生物学和生物信息学的开源软件工具 [Rice2000]。它包括软件，如 needle 和 water，用于成对序列比对。这是一个由 water 程序生成的输出示例，用于 Smith-Waterman 局部成对序列比对（在 Biopython 分发包的 Tests/Emboss 目录中以 water.txt 的形式提供）。

########################################
# Program:  water
# Rundate:  Wed Jan 16 17:23:19 2002
# Report_file: stdout
########################################
#=======================================
#
# Aligned_sequences: 2
# 1: IXI_234
# 2: IXI_235
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 131
# Identity:     112/131 (85.5%)
# Similarity:   112/131 (85.5%)
# Gaps:          19/131 (14.5%)
# Score: 591.5
#
#
#=======================================

IXI_234            1 TSPASIRPPAGPSSRPAMVSSRRTRPSPPGPRRPTGRPCCSAAPRRPQAT     50
                     |||||||||||||||         ||||||||||||||||||||||||||
IXI_235            1 TSPASIRPPAGPSSR---------RPSPPGPRRPTGRPCCSAAPRRPQAT     41

IXI_234           51 GGWKTCSGTCTTSTSTRHRGRSGWSARTTTAACLRASRKSMRAACSRSAG    100
                     ||||||||||||||||||||||||          ||||||||||||||||
IXI_235           42 GGWKTCSGTCTTSTSTRHRGRSGW----------RASRKSMRAACSRSAG     81

IXI_234          101 SRPNRFAPTLMSSCITSTTGPPAWAGDRSHE    131
                     |||||||||||||||||||||||||||||||
IXI_235           82 SRPNRFAPTLMSSCITSTTGPPAWAGDRSHE    112


#---------------------------------------
#---------------------------------------

由于此输出文件只包含一个比对，我们可以使用 Align.read 直接提取它。这里，我们将使用 Align.parse，以便我们可以看到此 water 运行的元数据。

>>> from Bio import Align
>>> alignments = Align.parse("water.txt", "emboss")

alignments 的 metadata 属性存储文件中标题中显示的信息，包括用于生成输出的程序、程序运行的日期和时间、输出文件名以及使用的特定比对文件格式（默认情况下假定为 srspair）。

>>> alignments.metadata
{'Align_format': 'srspair', 'Program': 'water', 'Rundate': 'Wed Jan 16 17:23:19 2002', 'Report_file': 'stdout'}

要提取比对，我们使用

>>> alignment = next(alignments)
>>> alignment  
<Alignment object (2 rows x 131 columns) at ...>
>>> alignment.shape
(2, 131)
>>> print(alignment)
IXI_234           0 TSPASIRPPAGPSSRPAMVSSRRTRPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTC
                  0 |||||||||||||||---------||||||||||||||||||||||||||||||||||||
IXI_235           0 TSPASIRPPAGPSSR---------RPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTC

IXI_234          60 TTSTSTRHRGRSGWSARTTTAACLRASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTTG
                 60 ||||||||||||||----------||||||||||||||||||||||||||||||||||||
IXI_235          51 TTSTSTRHRGRSGW----------RASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTTG

IXI_234         120 PPAWAGDRSHE 131
                120 ||||||||||| 131
IXI_235         101 PPAWAGDRSHE 112

>>> print(alignment.coordinates)
[[  0  15  24  74  84 131]
 [  0  15  15  65  65 112]]

我们可以使用索引来提取比对的特定部分。

>>> alignment[0]
'TSPASIRPPAGPSSRPAMVSSRRTRPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTCTTSTSTRHRGRSGWSARTTTAACLRASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTTGPPAWAGDRSHE'
>>> alignment[1]
'TSPASIRPPAGPSSR---------RPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTCTTSTSTRHRGRSGW----------RASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTTGPPAWAGDRSHE'
>>> alignment[1, 10:30]
'GPSSR---------RPSPPG'

alignment 的 annotations 属性存储与该比对相关联的特定信息。

>>> alignment.annotations
{'Matrix': 'EBLOSUM62', 'Gap_penalty': 10.0, 'Extend_penalty': 0.5, 'Identity': 112, 'Similarity': 112, 'Gaps': 19, 'Score': 591.5}

间隙数、相同数和错配数也可以通过在 alignment 对象上调用 counts 方法来获得。

>>> alignment.counts()
AlignmentCounts(gaps=19, identities=112, mismatches=0)

其中 AlignmentCounts 是 Python 标准库中 collections 模块中的 namedtuple。

两个序列之间显示的共识行存储在 column_annotations 属性中。

>>> alignment.column_annotations
{'emboss_consensus': '|||||||||||||||         ||||||||||||||||||||||||||||||||||||||||||||||||||          |||||||||||||||||||||||||||||||||||||||||||||||'}

使用 format 函数（或 format 方法）以其他格式打印比对，例如 PHYLIP 格式（参见第 PHYLIP 输出文件节）。

>>> print(format(alignment, "phylip"))
2 131
IXI_234   TSPASIRPPAGPSSRPAMVSSRRTRPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTCTTSTSTRHRGRSGWSARTTTAACLRASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTTGPPAWAGDRSHE
IXI_235   TSPASIRPPAGPSSR---------RPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTCTTSTSTRHRGRSGW----------RASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTTGPPAWAGDRSHE

我们可以使用 alignment.sequences 获取单个序列。但是，由于这是一个成对比对，我们也可以使用 alignment.target 和 alignment.query 获取目标序列和查询序列。

>>> alignment.target
SeqRecord(seq=Seq('TSPASIRPPAGPSSRPAMVSSRRTRPSPPGPRRPTGRPCCSAAPRRPQATGGWK...SHE'), id='IXI_234', name='<unknown name>', description='<unknown description>', dbxrefs=[])
>>> alignment.query
SeqRecord(seq=Seq('TSPASIRPPAGPSSRRPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGTCTTS...SHE'), id='IXI_235', name='<unknown name>', description='<unknown description>', dbxrefs=[])

目前，Biopython 不支持以 EMBOSS 定义的输出格式写入序列比对。

GCG 多序列格式 (MSF)

多序列格式 (MSF) 是为了存储由 GCG（遗传学计算机组）程序集生成的多个序列比对而创建的。 Biopython 分发包的 Tests/msf 目录中的 W_prot.msf 文件是 MSF 格式的序列比对文件示例。该文件显示了 11 个蛋白质序列的比对。

!!AA_MULTIPLE_ALIGNMENT

   MSF: 99  Type: P  Oct 18, 2017  11:35  Check: 0 ..

 Name: W*01:01:01:01    Len:    99  Check: 7236  Weight:  1.00
 Name: W*01:01:01:02    Len:    99  Check: 7236  Weight:  1.00
 Name: W*01:01:01:03    Len:    99  Check: 7236  Weight:  1.00
 Name: W*01:01:01:04    Len:    99  Check: 7236  Weight:  1.00
 Name: W*01:01:01:05    Len:    99  Check: 7236  Weight:  1.00
 Name: W*01:01:01:06    Len:    99  Check: 7236  Weight:  1.00
 Name: W*02:01          Len:    93  Check: 9483  Weight:  1.00
 Name: W*03:01:01:01    Len:    93  Check: 9974  Weight:  1.00
 Name: W*03:01:01:02    Len:    93  Check: 9974  Weight:  1.00
 Name: W*04:01          Len:    93  Check: 9169  Weight:  1.00
 Name: W*05:01          Len:    99  Check: 7331  Weight:  1.00
//

  W*01:01:01:01  GLTPFNGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
  W*01:01:01:02  GLTPFNGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
  W*01:01:01:03  GLTPFNGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
  W*01:01:01:04  GLTPFNGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
  W*01:01:01:05  GLTPFNGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
  W*01:01:01:06  GLTPFNGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
        W*02:01  GLTPSNGYTA ATWTRTAASS VGMNIPYDGA SYLVRNQELR SWTAADKAAQ
  W*03:01:01:01  GLTPSSGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
  W*03:01:01:02  GLTPSSGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ
        W*04:01  GLTPSNGYTA ATWTRTAASS VGMNIPYDGA SYLVRNQELR SWTAADKAAQ
        W*05:01  GLTPSSGYTA ATWTRTAVSS VGMNIPYHGA SYLVRNQELR SWTAADKAAQ

  W*01:01:01:01  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL
  W*01:01:01:02  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL
  W*01:01:01:03  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL
  W*01:01:01:04  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL
  W*01:01:01:05  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL
  W*01:01:01:06  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL
        W*02:01  MPWRRNMQSC SKPTCREGGR SGSAKSLRMG RRRCTAQNPK RLT
  W*03:01:01:01  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK RLT
  W*03:01:01:02  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK RLT
        W*04:01  MPWRRNMQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK RLT
        W*05:01  MPWRRNRQSC SKPTCREGGR SGSAKSLRMG RRGCSAQNPK DSHDPPPHL

要使用 Biopython 解析此文件，请使用

>>> from Bio import Align
>>> alignment = Align.read("W_prot.msf", "msf")

解析器会跳过所有行，直到并包括以“MSF:”开头的行。以下行（直到“//”分隔符）由解析器读取以验证每个序列的长度。比对部分（在“//”分隔符之后）由解析器读取并存储为 Alignment 对象。

>>> alignment  
<Alignment object (11 rows x 99 columns) at ...>
>>> print(alignment)
W*01:01:0         0 GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*01:01:0         0 GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*01:01:0         0 GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*01:01:0         0 GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*01:01:0         0 GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*01:01:0         0 GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*02:01           0 GLTPSNGYTAATWTRTAASSVGMNIPYDGASYLVRNQELRSWTAADKAAQMPWRRNMQSC
W*03:01:0         0 GLTPSSGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*03:01:0         0 GLTPSSGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC
W*04:01           0 GLTPSNGYTAATWTRTAASSVGMNIPYDGASYLVRNQELRSWTAADKAAQMPWRRNMQSC
W*05:01           0 GLTPSSGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWRRNRQSC

W*01:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99
W*01:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99
W*01:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99
W*01:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99
W*01:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99
W*01:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99
W*02:01          60 SKPTCREGGRSGSAKSLRMGRRRCTAQNPKRLT------ 93
W*03:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKRLT------ 93
W*03:01:0        60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKRLT------ 93
W*04:01          60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKRLT------ 93
W*05:01          60 SKPTCREGGRSGSAKSLRMGRRGCSAQNPKDSHDPPPHL 99

序列及其名称存储在 alignment.sequences 属性中。

>>> len(alignment.sequences)
11
>>> alignment.sequences[0].id
'W*01:01:01:01'
>>> alignment.sequences[0].seq
Seq('GLTPFNGYTAATWTRTAVSSVGMNIPYHGASYLVRNQELRSWTAADKAAQMPWR...PHL')

比对坐标存储在 alignment.coordinates 属性中。

>>> print(alignment.coordinates)
[[ 0 93 99]
 [ 0 93 99]
 [ 0 93 99]
 [ 0 93 99]
 [ 0 93 99]
 [ 0 93 99]
 [ 0 93 93]
 [ 0 93 93]
 [ 0 93 93]
 [ 0 93 93]
 [ 0 93 99]]

目前，Biopython 不支持以 MSF 格式写入序列比对。

Exonerate

Exonerate 是一个用于成对序列比对的通用程序 [Slater2005]。 Exonerate 找到的序列比对可以以人类可读的形式输出，以“cigar”（紧凑的特殊间隙比对报告）格式输出，或以“vulgar”（详细的有用标记间隙比对报告）格式输出。用户可以请求在输出中包含一个或多个这些格式。 Bio.Align 中的解析器只能解析 cigar 或 vulgar 格式的比对，不会解析包含以人类可读格式输出的比对的输出。

Biopython 测试套件中的 exn_22_m_cdna2genome_vulgar.exn 文件是显示 vulgar 格式比对的 Exonerate 输出文件示例。

Command line: [exonerate -m cdna2genome ../scer_cad1.fa /media/Waterloo/Downloads/genomes/scer_s288c/scer_s288c.fa --bestn 3 --showalignment no --showcigar no --showvulgar yes]
Hostname: [blackbriar]
vulgar: gi|296143771|ref|NM_001180731.1| 0 1230 + gi|330443520|ref|NC_001136.10| 1319275 1318045 - 6146 M 1 1 C 3 3 M 1226 1226
vulgar: gi|296143771|ref|NM_001180731.1| 1230 0 - gi|330443520|ref|NC_001136.10| 1318045 1319275 + 6146 M 129 129 C 3 3 M 1098 1098
vulgar: gi|296143771|ref|NM_001180731.1| 0 516 + gi|330443688|ref|NC_001145.3| 85010 667216 + 518 M 11 11 G 1 0 M 15 15 G 2 0 M 4 4 G 1 0 M 1 1 G 1 0 M 8 8 G 4 0 M 17 17 5 0 2 I 0 168904 3 0 2 M 4 4 G 0 1 M 8 8 G 2 0 M 3 3 G 1 0 M 33 33 G 0 2 M 7 7 G 0 1 M 102 102 5 0 2 I 0 96820 3 0 2 M 14 14 G 0 2 M 10 10 G 2 0 M 5 5 G 0 2 M 10 10 G 2 0 M 4 4 G 0 1 M 20 20 G 1 0 M 15 15 G 0 1 M 5 5 G 3 0 M 4 4 5 0 2 I 0 122114 3 0 2 M 20 20 G 0 5 M 6 6 5 0 2 I 0 193835 3 0 2 M 12 12 G 0 2 M 5 5 G 1 0 M 7 7 G 0 2 M 1 1 G 0 1 M 12 12 C 75 75 M 6 6 G 1 0 M 4 4 G 0 1 M 2 2 G 0 1 M 3 3 G 0 1 M 41 41
-- completed exonerate analysis

该文件包含三个比对。要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("exn_22_m_cdna2genome_vulgar.exn", "exonerate")

字典 alignments.metadata 存储有关这些比对的通用信息，显示在输出文件顶部。

>>> alignments.metadata  
{'Program': 'exonerate',
 'Command line': 'exonerate -m cdna2genome ../scer_cad1.fa /media/Waterloo/Downloads/genomes/scer_s288c/scer_s288c.fa --bestn 3 --showalignment no --showcigar no --showvulgar yes',
 'Hostname': 'blackbriar'}

现在我们可以迭代比对。第一个比对，比对分数为 6146.0，没有间隙。

>>> alignment = next(alignments)
>>> alignment.score
6146.0
>>> print(alignment.coordinates)
[[1319275 1319274 1319271 1318045]
 [      0       1       4    1230]]
>>> print(alignment)  
gi|330443   1319275 ????????????????????????????????????????????????????????????
                  0 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
gi|296143         0 ????????????????????????????????????????????????????????????
...
gi|330443   1318075 ?????????????????????????????? 1318045
               1200 ||||||||||||||||||||||||||||||    1230
gi|296143      1200 ??????????????????????????????    1230

请注意，打印比对中的目标（第一个序列）位于反向链上，而查询（第二个序列）位于正向链上，目标坐标递减，查询坐标递增。使用 Python 的内置 format 函数以 exonerate 格式打印此比对将写入一个 vulgar 行。

>>> print(format(alignment, "exonerate"))
vulgar: gi|296143771|ref|NM_001180731.1| 0 1230 + gi|330443520|ref|NC_001136.10| 1319275 1318045 - 6146 M 1 1 C 3 3 M 1226 1226

使用 format 方法允许我们请求 vulgar 行（默认）或 cigar 行。

>>> print(alignment.format("exonerate", "vulgar"))
vulgar: gi|296143771|ref|NM_001180731.1| 0 1230 + gi|330443520|ref|NC_001136.10| 1319275 1318045 - 6146 M 1 1 C 3 3 M 1226 1226

>>> print(alignment.format("exonerate", "cigar"))
cigar: gi|296143771|ref|NM_001180731.1| 0 1230 + gi|330443520|ref|NC_001136.10| 1319275 1318045 - 6146 M 1 M 3 M 1226

vulgar 行包含有关比对的信息（在 M 1 1 C 3 3 M 1226 部分），这些信息在 cigar 行 M 1 M 3 M 1226 中缺失。 vulgar 行指定比对以单个比对的核苷酸开头，后面跟着形成密码子的三个比对的核苷酸 (C)，后面跟着 1226 个比对的核苷酸。在 cigar 行中，我们看到一个比对的核苷酸，后面跟着三个比对的核苷酸，后面跟着 1226 个比对的核苷酸；它没有指定三个比对的核苷酸形成密码子。这些来自 vulgar 行的信息存储在 operations 属性中。

>>> alignment.operations
bytearray(b'MCM')

有关其他操作码的定义，请参见 Exonerate 文档。

类似地，当调用 Bio.Align.write 以 vulgar 或 cigar 比对行为文件写入文件时，可以使用 "vulgar" 或 "cigar" 参数。

我们可以以 BED 和 PSL 格式打印比对。

>>> print(format(alignment, "bed"))  
gi|330443520|ref|NC_001136.10|  1318045 1319275 gi|296143771|ref|NM_001180731.1| 6146   -   1318045 1319275 0   3   1226,3,1,   0,1226,1229,

>>> print(format(alignment, "psl"))  
1230    0   0   0   0   0   0   0   -   gi|296143771|ref|NM_001180731.1|    1230    0   1230    gi|330443520|ref|NC_001136.10|  1319275 1318045 1319275 3   1226,3,1,   0,1226,1229,    1318045,1319271,1319274,

SAM 格式解析器定义了自己的（可选）operations 属性（第序列比对/映射 (SAM) 节），这与 Exonerate 格式解析器中定义的 operations 属性略有不同。由于 operations 属性是可选的，因此我们在以 SAM 格式打印比对之前将其删除。

>>> del alignment.operations
>>> print(format(alignment, "sam"))  
gi|296143771|ref|NM_001180731.1|    16  gi|330443520|ref|NC_001136.10|  1318046 255 1226M3M1M   *   0   0   *   *   AS:i:6146

第三个比对包含四个长间隙。

>>> alignment = next(alignments)  # second alignment
>>> alignment = next(alignments)  # third alignment
>>> print(alignment)  
gi|330443     85010 ???????????-???????????????--????-?-????????----????????????
                  0 |||||||||||-|||||||||||||||--||||-|-||||||||----||||||||||||
gi|296143         0 ????????????????????????????????????????????????????????????

gi|330443     85061 ????????????????????????????????????????????????????????????
                 60 |||||-------------------------------------------------------
gi|296143        60 ?????-------------------------------------------------------
...
gi|330443    666990 ????????????????????????????????????????????????????????????
             582000 --------------------------------------------------||||||||||
gi|296143       346 --------------------------------------------------??????????

gi|330443    667050 ?????????-??????????????????????????????????????????????????
             582060 ||--|||||-|||||||--|-|||||||||||||||||||||||||||||||||||||||
gi|296143       356 ??--?????????????--?-???????????????????????????????????????

gi|330443    667109 ??????????????????????????????????????????????????????-?????
             582120 ||||||||||||||||||||||||||||||||||||||||||||||||||||||-||||-
gi|296143       411 ???????????????????????????????????????????????????????????-

gi|330443    667168 ???????????????????????????????????????????????? 667216
             582180 ||-|||-||||||||||||||||||||||||||||||||||||||||| 582228
gi|296143       470 ??-???-?????????????????????????????????????????    516

>>> print(format(alignment, "exonerate"))  
vulgar: gi|296143771|ref|NM_001180731.1| 0 516 + gi|330443688|ref|NC_001145.3|
85010 667216 + 518 M 11 11 G 1 0 M 15 15 G 2 0 M 4 4 G 1 0 M 1 1 G 1 0 M 8 8
 G 4 0 M 17 17 5 0 2 I 0 168904 3 0 2 M 4 4 G 0 1 M 8 8 G 2 0 M 3 3 G 1 0
 M 33 33 G 0 2 M 7 7 G 0 1 M 102 102 5 0 2 I 0 96820 3 0 2 M 14 14 G 0 2 M 10 10
 G 2 0 M 5 5 G 0 2 M 10 10 G 2 0 M 4 4 G 0 1 M 20 20 G 1 0 M 15 15 G 0 1 M 5 5
 G 3 0 M 4 4 5 0 2 I 0 122114 3 0 2 M 20 20 G 0 5 M 6 6 5 0 2 I 0 193835 3 0 2
 M 12 12 G 0 2 M 5 5 G 1 0 M 7 7 G 0 2 M 1 1 G 0 1 M 12 12 C 75 75 M 6 6 G 1 0
 M 4 4 G 0 1 M 2 2 G 0 1 M 3 3 G 0 1 M 41 41

NEXUS

NEXUS 文件格式 [Maddison1997] 由几个程序用来存储系统发育信息。这是一个 NEXUS 格式的文件示例（在 Biopython 分发包的 Tests/Nexus 子目录中以 codonposset.nex 的形式提供）。

#NEXUS
[MacClade 4.05 registered to Computational Biologist, University]


BEGIN DATA;
       DIMENSIONS  NTAX=2 NCHAR=22;
       FORMAT DATATYPE=DNA  MISSING=? GAP=- ;
MATRIX
[                           10        20 ]
[                           .         .  ]

Aegotheles         AAAAAGGCATTGTGGTGGGAAT   [22]
Aerodramus         ?????????TTGTGGTGGGAAT   [13]
;
END;


BEGIN CODONS;
       CODONPOSSET * CodonPositions =
               N: 1-10,
               1: 11-22\3,
               2: 12-22\3,
               3: 13-22\3;
       CODESET  * UNTITLED = Universal: all ;
END;

通常，NEXUS 格式的文件可能要复杂得多。 Bio.Align 很大程度上依赖于 Bio.Nexus 中的 NEXUS 解析器（参见第使用 Bio.Phylo 进行系统发育分析章）从 NEXUS 文件中提取 Alignment 对象。

要读取此 NEXUS 文件中的比对，请使用

>>> from Bio import Align
>>> alignment = Align.read("codonposset.nex", "nexus")
>>> print(alignment)
Aegothele         0 AAAAAGGCATTGTGGTGGGAAT 22
                  0 .........||||||||||||| 22
Aerodramu         0 ?????????TTGTGGTGGGAAT 22

>>> alignment.shape
(2, 22)

序列存储在 sequences 属性下。

>>> alignment.sequences[0].id
'Aegotheles'
>>> alignment.sequences[0].seq
Seq('AAAAAGGCATTGTGGTGGGAAT')
>>> alignment.sequences[0].annotations
{'molecule_type': 'DNA'}
>>> alignment.sequences[1].id
'Aerodramus'
>>> alignment.sequences[1].seq
Seq('?????????TTGTGGTGGGAAT')
>>> alignment.sequences[1].annotations
{'molecule_type': 'DNA'}

要以 NEXUS 格式打印此比对，请使用

>>> print(format(alignment, "nexus"))
#NEXUS
begin data;
dimensions ntax=2 nchar=22;
format datatype=dna missing=? gap=-;
matrix
Aegotheles AAAAAGGCATTGTGGTGGGAAT
Aerodramus ?????????TTGTGGTGGGAAT
;
end;

类似地，您可以使用 Align.write(alignment, "myfilename.nex", "nexus") 以 NEXUS 格式将比对写入文件 myfilename.nex。

来自 BLAST 或 FASTA 的表格输出

表格输出中的比对输出由 FASTA 比对器 [Pearson1988] 使用 -m 8CB 或 -m 8CC 参数运行，或由 BLAST [Altschul1990] 使用 -outfmt 7 参数运行生成。

Biopython 源代码分发包的 Tests/Fasta 子目录中的 nucleotide_m8CC.txt 文件是使用 -m 8CC 参数由 FASTA 生成的输出文件示例。

# fasta36 -m 8CC seq/mgstm1.nt seq/gst.nlib
# FASTA 36.3.8h May, 2020
# Query: pGT875  - 657 nt
# Database: seq/gst.nlib
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score, aln_code
# 12 hits found
pGT875  pGT875  100.00  657 0   0   1   657 38  694 4.6e-191    655.6   657M
pGT875  RABGLTR 79.10   646 135 0   1   646 34  679 1.6e-116    408.0   646M
pGT875  BTGST   59.56   413 167 21  176 594 228 655 1.9e-07 45.7    149M1D7M1I17M3D60M5I6M1I13M2I13M4I30M2I6M2D112M
pGT875  RABGSTB 66.93   127 42  8   159 289 157 287 3.2e-07 45.0    15M2I17M2D11M1I58M1I11M1D7M1D8M
pGT875  OCDHPR  91.30   23  2   1   266 289 2303    2325    0.012   29.7    17M1D6M
...
# FASTA processed 1 queries

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("nucleotide_m8CC.txt", "tabular")

文件中头部的信息存储在alignments属性的metadata属性中。

>>> alignments.metadata  
{'Command line': 'fasta36 -m 8CC seq/mgstm1.nt seq/gst.nlib',
 'Program': 'FASTA',
 'Version': '36.3.8h May, 2020',
 'Database': 'seq/gst.nlib'}

通过遍历alignments提取特定的比对。例如，让我们转到第四个比对。

>>> alignment = next(alignments)
>>> alignment = next(alignments)
>>> alignment = next(alignments)
>>> alignment = next(alignments)
>>> print(alignment)
RABGSTB         156 ??????????????????????????????????--????????????????????????
                  0 |||||||||||||||--|||||||||||||||||--|||||||||||-||||||||||||
pGT875          158 ???????????????--??????????????????????????????-????????????

RABGSTB         214 ??????????????????????????????????????????????????????????-?
                 60 ||||||||||||||||||||||||||||||||||||||||||||||-|||||||||||-|
pGT875          215 ??????????????????????????????????????????????-?????????????

RABGSTB         273 ??????-???????? 287
                120 ||||||-|||||||| 135
pGT875          274 ??????????????? 289

>>> print(alignment.coordinates)
[[156 171 173 190 190 201 202 260 261 272 272 279 279 287]
 [158 173 173 190 192 203 203 261 261 272 273 280 281 289]]
>>> alignment.aligned
array([[[156, 171],
        [173, 190],
        [190, 201],
        [202, 260],
        [261, 272],
        [272, 279],
        [279, 287]],

       [[158, 173],
        [173, 190],
        [192, 203],
        [203, 261],
        [261, 272],
        [273, 280],
        [281, 289]]])

目标和查询序列的序列信息存储在target和query属性中（以及在alignment.sequences下）。

>>> alignment.target
SeqRecord(seq=Seq(None, length=287), id='RABGSTB', name='<unknown name>', description='<unknown description>', dbxrefs=[])
>>> alignment.query
SeqRecord(seq=Seq(None, length=657), id='pGT875', name='<unknown name>', description='<unknown description>', dbxrefs=[])

比对的信息存储在alignment属性的annotations属性下。

>>> alignment.annotations  
{'% identity': 66.93,
 'mismatches': 42,
 'gap opens': 8,
 'evalue': 3.2e-07,
 'bit score': 45.0}

特别是BLAST提供了许多选项，可以通过包含或排除特定列来自定义表格输出；有关详细信息，请参阅BLAST文档。此信息存储在字典alignment.annotations，alignment.target.annotations或alignment.query.annotations中，具体取决于情况。

HH-suite输出文件

HH-suite [Steinegger2019]中的hhsearch或hhblits生成了hhr格式的比对文件。例如，请参见Biopython测试套件中的文件2uvo_hhblits.hhr。

Query         2UVO:A|PDBID|CHAIN|SEQUENCE
Match_columns 171
No_of_seqs    1560 out of 4005
Neff          8.3
Searched_HMMs 34
Date          Fri Feb 15 16:34:13 2019
Command       hhblits -i 2uvoAh.fasta -d /pdb70

 No Hit                             Prob E-value P-value  Score    SS Cols Query HMM  Template HMM
  1 2uvo_A Agglutinin isolectin 1; 100.0 3.7E-34 4.8E-38  210.3   0.0  171    1-171     1-171 (171)
  2 2wga   ; lectin (agglutinin);   99.9 1.1E-30 1.4E-34  190.4   0.0  162    2-169     2-163 (164)
  3 1ulk_A Lectin-C; chitin-bindin  99.8 5.2E-24 6.6E-28  148.2   0.0  120    1-124     2-121 (126)
...
 31 4z8i_A BBTPGRP3, peptidoglycan  79.6    0.12 1.5E-05   36.1   0.0   37    1-37      9-54  (236)
 32 1wga   ; lectin (agglutinin);   40.4     2.6 0.00029   25.9   0.0  106   54-163    11-116 (164)

No 1
>2uvo_A Agglutinin isolectin 1; carbohydrate-binding protein, hevein domain, chitin-binding, GERM agglutinin, chitin-binding protein; HET: NDG NAG GOL; 1.40A {Triticum aestivum} PDB: 1wgc_A* 2cwg_A* 2x3t_A* 4aml_A* 7wga_A 9wga_A 2wgc_A 1wgt_A 1k7t_A* 1k7v_A* 1k7u_A 2x52_A* 1t0w_A*
Probab=99.95  E-value=3.7e-34  Score=210.31  Aligned_cols=171  Identities=100%  Similarity=2.050  Sum_probs=166.9

Q 2UVO:A|PDBID|C    1 ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRCGSQAGGATCTNNQCCSQYGYCGFGAEYCGAGCQG   80 (171)
Q Consensus         1 ~~cg~~~~~~~c~~~~CCs~~g~CG~~~~~c~~~c~~~~c~~~~~Cg~~~~~~~c~~~~CCs~~g~CG~~~~~c~~~c~~   80 (171)
                      ||||++.++..||++.|||+|+|||.+.+||+++||.+.|++..+|+++++.++|....|||.++||+.+.+||+.+||.
T Consensus         1 ~~cg~~~~~~~c~~~~CCS~~g~Cg~~~~~Cg~gC~~~~c~~~~~cg~~~~~~~c~~~~CCs~~g~Cg~~~~~c~~~c~~   80 (171)
T 2uvo_A            1 ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRCGSQAGGATCTNNQCCSQYGYCGFGAEYCGAGCQG   80 (171)
T ss_dssp             CBCBGGGTTBBCGGGCEECTTSBEEBSHHHHSTTCCBSSCSSCCBCBGGGTTBCCSTTCEECTTSBEEBSHHHHSTTCCB
T ss_pred             CCCCCCCCCcCCCCCCeeCCCCeECCCcccccCCccccccccccccCcccCCcccCCccccCCCceeCCCccccCCCccc
Confidence            79999999999999999999999999999999999999999999999999999999999999999999999999999999


Q 2UVO:A|PDBID|C   81 GPCRADIKCGSQAGGKLCPNNLCCSQWGFCGLGSEFCGGGCQSGACSTDKPCGKDAGGRVCTNNYCCSKWGSCGIGPGYC  160 (171)
Q Consensus        81 ~~~~~~~~Cg~~~~~~~c~~~~CCS~~G~CG~~~~~C~~~Cq~~~c~~~~~Cg~~~~~~~c~~~~CCS~~G~CG~~~~~C  160 (171)
                      +++++|+.|+...+++.||++.|||.|||||...+||+.+||+++|++|.+|++.+++++|..+.|||+++-||+...||
T Consensus        81 ~~~~~~~~cg~~~~~~~c~~~~CCs~~g~CG~~~~~C~~gCq~~~c~~~~~cg~~~~~~~c~~~~ccs~~g~Cg~~~~~C  160 (171)
T 2uvo_A           81 GPCRADIKCGSQAGGKLCPNNLCCSQWGFCGLGSEFCGGGCQSGACSTDKPCGKDAGGRVCTNNYCCSKWGSCGIGPGYC  160 (171)
T ss_dssp             SSCSSCCBCBGGGTTBCCGGGCEECTTSBEEBSHHHHSTTCCBSSCSSCCCCBTTTTTBCCSTTCEECTTSCEEBSHHHH
T ss_pred             ccccccccccccccCCCCCCCcccCCCCccCCCcccccCCCcCCccccccccccccccccCCCCCCcCCCCEecCchhhc
Confidence            99999999999988999999999999999999999999999999999999999999999999999999999999999999


Q 2UVO:A|PDBID|C  161 GAGCQSGGCDG  171 (171)
Q Consensus       161 ~~gCq~~~c~~  171 (171)
                      +++||++.|||
T Consensus       161 ~~~cq~~~~~~  171 (171)
T 2uvo_A          161 GAGCQSGGCDG  171 (171)
T ss_dssp             STTCCBSSCC-
T ss_pred             ccccccCCCCC
Confidence            99999999986


No 2
...


No 32
>1wga   ; lectin (agglutinin); NMR {}
Probab=40.43  E-value=2.6  Score=25.90  Aligned_cols=106  Identities=20%  Similarity=0.652  Sum_probs=54.7

Q 2UVO:A|PDBID|C   54 TCTNNQCCSQYGYCGFGAEYCGAGCQGGPCRADIKCGSQAGGKLCPNNLCCSQWGFCGLGSEFCGGGCQSGACSTDKPCG  133 (171)
Q Consensus        54 ~c~~~~CCs~~g~CG~~~~~c~~~c~~~~~~~~~~Cg~~~~~~~c~~~~CCS~~G~CG~~~~~C~~~Cq~~~c~~~~~Cg  133 (171)
                      .|....||.....|......|...|....|.....|...  ...|....||.....|......|...|....+.....|.
T Consensus        11 ~c~~~~cc~~~~~c~~~~~~c~~~c~~~~c~~~~~c~~~--~~~c~~~~cc~~~~~c~~~~~~c~~~c~~~~c~~~~~c~   88 (164)
T 1wga             11 XCXXXXCCXXXXXCXXXXXXCXXXCXXXXCXXXXXCXXX--XXXCXXXXCCXXXXXCXXXXXXCXXXCXXXXCXXXXXCX   88 (164)
T ss_pred             ccccccccccccccccccccccccccccccccccccccc--ccccccccccccccccccccccccccccccccccccccc
Confidence            344556666666666666566555543333223333321  234666677777777777766666655544332223333


Q 2UVO:A|PDBID|C  134 KDAGGRVCTNNYCCSKWGSCGIGPGYCGAG  163 (171)
Q Consensus       134 ~~~~~~~c~~~~CCS~~G~CG~~~~~C~~g  163 (171)
                      ..  ...|....||.....|......|...
T Consensus        89 ~~--~~~c~~~~cc~~~~~c~~~~~~c~~~  116 (164)
T 1wga             89 XX--XXXCXXXXCCXXXXXCXXXXXXCXXX  116 (164)
T ss_pred             cc--cccccccccccccccccccccccccc
Confidence            22  23344455555555555555544433


Done!

该文件包含三个部分。

包含有关比对的常规信息的头部；
包含每个获得的比对的单行摘要；
详细显示的连续比对。

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("2uvo_hhblits.hhr", "hhr")

大部分头部信息存储在alignments属性的metadata属性中。

>>> alignments.metadata  
{'Match_columns': 171,
 'No_of_seqs': (1560, 4005),
 'Neff': 8.3,
 'Searched_HMMs': 34,
 'Rundate': 'Fri Feb 15 16:34:13 2019',
 'Command line': 'hhblits -i 2uvoAh.fasta -d /pdb70'}

除了查询名称，它作为属性存储

>>> alignments.query_name
'2UVO:A|PDBID|CHAIN|SEQUENCE'

因为它将在每个比对中再次出现。

遍历比对。

>>> for alignment in alignments:
...     print(alignment.target.id)  
...
2uvo_A
2wga
1ulk_A
...
4z8i_A
1wga

让我们更详细地看一下第一个比对。

>>> alignments = iter(alignments)
>>> alignment = next(alignments)
>>> alignment  
<Alignment object (2 rows x 171 columns) at ...>
>>> print(alignment)
2uvo_A            0 ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRCGSQAGGATCTNNQC
                  0 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2UVO:A|PD         0 ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRCGSQAGGATCTNNQC

2uvo_A           60 CSQYGYCGFGAEYCGAGCQGGPCRADIKCGSQAGGKLCPNNLCCSQWGFCGLGSEFCGGG
                 60 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2UVO:A|PD        60 CSQYGYCGFGAEYCGAGCQGGPCRADIKCGSQAGGKLCPNNLCCSQWGFCGLGSEFCGGG

2uvo_A          120 CQSGACSTDKPCGKDAGGRVCTNNYCCSKWGSCGIGPGYCGAGCQSGGCDG 171
                120 ||||||||||||||||||||||||||||||||||||||||||||||||||| 171
2UVO:A|PD       120 CQSGACSTDKPCGKDAGGRVCTNNYCCSKWGSCGIGPGYCGAGCQSGGCDG 171

目标和查询序列存储在alignment.sequences中。由于这些是成对比对，我们也可以通过alignment.target和alignment.query访问它们。

>>> alignment.target is alignment.sequences[0]
True
>>> alignment.query is alignment.sequences[1]
True

查询的ID从alignments.query_name设置（请注意，hhr文件中比对中打印的查询ID是缩写的）。

>>> alignment.query.id
'2UVO:A|PDBID|CHAIN|SEQUENCE'

目标的ID取自序列比对块（以T 2uvo_A开头的行）。

>>> alignment.target.id
'2uvo_A'

目标和查询的序列内容根据此比对中可用的信息填充。

>>> alignment.target.seq
Seq('ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRCGSQAGGAT...CDG')
>>> alignment.query.seq
Seq('ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRCGSQAGGAT...CDG')

如果比对没有扩展到整个序列，则序列内容将不完整（部分定义的序列；请参见第序列内容部分定义的序列节）。

此比对块的第二行，以“>”开头，显示了隐藏马尔可夫模型的名称和描述，目标序列取自该模型。这些存储在alignment.target.annotations字典中的"hmm_name"和"hmm_description"键下。

>>> alignment.target.annotations  
{'hmm_name': '2uvo_A',
 'hmm_description': 'Agglutinin isolectin 1; carbohydrate-binding protein, hevein domain, chitin-binding, GERM agglutinin, chitin-binding protein; HET: NDG NAG GOL; 1.40A {Triticum aestivum} PDB: 1wgc_A* 2cwg_A* 2x3t_A* 4aml_A* 7wga_A 9wga_A 2wgc_A 1wgt_A 1k7t_A* 1k7v_A* 1k7u_A 2x52_A* 1t0w_A*'}

字典alignment.target.letter_annotations存储目标比对的共识序列，PSIPRED预测的二级结构以及DSSP确定的目标二级结构。

>>> alignment.target.letter_annotations  
{'Consensus': '~~cg~~~~~~~c~~~~CCS~~g~Cg~~~~~Cg~gC~~~~c~~~~~cg~~~~~~~c~~~~CCs~~g~Cg~~~~~c~~~c~~~~~~~~~~cg~~~~~~~c~~~~CCs~~g~CG~~~~~C~~gCq~~~c~~~~~cg~~~~~~~c~~~~ccs~~g~Cg~~~~~C~~~cq~~~~~~',
 'ss_pred': 'CCCCCCCCCcCCCCCCeeCCCCeECCCcccccCCccccccccccccCcccCCcccCCccccCCCceeCCCccccCCCcccccccccccccccccCCCCCCCcccCCCCccCCCcccccCCCcCCccccccccccccccccCCCCCCcCCCCEecCchhhcccccccCCCCC',
 'ss_dssp': 'CBCBGGGTTBBCGGGCEECTTSBEEBSHHHHSTTCCBSSCSSCCBCBGGGTTBCCSTTCEECTTSBEEBSHHHHSTTCCBSSCSSCCBCBGGGTTBCCGGGCEECTTSBEEBSHHHHSTTCCBSSCSSCCCCBTTTTTBCCSTTCEECTTSCEEBSHHHHSTTCCBSSCC '}

在本示例中，对于查询序列，只有共识序列可用。

>>> alignment.query.letter_annotations
{'Consensus': '~~cg~~~~~~~c~~~~CCs~~g~CG~~~~~c~~~c~~~~c~~~~~Cg~~~~~~~c~~~~CCs~~g~CG~~~~~c~~~c~~~~~~~~~~Cg~~~~~~~c~~~~CCS~~G~CG~~~~~C~~~Cq~~~c~~~~~Cg~~~~~~~c~~~~CCS~~G~CG~~~~~C~~gCq~~~c~~'}

字典alignment.annotations存储有关比对块第三行中显示的比对的信息。

>>> alignment.annotations  
{'Probab': 99.95,
 'E-value': 3.7e-34,
 'Score': 210.31,
 'Identities': 100.0,
 'Similarity': 2.05,
 'Sum_probs': 166.9}

成对比对的置信度值存储在alignment.column_annotations字典的"Confidence"键下。此字典还存储每个列的分数，该分数显示在每个比对块的查询和目标部分之间。

>>> alignment.column_annotations  
{'column score': '||||++.++..||++.|||+|+|||.+.+||+++||.+.|++..+|+++++.++|....|||.++||+.+.+||+.+||.+++++|+.|+...+++.||++.|||.|||||...+||+.+||+++|++|.+|++.+++++|..+.|||+++-||+...||+++||++.|||',
 'Confidence': '799999999999999999999999999999999999999999999999999999999999999999999999999999999999999999998899999999999999999999999999999999999999999999999999999999999999999999999999986'}

A2M

A2M文件是SAM序列比对和建模软件系统[Krogh1994]，[Hughey1996]中的align2model或hmmscore创建的比对文件。A2M文件包含一个多重比对。A2M文件格式类似于对齐的FASTA（请参见对齐的FASTA部分）。但是，为了区分插入和删除，A2M使用破折号和句点来表示间隙，并在对齐序列中使用大小写字母。匹配用大写字母表示，删除用破折号表示，在仅包含匹配或删除的对齐列中。插入用小写字母表示，将间隙对齐到显示为句点的插入。

Biopython测试套件中的文件probcons.a2m是A2M文件的示例（请参见对齐的FASTA部分，了解以对齐的FASTA格式显示的相同比对）。

>plas_horvu
D.VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG.VD.VSKISQEEYLTAPGETFSVTLTV...PGTYGFYCEPHAGAGMVGKVT
V
>plas_chlre
-.VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG.VN.ADAISRDDYLNAPGETYSVKLTA...AGEYGYYCEPHQGAGMVGKII
V
>plas_anava
-.VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKsADlAKSLSHKQLLMSPGQSTSTTFPAdapAGEYTFYCEPHRGAGMVGKIT
V
>plas_proho
VqIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG.ES.APALSNTKLRIAPGSFYSVTLGT...PGTYSFYCTPHRGAGMVGTIT
V
>azup_achcy
VhMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG.AE.A-------FKSKINENYKVTFTA...PGVYGVKCTPHYGMGMVGVVE
V

要解析此比对，请使用

>>> from Bio import Align
>>> alignment = Align.read("probcons.a2m", "a2m")
>>> alignment  
<Alignment object (5 rows x 101 columns) at ...>
>>> print(alignment)
plas_horv         0 D-VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG-VD-VSKISQE
plas_chlr         0 --VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG-VN-ADAISRD
plas_anav         0 --VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKSADLAKSLSHK
plas_proh         0 VQIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG-ES-APALSNT
azup_achc         0 VHMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG-AE-A------

plas_horv        57 EYLTAPGETFSVTLTV---PGTYGFYCEPHAGAGMVGKVTV 95
plas_chlr        56 DYLNAPGETYSVKLTA---AGEYGYYCEPHQGAGMVGKIIV 94
plas_anav        58 QLLMSPGQSTSTTFPADAPAGEYTFYCEPHRGAGMVGKITV 99
plas_proh        56 KLRIAPGSFYSVTLGT---PGTYSFYCTPHRGAGMVGTITV 94
azup_achc        51 -FKSKINENYKVTFTA---PGVYGVKCTPHYGMGMVGVVEV 88

解析器分析A2M文件中破折号、句点以及大小写字母的模式，以确定一列是匹配/不匹配/删除（“D”）还是插入（“I”）。此信息存储在alignment.column_annotations字典的match键下。

>>> alignment.column_annotations
{'state': 'DIDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDIDDIDDDDDDDDDDDDDDDDDDDDDDDIIIDDDDDDDDDDDDDDDDDDDDDD'}

由于状态信息存储在alignment中，因此我们可以以A2M格式打印比对。

>>> print(format(alignment, "a2m"))
>plas_horvu
D.VLLGANGGVLVFEPNDFSVKAGETITFKNNAGYPHNVVFDEDAVPSG.VD.VSKISQEEYLTAPGETFSVTLTV...PGTYGFYCEPHAGAGMVGKVTV
>plas_chlre
-.VKLGADSGALEFVPKTLTIKSGETVNFVNNAGFPHNIVFDEDAIPSG.VN.ADAISRDDYLNAPGETYSVKLTA...AGEYGYYCEPHQGAGMVGKIIV
>plas_anava
-.VKLGSDKGLLVFEPAKLTIKPGDTVEFLNNKVPPHNVVFDAALNPAKsADlAKSLSHKQLLMSPGQSTSTTFPAdapAGEYTFYCEPHRGAGMVGKITV
>plas_proho
VqIKMGTDKYAPLYEPKALSISAGDTVEFVMNKVGPHNVIFDK--VPAG.ES.APALSNTKLRIAPGSFYSVTLGT...PGTYSFYCTPHRGAGMVGTITV
>azup_achcy
VhMLNKGKDGAMVFEPASLKVAPGDTVTFIPTDK-GHNVETIKGMIPDG.AE.A-------FKSKINENYKVTFTA...PGVYGVKCTPHYGMGMVGVVEV

类似地，可以使用Align.write将比对以A2M格式写入输出文件（请参见写入比对部分）。

Mauve扩展的多FASTA（xmfa）格式

Mauve [Darling2004]是一个用于构建多基因组比对的软件包。这些比对存储在扩展的多FASTA（xmfa）格式中。根据progressiveMauve（Mauve中的比对程序）的调用方式，xmfa格式略有不同。

如果使用单个序列输入文件调用progressiveMauve，如

progressiveMauve combined.fasta  --output=combined.xmfa ...

其中combined.fasta包含基因组序列

>equCab1
GAAAAGGAAAGTACGGCCCGGCCACTCCGGGTGTGTGCTAGGAGGGCTTA
>mm9
GAAGAGGAAAAGTAGATCCCTGGCGTCCGGAGCTGGGACGT
>canFam2
CAAGCCCTGCGCGCTCAGCCGGAGTGTCCCGGGCCCTGCTTTCCTTTTC

那么输出文件combined.xmfa如下所示

#FormatVersion Mauve1
#Sequence1File  combined.fa
#Sequence1Entry 1
#Sequence1Format    FastA
#Sequence2File  combined.fa
#Sequence2Entry 2
#Sequence2Format    FastA
#Sequence3File  combined.fa
#Sequence3Entry 3
#Sequence3Format    FastA
#BackboneFile   combined.xmfa.bbcols
> 1:2-49 - combined.fa
AAGCCCTCCTAGCACACACCCGGAGTGG-CCGGGCCGTACTTTCCTTTT
> 2:0-0 + combined.fa
-------------------------------------------------
> 3:2-48 + combined.fa
AAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGCTTTCCTTTT
=
> 1:1-1 + combined.fa
G
=
> 1:50-50 + combined.fa
A
=
> 2:1-41 + combined.fa
GAAGAGGAAAAGTAGATCCCTGGCGTCCGGAGCTGGGACGT
=
> 3:1-1 + combined.fa
C
=
> 3:49-49 + combined.fa
C
=

数字（1、2、3）分别指代马匹（equCab1）、小鼠（mm9）和狗（canFam2）的输入基因组序列。此xmfa文件包含六个比对块，由=字符隔开。使用Align.parse提取这些比对。

>>> from Bio import Align
>>> alignments = Align.parse("combined.xmfa", "mauve")

文件头部数据存储在metadata属性中。

>>> alignments.metadata  
{'FormatVersion': 'Mauve1',
 'BackboneFile': 'combined.xmfa.bbcols',
 'File': 'combined.fa'}

属性identifiers存储三个序列的序列标识符，在本例中是三个数字。

>>> alignments.identifiers
['0', '1', '2']

这些标识符用于各个比对中。

>>> for alignment in alignments:
...     print([record.id for record in alignment.sequences])
...     print(alignment)
...     print("******")
...
['0', '1', '2']
0                49 AAGCCCTCCTAGCACACACCCGGAGTGG-CCGGGCCGTACTTTCCTTTT  1
1                 0 -------------------------------------------------  0
2                 1 AAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGCTTTCCTTTT 48

******
['0']
0                 0 G 1

******
['0']
0                49 A 50

******
['1']
1                 0 GAAGAGGAAAAGTAGATCCCTGGCGTCCGGAGCTGGGACGT 41

******
['2']
2                 0 C 1

******
['2']
2                48 C 49

******

请注意，只有第一个块是真正的比对；其他块只包含一个序列。通过包含这些块，xmfa文件包含了在combined.fa输入文件中提供的完整序列。

如果使用单独的输入文件为每个基因组调用progressiveMauve，如

progressiveMauve equCab1.fa canFam2.fa mm9.fa --output=separate.xmfa ...

其中每个Fasta文件仅包含一个物种的基因组序列，那么输出文件separate.xmfa如下所示

#FormatVersion Mauve1
#Sequence1File  equCab1.fa
#Sequence1Format    FastA
#Sequence2File  canFam2.fa
#Sequence2Format    FastA
#Sequence3File  mm9.fa
#Sequence3Format    FastA
#BackboneFile   separate.xmfa.bbcols
> 1:1-50 - equCab1.fa
TAAGCCCTCCTAGCACACACCCGGAGTGGCC-GGGCCGTAC-TTTCCTTTTC
> 2:1-49 + canFam2.fa
CAAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGC-TTTCCTTTTC
> 3:1-19 - mm9.fa
---------------------------------GGATCTACTTTTCCTCTTC
=
> 3:20-41 + mm9.fa
CTGGCGTCCGGAGCTGGGACGT
=

马匹的标识符equCab1，小鼠的标识符mm9以及狗的标识符canFam2现在在输出文件中显式显示。此xmfa文件包含两个比对块，由=字符隔开。使用Align.parse提取这些比对。

>>> from Bio import Align
>>> alignments = Align.parse("separate.xmfa", "mauve")

文件头部数据现在不包含输入文件名。

>>> alignments.metadata  
{'FormatVersion': 'Mauve1',
 'BackboneFile': 'separate.xmfa.bbcols'}

属性identifiers存储三个序列的序列标识符。

>>> alignments.identifiers
['equCab1.fa', 'canFam2.fa', 'mm9.fa']

这些标识符用于各个比对中。

>>> for alignment in alignments:
...     print([record.id for record in alignment.sequences])
...     print(alignment)
...     print("******")
...
['equCab1.fa', 'canFam2.fa', 'mm9.fa']
equCab1.f        50 TAAGCCCTCCTAGCACACACCCGGAGTGGCC-GGGCCGTAC-TTTCCTTTTC  0
canFam2.f         0 CAAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGC-TTTCCTTTTC 49
mm9.fa           19 ---------------------------------GGATCTACTTTTCCTCTTC  0

******
['mm9.fa']
mm9.fa           19 CTGGCGTCCGGAGCTGGGACGT 41

******

要以Mauve格式输出比对，请使用Align.write。

>>> from io import StringIO
>>> stream = StringIO()
>>> alignments = Align.parse("separate.xmfa", "mauve")
>>> Align.write(alignments, stream, "mauve")
2
>>> print(stream.getvalue())  
#FormatVersion Mauve1
#Sequence1File  equCab1.fa
#Sequence1Format    FastA
#Sequence2File  canFam2.fa
#Sequence2Format    FastA
#Sequence3File  mm9.fa
#Sequence3Format    FastA
#BackboneFile   separate.xmfa.bbcols
> 1:1-50 - equCab1.fa
TAAGCCCTCCTAGCACACACCCGGAGTGGCC-GGGCCGTAC-TTTCCTTTTC
> 2:1-49 + canFam2.fa
CAAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGC-TTTCCTTTTC
> 3:1-19 - mm9.fa
---------------------------------GGATCTACTTTTCCTCTTC
=
> 3:20-41 + mm9.fa
CTGGCGTCCGGAGCTGGGACGT
=

在这里，编写器利用存储在alignments.metadata和alignments.identifiers中的信息来创建此格式。如果您的alignments对象没有这些属性，则可以将它们作为关键字参数提供给Align.write。

>>> stream = StringIO()
>>> alignments = Align.parse("separate.xmfa", "mauve")
>>> metadata = alignments.metadata
>>> identifiers = alignments.identifiers
>>> alignments = list(alignments)  # this drops the attributes
>>> alignments.metadata  
Traceback (most recent call last):
 ...
AttributeError: 'list' object has no attribute 'metadata'
>>> alignments.identifiers  
Traceback (most recent call last):
 ...
AttributeError: 'list' object has no attribute 'identifiers'
>>> Align.write(alignments, stream, "mauve", metadata=metadata, identifiers=identifiers)
2
>>> print(stream.getvalue())  
#FormatVersion Mauve1
#Sequence1File  equCab1.fa
#Sequence1Format    FastA
#Sequence2File  canFam2.fa
#Sequence2Format    FastA
#Sequence3File  mm9.fa
#Sequence3Format    FastA
#BackboneFile   separate.xmfa.bbcols
> 1:1-50 - equCab1.fa
TAAGCCCTCCTAGCACACACCCGGAGTGGCC-GGGCCGTAC-TTTCCTTTTC
> 2:1-49 + canFam2.fa
CAAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGC-TTTCCTTTTC
> 3:1-19 - mm9.fa
---------------------------------GGATCTACTTTTCCTCTTC
=
> 3:20-41 + mm9.fa
CTGGCGTCCGGAGCTGGGACGT
=

Python不允许您直接将这些属性添加到alignments对象中，因为在本示例中它被转换为一个普通列表。但是，您可以构造一个Alignments对象并向其添加属性（请参见第Alignments类节）。

>>> alignments = Align.Alignments(alignments)
>>> alignments.metadata = metadata
>>> alignments.identifiers = identifiers
>>> stream = StringIO()
>>> Align.write(alignments, stream, "mauve", metadata=metadata, identifiers=identifiers)
2
>>> print(stream.getvalue())  
#FormatVersion Mauve1
#Sequence1File  equCab1.fa
#Sequence1Format    FastA
#Sequence2File  canFam2.fa
#Sequence2Format    FastA
#Sequence3File  mm9.fa
#Sequence3Format    FastA
#BackboneFile   separate.xmfa.bbcols
> 1:1-50 - equCab1.fa
TAAGCCCTCCTAGCACACACCCGGAGTGGCC-GGGCCGTAC-TTTCCTTTTC
> 2:1-49 + canFam2.fa
CAAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGC-TTTCCTTTTC
> 3:1-19 - mm9.fa
---------------------------------GGATCTACTTTTCCTCTTC
=
> 3:20-41 + mm9.fa
CTGGCGTCCGGAGCTGGGACGT
=

当以Mauve格式打印单个比对时，请使用关键字参数提供元数据和标识符。

>>> alignment = alignments[0]
>>> print(alignment.format("mauve", metadata=metadata, identifiers=identifiers))
> 1:1-50 - equCab1.fa
TAAGCCCTCCTAGCACACACCCGGAGTGGCC-GGGCCGTAC-TTTCCTTTTC
> 2:1-49 + canFam2.fa
CAAGCCCTGC--GCGCTCAGCCGGAGTGTCCCGGGCCCTGC-TTTCCTTTTC
> 3:1-19 - mm9.fa
---------------------------------GGATCTACTTTTCCTCTTC
=

序列比对/映射（SAM）

序列比对/映射（SAM）格式[Li2009]的文件存储成对序列比对，通常是下一代测序数据与参考基因组的比对。Biopython测试套件中的文件ex1.sam是SAM格式的最小文件的示例。它的前几行如下所示

EAS56_57:6:190:289:82   69      chr1    100     0       *       =       100     0       CTCAAGGTTGTTGCAAGGGGGTCTATGTGAACAAA     <<<7<<<;<<<<<<<<8;;<7;4<;<;;;;;94<;     MF:i:192
EAS56_57:6:190:289:82   137     chr1    100     73      35M     =       100     0       AGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCAC     <<<<<<;<<<<<<<<<<;<<;<<<<;8<6;9;;2;     MF:i:64 Aq:i:0  NM:i:0  UQ:i:0  H0:i:1  H1:i:0
EAS51_64:3:190:727:308  99      chr1    103     99      35M     =       263     195     GGTGCAGAGCCGAGTCACGGGGTTGCCAGCACAGG     <<<<<<<<<<<<<<<<<<<<<<<<<<<::<<<844     MF:i:18 Aq:i:73 NM:i:0  UQ:i:0  H0:i:1  H1:i:0
...

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("ex1.sam", "sam")
>>> alignment = next(alignments)

第一行的flag是69。根据SAM/BAM文件格式规范，对于标志包含位标志4的行，它们未映射。由于69在对应于此位置的位被设置为True，因此该序列未映射，并且未与基因组比对（尽管第一行显示了chr1）。因此，此比对的目标（或alignment.sequences中的第一个项目）是None。

>>> alignment.flag
69
>>> bin(69)
'0b1000101'
>>> bin(4)
'0b100'
>>> if alignment.flag & 4:
...     print("unmapped")
... else:
...     print("mapped")
...
unmapped
>>> alignment.sequences
[None, SeqRecord(seq=Seq('CTCAAGGTTGTTGCAAGGGGGTCTATGTGAACAAA'), id='EAS56_57:6:190:289:82', name='<unknown name>', description='', dbxrefs=[])]
>>> alignment.target is None
True

第二行表示与染色体1的比对。

>>> alignment = next(alignments)
>>> if alignment.flag & 4:
...     print("unmapped")
... else:
...     print("mapped")
...
mapped
>>> alignment.target
SeqRecord(seq=None, id='chr1', name='<unknown name>', description='', dbxrefs=[])

由于此 SAM 文件未存储每个比对的基因组序列信息，因此我们无法打印比对结果。但是，我们可以以 SAM 格式或任何其他格式（例如 BED，请参见第浏览器扩展数据 (BED)节）打印比对信息，这些格式不需要目标序列信息。

>>> format(alignment, "sam")
'EAS56_57:6:190:289:82\t137\tchr1\t100\t73\t35M\t=\t100\t0\tAGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCAC\t<<<<<<;<<<<<<<<<<;<<;<<<<;8<6;9;;2;\tMF:i:64\tAq:i:0\tNM:i:0\tUQ:i:0\tH0:i:1\tH1:i:0\n'
>>> format(alignment, "bed")
'chr1\t99\t134\tEAS56_57:6:190:289:82\t0\t+\t99\t134\t0\t1\t35,\t0,\n'

但是，我们无法以 PSL 格式（请参见第模式空间布局 (PSL)节）打印比对结果，因为这将需要知道目标序列 chr1 的大小。

>>> format(alignment, "psl")  
Traceback (most recent call last):
 ...
TypeError: ...

如果您知道目标序列的大小，您可以手动设置它们。

>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> target = SeqRecord(Seq(None, length=1575), id="chr1")
>>> alignment.target = target
>>> format(alignment, "psl")  
'35\t0\t0\t0\t0\t0\t0\t0\t+\tEAS56_57:6:190:289:82\t35\t0\t35\tchr1\t1575\t99\t134\t1\t35,\t0,\t99,\n'

Biopython 测试套件中的 ex1_header.sam 文件包含相同的比对，但现在还包括一个标题。其前几行如下所示：

@HD\tVN:1.3\tSO:coordinate
@SQ\tSN:chr1\tLN:1575
@SQ\tSN:chr2\tLN:1584
EAS56_57:6:190:289:82   69      chr1    100     0       *       =       100     0       CTCAAGGTTGTTGCAAGGGGGTCTATGTGAACAAA     <<<7<<<;<<<<<<<<8;;<7;4<;<;;;;;94<;     MF:i:192
...

标题存储有关比对的一般信息，包括目标染色体的大小。目标信息存储在 alignments 对象的 targets 属性中。

>>> from Bio import Align
>>> alignments = Align.parse("ex1_header.sam", "sam")
>>> len(alignments.targets)
2
>>> alignments.targets[0]
SeqRecord(seq=Seq(None, length=1575), id='chr1', name='<unknown name>', description='', dbxrefs=[])
>>> alignments.targets[1]
SeqRecord(seq=Seq(None, length=1584), id='chr2', name='<unknown name>', description='', dbxrefs=[])

标题中提供的其他信息存储在 metadata 属性中。

>>> alignments.metadata
{'HD': {'VN': '1.3', 'SO': 'coordinate'}}

有了目标信息，我们现在也可以以 PSL 格式打印比对结果。

>>> alignment = next(alignments)  # the unmapped sequence; skip it
>>> alignment = next(alignments)
>>> format(alignment, "psl")
'35\t0\t0\t0\t0\t0\t0\t0\t+\tEAS56_57:6:190:289:82\t35\t0\t35\tchr1\t1575\t99\t134\t1\t35,\t0,\t99,\n'

我们现在也可以以人类可读的形式打印比对结果，但请注意，目标序列内容无法从该文件获得。

>>> print(alignment)
chr1             99 ??????????????????????????????????? 134
                  0 ...................................  35
EAS56_57:         0 AGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCAC  35

Biopython 测试套件中 sam1.sam 文件中的比对包含一个额外的 MD 标签，该标签显示查询序列与目标序列的不同之处。

@SQ     SN:1    LN:239940
@PG     ID:bwa  PN:bwa  VN:0.6.2-r126
HWI-1KL120:88:D0LRBACXX:1:1101:1780:2146        77      *       0       0       *       *       0       0       GATGGGAAACCCATGGCCGAGTGGGAAGAAACCAGCTGAGGTCACATCACCAGAGGAGGGAGAGTGTGGCCCCTGACTCAGTCCATCAGCTTGTGGAGCTG   @=?DDDDBFFFF7A;E?GGEGE8BB?FF?F>G@F=GIIDEIBCFF<FEFEC@EEEE2?8B8/=@((-;?@2<B9@##########################
...
HWI-1KL120:88:D0LRBACXX:1:1101:2852:2134        137     1       136186  25      101M    =       136186  0       TCACGGTGGCCTGTTGAGGCAGGGGCTCACGCTGACCTCTCTCGGCGTGGGAGGGGCCGGTGTGAGGCAAGGGCTCACGCTGACCTCTCTCGGCGTGGGAG   @C@FFFDFHGHHHJJJIJJJJIJJJGEDHHGGHGBGIIGIIAB@GEE=BDBBCCDD@D@B7@;@DDD?<A?DD728:>8()009>:>>C@>5??B######   XT:A:U  NM:i:5  SM:i:25 AM:i:0  X0:i:1  X1:i:0  XM:i:5  XO:i:0  XG:i:0  MD:Z:25G14G2C34A12A9

解析器从 MD 标签重建局部基因组序列，使我们能够在打印比对结果时显式查看目标序列。

>>> from Bio import Align
>>> alignments = Align.parse("sam1.sam", "sam")
>>> for alignment in alignments:
...     if not alignment.flag & 4:  # Skip the unmapped lines
...         break
...
>>> alignment  
<Alignment object (2 rows x 101 columns) at ...>
>>> print(alignment)
1            136185 TCACGGTGGCCTGTTGAGGCAGGGGGTCACGCTGACCTCTGTCCGCGTGGGAGGGGCCGG
                  0 |||||||||||||||||||||||||.||||||||||||||.||.||||||||||||||||
HWI-1KL12         0 TCACGGTGGCCTGTTGAGGCAGGGGCTCACGCTGACCTCTCTCGGCGTGGGAGGGGCCGG

1            136245 TGTGAGGCAAGGGCTCACACTGACCTCTCTCAGCGTGGGAG 136286
                 60 ||||||||||||||||||.||||||||||||.|||||||||    101
HWI-1KL12        60 TGTGAGGCAAGGGCTCACGCTGACCTCTCTCGGCGTGGGAG    101

SAM 文件可能包含其他信息，以区分简单的序列插入和删除与基因组的跳过区域（例如内含子）、硬剪切和软剪切以及填充序列区域。由于这些信息无法存储在 Alignment 对象的 coordinates 属性中，而是存储在专用的 operations 属性中。让我们使用此 SAM 文件中的第三个比对作为示例。

>>> from Bio import Align
>>> alignments = Align.parse("dna_rna.sam", "sam")
>>> alignment = next(alignments)
>>> alignment = next(alignments)
>>> alignment = next(alignments)
>>> print(format(alignment, "SAM"))  
NR_111921.1 0   chr3    48663768    0   46M1827N82M3376N76M12H  *   0   0   CACGAGAGGAGCGGAGGCGAGGGGTGAACGCGGAGCACTCCAATCGCTCCCAACTAGAGGTCCACCCAGGACCCAGAGACCTGGATTTGAGGCTGCTGGGCGGCAGATGGAGCGATCAGAAGACCAGGAGACGGGAGCTGGAGTGCAGTGGCTGTTCACAAGCGTGAAAGCAAAGATTAAAAAATTTGTTTTTATATTAAAAAA    *   AS:i:1000   NM:i:0

>>> print(alignment.coordinates)
[[48663767 48663813 48665640 48665722 48669098 48669174]
 [       0       46       46      128      128      204]]
>>> alignment.operations
bytearray(b'MNMNM')
>>> alignment.query.annotations["hard_clip_right"]
12

在此比对中，cigar 字符串 63M1062N75M468N43M 定义了 46 个比对的核苷酸、1827 个核苷酸的内含子、82 个比对的核苷酸、3376 个核苷酸的内含子、76 个比对的核苷酸和 12 个硬剪切的核苷酸。这些操作显示在 operations 属性中，除了硬剪切，硬剪切存储在 alignment.query.annotations["hard_clip_right"]（或 alignment.query.annotations["hard_clip_left"]，如果适用）中。

要写入包含从头创建的比对的 SAM 文件，请使用 Alignments（复数）对象（请参见第Alignments 类节）来存储比对以及元数据和目标。

>>> from io import StringIO
>>> import numpy as np

>>> from Bio import Align
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord

>>> alignments = Align.Alignments()

>>> seq1 = Seq(None, length=10000)
>>> target1 = SeqRecord(seq1, id="chr1")
>>> seq2 = Seq(None, length=15000)
>>> target2 = SeqRecord(seq2, id="chr2")
>>> alignments.targets = [target1, target2]
>>> alignments.metadata = {"HD": {"VN": "1.3", "SO": "coordinate"}}

>>> seqA = Seq(None, length=20)
>>> queryA = SeqRecord(seqA, id="readA")
>>> sequences = [target1, queryA]
>>> coordinates = np.array([[4300, 4320], [0, 20]])
>>> alignment = Align.Alignment(sequences, coordinates)
>>> alignments.append(alignment)

>>> seqB = Seq(None, length=25)
>>> queryB = SeqRecord(seqB, id="readB")
>>> sequences = [target1, queryB]
>>> coordinates = np.array([[5900, 5925], [25, 0]])
>>> alignment = Align.Alignment(sequences, coordinates)
>>> alignments.append(alignment)

>>> seqC = Seq(None, length=40)
>>> queryC = SeqRecord(seqC, id="readC")
>>> sequences = [target2, queryC]
>>> coordinates = np.array([[12300, 12318], [0, 18]])
>>> alignment = Align.Alignment(sequences, coordinates)
>>> alignments.append(alignment)

>>> stream = StringIO()
>>> Align.write(alignments, stream, "sam")
3
>>> print(stream.getvalue())  
@HD VN:1.3  SO:coordinate
@SQ SN:chr1 LN:10000
@SQ SN:chr2 LN:15000
readA   0   chr1    4301    255 20M *   0   0   *   *
readB   16  chr1    5901    255 25M *   0   0   *   *
readC   0   chr2    12301   255 18M22S  *   0   0   *       *

浏览器扩展数据 (BED)

BED（浏览器扩展数据）文件通常用于存储基因转录本与基因组的比对结果。有关 BED 格式的完整说明，请参见UCSC 的说明。

BED 文件具有三个必填字段和九个可选字段。子目录 Tests/Blat 中的 bed12.bed 文件是具有 12 个字段的 BED 文件的示例。

chr22   1000    5000    mRNA1   960 +   1200    4900    255,0,0 2   567,488,    0,3512,
chr22   2000    6000    mRNA2   900 -   2300    5960    0,255,0 2   433,399,    0,3601,

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("bed12.bed", "bed")
>>> len(alignments)
2
>>> for alignment in alignments:
...     print(alignment.coordinates)
...
[[1000 1567 4512 5000]
 [   0  567  567 1055]]
[[2000 2433 5601 6000]
 [ 832  399  399    0]]

请注意，第一个序列（“mRNA1”）被映射到正向链，而第二个序列（“mRNA2”）被映射到反向链。

由于 BED 文件未存储每条染色体的长度，因此目标序列的长度被设置为其最大值。

>>> alignment.target
SeqRecord(seq=Seq(None, length=9223372036854775807), id='chr22', name='<unknown name>', description='', dbxrefs=[])

查询序列的长度可以从其比对信息中推断出来。

>>> alignment.query
SeqRecord(seq=Seq(None, length=832), id='mRNA2', name='<unknown name>', description='', dbxrefs=[])

比对分数（第 5 列）和存储在第 7-9 列中的信息（在 BED 格式规范中称为 thickStart、thickEnd 和 itemRgb）存储在 alignment 对象上的属性中。

>>> alignment.score
900.0
>>> alignment.thickStart
2300
>>> alignment.thickEnd
5960
>>> alignment.itemRgb
'0,255,0'

要以 BED 格式打印比对结果，您可以使用 Python 的内置 format 函数：

>>> print(format(alignment, "bed"))  
chr22   2000    6000    mRNA2   900 -   2300    5960    0,255,0 2   433,399,    0,3601,

或者，您可以使用 alignment 对象的 format 方法。这使您可以指定要写入的字段数作为 bedN 关键字参数：

>>> print(alignment.format("bed"))  
chr22   2000    6000    mRNA2   900 -   2300    5960    0,255,0 2   433,399,    0,3601,

>>> print(alignment.format("bed", 3))  
chr22   2000    6000

>>> print(alignment.format("bed", 6))  
chr22   2000    6000    mRNA2   900 -

相同的关键字参数可与 Align.write 一起使用：

>>> Align.write(alignments, "mybed3file.bed", "bed", bedN=3)
2
>>> Align.write(alignments, "mybed6file.bed", "bed", bedN=6)
2
>>> Align.write(alignments, "mybed12file.bed", "bed")
2

bigBed

bigBed 文件格式是 BED 文件浏览器扩展数据 (BED)的索引二进制版本。要创建 bigBed 文件，您可以使用 UCSC 的 bedToBigBed 程序（`) <https://genome.ucsc.edu/goldenPath/help/bigBed.html>`__. 或您可以通过调用 Bio.Align.write 函数并使用 fmt="bigbed" 来使用 Biopython 创建它。虽然这两种方法应该生成相同的 bigBed 文件，但使用 bedToBigBed 速度更快，可能更可靠，因为它是最权威的标准。由于 bigBed 文件带有内置索引，因此它允许您快速搜索特定基因组区域。

例如，让我们解析 bigBed 文件 dna_rna.bb，该文件位于 Biopython 分发版中的 Tests/Blat 子目录中。

>>> from Bio import Align
>>> alignments = Align.parse("dna_rna.bb", "bigbed")
>>> len(alignments)
4
>>> print(alignments.declaration)  
table bed
"Browser Extensible Data"
(
   string          chrom;          "Reference sequence chromosome or scaffold"
   uint            chromStart;     "Start position in chromosome"
   uint            chromEnd;       "End position in chromosome"
   string          name;           "Name of item."
   uint            score;          "Score (0-1000)"
   char[1]         strand;         "+ or - for strand"
   uint            thickStart;     "Start of where display should be thick (start codon)"
   uint            thickEnd;       "End of where display should be thick (stop codon)"
   uint            reserved;       "Used as itemRgb as of 2004-11-22"
   int             blockCount;     "Number of blocks"
   int[blockCount] blockSizes;     "Comma separated list of block sizes"
   int[blockCount] chromStarts;    "Start positions relative to chromStart"
)

该 declaration 包含用于创建 bigBed 文件的列规范（以 AutoSql 格式表示）。目标序列（通常是与序列比对的染色体）存储在 targets 属性中。在 bigBed 格式中，只存储每个目标的标识符和大小。在此示例中，只有一个染色体。

>>> alignments.targets
[SeqRecord(seq=Seq(None, length=198295559), id='chr3', name='<unknown name>', description='<unknown description>', dbxrefs=[])]

让我们看看单个比对。比对信息以与 BED 文件（请参见第浏览器扩展数据 (BED)节）相同的方式存储。

>>> alignment = next(alignments)
>>> alignment.target.id
'chr3'
>>> alignment.query.id
'NR_046654.1'
>>> alignment.coordinates
array([[42530895, 42530958, 42532020, 42532095, 42532563, 42532606],
       [     181,      118,      118,       43,       43,        0]])
>>> alignment.thickStart
42530895
>>> alignment.thickEnd
42532606
>>> print(alignment)  
chr3       42530895 ????????????????????????????????????????????????????????????
                  0 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
NR_046654       181 ????????????????????????????????????????????????????????????

chr3       42530955 ????????????????????????????????????????????????????????????
                 60 |||---------------------------------------------------------
NR_046654       121 ???---------------------------------------------------------
...
chr3       42532515 ????????????????????????????????????????????????????????????
               1620 ------------------------------------------------||||||||||||
NR_046654        43 ------------------------------------------------????????????

chr3       42532575 ??????????????????????????????? 42532606
               1680 |||||||||||||||||||||||||||||||     1711
NR_046654        31 ???????????????????????????????        0

默认 bigBed 格式不存储目标和查询的序列内容。如果这些内容在其他地方可用（例如，Fasta 文件），您可以设置 alignment.target.seq 和 alignment.query.seq，以便在打印比对结果时显示序列内容，或者以需要序列内容的格式（例如 Clustal，请参见第ClustalW节）写入比对结果。Biopython 分发版中 Tests 子目录中的 test_Align_bigbed.py 测试脚本提供了一些有关如何执行此操作的示例。

现在让我们看看如何搜索序列区域。这些是在 bigBed 文件中存储的序列，以 BED 格式打印（请参见第浏览器扩展数据 (BED)节）。

>>> for alignment in alignments:
...     print(format(alignment, "bed"))  
...
chr3    42530895    42532606    NR_046654.1 1000    -   42530895    42532606    0   3   63,75,43,   0,1125,1668,

chr3    42530895    42532606    NR_046654.1_modified    978 -   42530895    42532606    0   5   27,36,17,56,43, 0,27,1125,1144,1668,

chr3    48663767    48669174    NR_111921.1 1000    +   48663767    48669174    0   3   46,82,76,   0,1873,5331,

chr3    48663767    48669174    NR_111921.1_modified    972 +   48663767    48669174    0   5   28,17,76,6,76,  0,29,1873,1949,5331,

对 alignments 对象使用 search 方法，以在 chr3 上的 48000000 到 49000000 位置之间查找区域。此方法返回一个迭代器。

>>> selected_alignments = alignments.search("chr3", 48000000, 49000000)
>>> for alignment in selected_alignments:
...     print(alignment.query.id)
...
NR_111921.1
NR_111921.1_modified

染色体名称可以是 None，以包含所有染色体，并且开始和结束位置可以是 None，以分别从位置 0 开始搜索或继续搜索到染色体的末尾。

以 bigBed 格式写入比对结果与调用 Bio.Align.write 一样简单。

>>> Align.write(alignments, "output.bb", "bigbed")

您可以指定要包含在 bigBed 文件中的 BED 字段数。例如，要写入 BED6 文件，请使用：

>>> Align.write(alignments, "output.bb", "bigbed", bedN=6)

与 bedToBigBed 相同，您可以在 bigBed 输出中包含其他列。假设文件 bedExample2.as（位于 Biopython 分发版中 Tests/Blat 子目录中）存储了以 AutoSql 格式表示的包含的 BED 字段的声明。我们可以按如下方式读取此声明：

>>> from Bio.Align import bigbed
>>> with open("bedExample2.as") as stream:
...     autosql_data = stream.read()
...
>>> declaration = bigbed.AutoSQLTable.from_string(autosql_data)
>>> type(declaration)
<class 'Bio.Align.bigbed.AutoSQLTable'>
>>> print(declaration)
table hg18KGchr7
"UCSC Genes for chr7 with color plus GeneSymbol and SwissProtID"
(
   string  chrom;         "Reference sequence chromosome or scaffold"
   uint    chromStart;    "Start position of feature on chromosome"
   uint    chromEnd;      "End position of feature on chromosome"
   string  name;          "Name of gene"
   uint    score;         "Score"
   char[1] strand;        "+ or - for strand"
   uint    thickStart;    "Coding region start"
   uint    thickEnd;      "Coding region end"
   uint    reserved;      "Green on + strand, Red on - strand"
   string  geneSymbol;    "Gene Symbol"
   string  spID;          "SWISS-PROT protein Accession number"
)

现在，我们可以通过调用以下命令写入包含 9 个 BED 字段以及附加字段 geneSymbol 和 spID 的 bigBed 文件：

>>> Align.write(
...     alignments,
...     "output.bb",
...     "bigbed",
...     bedN=9,
...     declaration=declaration,
...     extraIndex=["name", "geneSymbol"],
... )

在这里，我们还要求在 bigBed 文件中包含对 name 和 geneSymbol 的附加索引。 Align.write 预计在 alignment.annotations 字典中找到键 geneSymbol 和 spID。请参阅 Biopython 分发版中 Tests 子目录中的 test_Align_bigbed.py 测试脚本，以获取有关以 bigBed 格式写入比对文件的更多示例。

可选参数是 compress（默认值为 True）、blockSize（默认值为 256）和 itemsPerSlot（默认值为 512）。有关这些参数的说明，请参见 UCSC 的 bedToBigBed 程序的文档。通过在创建 bigBed 文件时使用 compress=False 和 itemsPerSlot=1，可以加快对 bigBed 文件的搜索速度。

模式空间布局 (PSL)

PSL（模式空间布局）文件是由 BLAST 类比对工具 BLAT [Kent2002] 生成的。与 BED 文件（请参见第浏览器扩展数据 (BED)节）类似，PSL 文件通常用于存储转录本与基因组的比对结果。这是一个简短的 BLAT 文件的示例（位于 Biopython 分发版中 Tests/Blat 子目录中的 dna_rna.psl），标准 PSL 标题由 5 行组成。

psLayout version 3

match   mis-    rep.    N's Q gap   Q gap   T gap   T gap   strand  Q           Q       Q       Q   T           T       T       T   block   blockSizes  qStarts  tStarts
        match   match       count   bases   count   bases           name        size    start   end name        size    start   end count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
165 0   39  0   0   0   2   5203    +   NR_111921.1 216 0   204 chr3    198295559   48663767    48669174    3   46,82,76,   0,46,128,   48663767,48665640,48669098,
175 0   6   0   0   0   2   1530    -   NR_046654.1 181 0   181 chr3    198295559   42530895    42532606    3   63,75,43,   0,63,138,   42530895,42532020,42532563,
162 2   39  0   1   2   3   5204    +   NR_111921.1_modified    220 3   208 chr3    198295559   48663767    48669174    5   28,17,76,6,76,  3,31,48,126,132,    48663767,48663796,48665640,48665716,48669098,
172 1   6   0   1   3   3   1532    -   NR_046654.1_modified    190 3   185 chr3    198295559   42530895    42532606    5   27,36,17,56,43, 5,35,71,88,144, 42530895,42530922,42532020,42532039,42532563,

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("dna_rna.psl", "psl")
>>> alignments.metadata
{'psLayout version': '3'}

遍历比对结果，以获得每行一个 Alignment 对象。

>>> for alignment in alignments:
...     print(alignment.target.id, alignment.query.id)
...
chr3 NR_046654.1
chr3 NR_046654.1_modified
chr3 NR_111921.1
chr3 NR_111921.1_modified

让我们更详细地看看最后一个比对。PSL 文件中的前四列显示匹配次数、不匹配次数、与重复区域比对的核苷酸数以及与 N（未知）字符比对的核苷酸数。这些值存储为 Alignment 对象的属性。

>>> alignment.matches
162
>>> alignment.misMatches
2
>>> alignment.repMatches
39
>>> alignment.nCount
0

由于目标和查询的序列数据未在 PSL 文件中明确存储，因此 alignment.target 和 alignment.query 的序列内容未定义。但是，它们的序列长度是已知的。

>>> alignment.target  
SeqRecord(seq=Seq(None, length=198295559), id='chr3', ...)
>>> alignment.query  
SeqRecord(seq=Seq(None, length=220), id='NR_111921.1_modified', ...)

我们可以以 BED 或 PSL 格式打印比对结果。

>>> print(format(alignment, "bed"))  
chr3    48663767    48669174    NR_111921.1_modified    0   +   48663767    48669174    0   5   28,17,76,6,76,  0,29,1873,1949,5331,

>>> print(format(alignment, "psl"))  
162 2   39  0   1   2   3   5204    +   NR_111921.1_modified    220 3   208 chr3    198295559   48663767    48669174    5   28,17,76,6,76,  3,31,48,126,132,    48663767,48663796,48665640,48665716,48669098,

这里，匹配数、不匹配数、重复区域匹配数以及与未知核苷酸的匹配数来自 Alignment 对象的相应属性。如果这些属性不可用，例如比对结果不是来自 PSL 文件，则这些数字将使用序列内容（如果可用）进行计算。目标序列中的重复区域通过将序列掩蔽为小写或大写字符来表示，掩蔽方式由 mask 关键字参数的以下值定义。

False（默认）：不要单独计算与掩蔽序列的匹配。
"lower"：计算并报告与小写字符的匹配，作为与重复区域的匹配。
"upper"：计算并报告与大写字符的匹配，作为与重复区域的匹配。

用于未知核苷酸的字符由 wildcard 参数定义。为了与 BLAT 保持一致，默认情况下通配符为 "N"。如果您不想单独计算与任何未知核苷酸的匹配，请使用 wildcard=None。

>>> import numpy
>>> from Bio import Align
>>> query = "GGTGGGGG"
>>> target = "AAAAAAAggggGGNGAAAAA"
>>> coordinates = numpy.array([[0, 7, 15, 20], [0, 0, 8, 8]])
>>> alignment = Align.Alignment([target, query], coordinates)
>>> print(alignment)
target            0 AAAAAAAggggGGNGAAAAA 20
                  0 -------....||.|----- 20
query             0 -------GGTGGGGG-----  8

>>> line = alignment.format("psl")
>>> print(line)  
6   1   0   1   0   0   0   0   +   query   8   0   8   target   20   7   15   1   8,   0,   7,
>>> line = alignment.format("psl", mask="lower")
>>> print(line)  
3   1   3   1   0   0   0   0   +   query   8   0   8   target   20   7   15   1   8,   0,   7,
>>> line = alignment.format("psl", mask="lower", wildcard=None)
>>> print(line)  
3   2   3   0   0   0   0   0   +   query   8   0   8   target   20   7   15   1   8,   0,   7,

使用 Bio.Align.write 将比对结果写入 PSL 格式的输出文件时，可以使用相同的参数。此函数还有一个额外的关键字 header（默认值为 True），用于指定是否应写入 PSL 标头。

除了 format 方法之外，还可以使用 Python 内置的 format 函数。

>>> print(format(alignment, "psl"))  
6   1   0   1   0   0   0   0   +   query   8   0   8   target   20   7   15   1   8,   0,   7,

允许在 Python 中的格式化 (f-) 字符串中使用 Alignment 对象。

>>> line = f"The alignment in PSL format is '{alignment:psl}'."
>>> print(line)  
The alignment in PSL format is '6   1   0   1   0   0   0   0   +   query   8   0   8   target   20   7   15   1   8,   0,   7,
'

请注意，可选关键字参数不能与 format 函数或格式化字符串一起使用。

bigPsl

bigPsl 文件是一个 bigBed 文件，具有 BED12+13 格式，包括 12 个预定义的 BED 字段和 13 个自定义字段，这些字段在 UCSC 提供的 AutoSql 文件 bigPsl.as 中定义，创建了 PSL 文件的索引二进制版本（参见第模式空间布局 (PSL) 节）。要创建 bigPsl 文件，可以使用 UCSC 的 pslToBigPsl 和 bedToBigBed 程序，或者通过调用 Bio.Align.write 函数并使用 fmt="bigpsl" 来使用 Biopython。虽然这两种方法应该生成相同的 bigPsl 文件，但 UCSC 工具速度更快，并且可能更可靠，因为它是金标准。由于 bigPsl 文件是 bigBed 文件，因此它们带有内置索引，允许您快速搜索特定基因组区域。

例如，让我们解析 Biopython 发行版中 Tests/Blat 子目录中提供的 bigBed 文件 dna_rna.psl.bb。此文件是 bigBed 文件 dna_rna.bb（参见第 bigBed 节）和 PSL 文件 dna_rna.psl（参见第模式空间布局 (PSL) 节）的 bigPsl 等效文件。

>>> from Bio import Align
>>> alignments = Align.parse("dna_rna.psl.bb", "bigpsl")
>>> len(alignments)
4
>>> print(alignments.declaration)  
table bigPsl
"bigPsl pairwise alignment"
(
   string          chrom;           "Reference sequence chromosome or scaffold"
   uint            chromStart;      "Start position in chromosome"
   uint            chromEnd;        "End position in chromosome"
   string          name;            "Name or ID of item, ideally both human readable and unique"
   uint            score;           "Score (0-1000)"
   char[1]         strand;          "+ or - indicates whether the query aligns to the + or - strand on the reference"
   uint            thickStart;      "Start of where display should be thick (start codon)"
   uint            thickEnd;        "End of where display should be thick (stop codon)"
   uint            reserved;        "RGB value (use R,G,B string in input file)"
   int             blockCount;      "Number of blocks"
   int[blockCount] blockSizes;      "Comma separated list of block sizes"
   int[blockCount] chromStarts;     "Start positions relative to chromStart"
   uint            oChromStart;     "Start position in other chromosome"
   uint            oChromEnd;       "End position in other chromosome"
   char[1]         oStrand;         "+ or -, - means that psl was reversed into BED-compatible coordinates"
   uint            oChromSize;      "Size of other chromosome."
   int[blockCount] oChromStarts;    "Start positions relative to oChromStart or from oChromStart+oChromSize depending on strand"
   lstring         oSequence;       "Sequence on other chrom (or edit list, or empty)"
   string          oCDS;            "CDS in NCBI format"
   uint            chromSize;       "Size of target chromosome"
   uint            match;           "Number of bases matched."
   uint            misMatch;        "Number of bases that don't match"
   uint            repMatch;        "Number of bases that match but are part of repeats"
   uint            nCount;          "Number of 'N' bases"
   uint            seqType;         "0=empty, 1=nucleotide, 2=amino_acid"
)

声明包含由 UCSC 的 bigPsl.as AutoSql 文件定义的列规范。目标序列（通常是将序列与其比对的染色体）存储在 targets 属性中。在 bigBed 格式中，只存储每个目标的标识符和大小。在本例中，只有一个染色体。

>>> alignments.targets
[SeqRecord(seq=Seq(None, length=198295559), id='chr3', name='<unknown name>', description='<unknown description>', dbxrefs=[])]

迭代比对结果会为每行提供一个 Alignment 对象。

>>> for alignment in alignments:
...     print(alignment.target.id, alignment.query.id)
...
chr3 NR_046654.1
chr3 NR_046654.1_modified
chr3 NR_111921.1
chr3 NR_111921.1_modified

让我们看一下各个比对结果。比对信息以与相应 PSL 文件相同的方式存储（参见第模式空间布局 (PSL) 节）。

>>> alignment.coordinates
array([[48663767, 48663795, 48663796, 48663813, 48665640, 48665716,
        48665716, 48665722, 48669098, 48669174],
       [       3,       31,       31,       48,       48,      124,
             126,      132,      132,      208]])
>>> alignment.thickStart
48663767
>>> alignment.thickEnd
48669174
>>> alignment.matches
162
>>> alignment.misMatches
2
>>> alignment.repMatches
39
>>> alignment.nCount
0

我们可以以 BED 或 PSL 格式打印比对结果。

>>> print(format(alignment, "bed"))  
chr3    48663767    48669174    NR_111921.1_modified    1000    +   48663767    48669174    0   5   28,17,76,6,76,  0,29,1873,1949,5331,

>>> print(format(alignment, "psl"))  
162 2   39  0   1   2   3   5204    +   NR_111921.1_modified    220 3   208 chr3    198295559   48663767    48669174    5   28,17,76,6,76,  3,31,48,126,132,    48663767,48663796,48665640,48665716,48669098,

由于 bigPsl 文件是 bigBed 文件的特例，因此可以在比对结果对象上使用 search 方法来查找对特定基因组区域的比对结果。例如，我们可以查找 chr3 上位置 48000000 和 49000000 之间的区域。

>>> selected_alignments = alignments.search("chr3", 48000000, 49000000)
>>> for alignment in selected_alignments:
...     print(alignment.query.id)
...
NR_111921.1
NR_111921.1_modified

染色体名称可以是 None，以包含所有染色体，并且开始和结束位置可以是 None，以分别从位置 0 开始搜索或继续搜索到染色体的末尾。

要使用 Biopython 编写 bigPsl 文件，请使用 Bio.Align.write(alignments, "myfilename.bb", fmt="bigpsl")，其中 myfilename.bb 是输出 bigPsl 文件的名称。或者，可以使用（二进制）流进行输出。其他选项是：

compress：如果为 True（默认），则使用 zlib 压缩数据；如果为 False，则不压缩数据。
extraIndex：要索引的额外列名称的字符串列表。
cds：如果为 True，则查找类型为 CDS 的查询特征并在 PSL 文件中以 NCBI 风格写入（默认：False）。
fa：如果为 True，则将查询序列包含在 PSL 文件中（默认：False）。
mask：指定目标序列中的重复区域是否被掩蔽，并应在 repMatches 字段中报告，而不是在 matches 字段中报告。可接受的值为：
- None：不掩蔽（默认）；
- "lower"：用小写字符掩蔽；
- "upper"：用大写字符掩蔽。
wildcard：在 nCount 字段中报告目标或查询序列中通配符（代表未知核苷酸）的比对结果，而不是在 matches、misMatches 或 repMatches 字段中报告。默认值为 "N"。

有关如何获取匹配数、不匹配数、重复区域匹配数以及与未知核苷酸的匹配数的说明，请参见第模式空间布局 (PSL) 节。

其他可选参数是 blockSize（默认值为 256）和 itemsPerSlot（默认值为 512）。有关这些参数的描述，请参阅 UCSC 的 bedToBigBed 程序的文档。通过在创建 bigPsl 文件时使用 compress=False 和 itemsPerSlot=1，可以加快搜索 bigPsl 文件的速度。

多重比对格式 (MAF)

MAF（多重比对格式）文件以人类可读的格式存储一系列多重序列比对结果。MAF 文件通常用于存储基因组彼此比对的结果。Biopython 发行版中 Tests/MAF 子目录中的文件 ucsc_test.maf 是一个简单 MAF 文件的示例。

track name=euArc visibility=pack mafDot=off frames="multiz28wayFrames" speciesOrder="hg16 panTro1 baboon mm4 rn3" description="A sample alignment"
##maf version=1 scoring=tba.v8
# tba.v8 (((human chimp) baboon) (mouse rat))
# multiz.v7
# maf_project.v5 _tba_right.maf3 mouse _tba_C
# single_cov2.v4 single_cov2 /dev/stdin

a score=23262.0
s hg16.chr7    27578828 38 + 158545518 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s panTro1.chr6 28741140 38 + 161576975 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s baboon         116834 38 +   4622798 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
s mm4.chr6     53215344 38 + 151104725 -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
s rn3.chr4     81344243 40 + 187371129 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG

a score=5062.0
s hg16.chr7    27699739 6 + 158545518 TAAAGA
s panTro1.chr6 28862317 6 + 161576975 TAAAGA
s baboon         241163 6 +   4622798 TAAAGA
s mm4.chr6     53303881 6 + 151104725 TAAAGA
s rn3.chr4     81444246 6 + 187371129 taagga

a score=6636.0
s hg16.chr7    27707221 13 + 158545518 gcagctgaaaaca
s panTro1.chr6 28869787 13 + 161576975 gcagctgaaaaca
s baboon         249182 13 +   4622798 gcagctgaaaaca
s mm4.chr6     53310102 13 + 151104725 ACAGCTGAAAATA

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("ucsc_test.maf", "maf")

在文件头（轨道行和之后以“#”开头的行）中显示的信息存储在 alignments 对象的 metadata 属性中。

>>> alignments.metadata  
{'name': 'euArc',
 'visibility': 'pack',
 'mafDot': 'off',
 'frames': 'multiz28wayFrames',
 'speciesOrder': ['hg16', 'panTro1', 'baboon', 'mm4', 'rn3'],
 'description': 'A sample alignment',
 'MAF Version': '1',
 'Scoring': 'tba.v8',
 'Comments': ['tba.v8 (((human chimp) baboon) (mouse rat))',
              'multiz.v7',
              'maf_project.v5 _tba_right.maf3 mouse _tba_C',
              'single_cov2.v4 single_cov2 /dev/stdin']}

通过迭代 alignments，我们为 MAF 文件中的每个比对块获取一个 Alignment 对象。

>>> alignment = next(alignments)
>>> alignment.score
23262.0
>>> {seq.id: len(seq) for seq in alignment.sequences}  
{'hg16.chr7': 158545518,
 'panTro1.chr6': 161576975,
 'baboon': 4622798,
 'mm4.chr6': 151104725,
 'rn3.chr4': 187371129}
>>> print(alignment.coordinates)
[[27578828 27578829 27578831 27578831 27578850 27578850 27578866]
 [28741140 28741141 28741143 28741143 28741162 28741162 28741178]
 [  116834   116835   116837   116837   116856   116856   116872]
 [53215344 53215344 53215346 53215347 53215366 53215366 53215382]
 [81344243 81344243 81344245 81344245 81344264 81344267 81344283]]
>>> print(alignment)
hg16.chr7  27578828 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG 27578866
panTro1.c  28741140 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG 28741178
baboon       116834 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG   116872
mm4.chr6   53215344 -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG 53215382
rn3.chr4   81344243 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG 81344283

>>> print(format(alignment, "phylip"))
5 42
hg16.chr7 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
panTro1.chAAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
baboon    AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
mm4.chr6  -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
rn3.chr4  -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG

除了“a”（比对块）和“s”（序列）行之外，MAF 文件可能包含带有关于此块之前和之后的基因组序列信息的“i”行，带有关于比对结果中空部分信息的“e”行，以及显示每个比对碱基质量的“q”行。这是一个包含这些行的比对块的示例。

a score=19159.000000
s mm9.chr10                         3014644 45 + 129993255 CCTGTACC---CTTTGGTGAGAATTTTTGTTTCAGTGTTAAAAGTTTG
s hg18.chr6                        15870786 46 - 170899992 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT
i hg18.chr6                        I 9085 C 0
s panTro2.chr6                     16389355 46 - 173908612 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT
q panTro2.chr6                                             99999999999999999999999-9999999999999999999-9999
i panTro2.chr6                     I 9106 C 0
s calJac1.Contig6394                   6182 46 +    133105 CCTATACCTTTCTTTCATGAGAA-TTTTGTTTGAATCCTAAAC-TTTT
i calJac1.Contig6394               N 0 C 0
s loxAfr1.scaffold_75566               1167 34 -     10574 ------------TTTGGTTAGAA-TTATGCTTTAATTCAAAAC-TTCC
q loxAfr1.scaffold_75566                                   ------------99999699899-9999999999999869998-9997
i loxAfr1.scaffold_75566           N 0 C 0
e tupBel1.scaffold_114895.1-498454   167376 4145 -    498454 I
e echTel1.scaffold_288249             87661 7564 +    100002 I
e otoGar1.scaffold_334.1-359464      181217 2931 -    359464 I
e ponAbe2.chr6                     16161448 8044 - 174210431 I

这是文件 ucsc_mm9_chr10.maf（Biopython 发行版中 Tests/MAF 子目录中提供）中的第 10 个比对块。

>>> from Bio import Align
>>> alignments = Align.parse("ucsc_mm9_chr10.maf", "maf")
>>> for i in range(10):
...     alignment = next(alignments)
...
>>> alignment.score
19159.0
>>> print(alignment)
mm9.chr10   3014644 CCTGTACC---CTTTGGTGAGAATTTTTGTTTCAGTGTTAAAAGTTTG   3014689
hg18.chr6 155029206 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT 155029160
panTro2.c 157519257 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT 157519211
calJac1.C      6182 CCTATACCTTTCTTTCATGAGAA-TTTTGTTTGAATCCTAAAC-TTTT      6228
loxAfr1.s      9407 ------------TTTGGTTAGAA-TTATGCTTTAATTCAAAAC-TTCC      9373

“i”行显示当前比对块中的序列与前一个和后一个比对块中的序列之间的关系。此信息存储在相应序列的 annotations 属性中。

>>> alignment.sequences[0].annotations
{}
>>> alignment.sequences[1].annotations
{'leftStatus': 'I', 'leftCount': 9085, 'rightStatus': 'C', 'rightCount': 0}

显示在此块和前一个块之间插入了 9085 个碱基（“I”），而此块与后一个块是连续的（“C”）。有关这些字段和状态字符的完整说明，请参见 UCSC 文档。

“q”行显示序列质量，序列质量存储在相应序列的 annotations 属性的“quality”字典键下。

>>> alignment.sequences[2].annotations["quality"]
'9999999999999999999999999999999999999999999999'
>>> alignment.sequences[4].annotations["quality"]
'9999969989999999999999998699989997'

“e” 行显示有关物种的信息，这些物种在该比对块之前和之后具有连续序列，但在该比对块中没有比对的核苷酸。这些信息存储在 alignment.annotations 字典的“empty” 键下。

>>> alignment.annotations["empty"]  
[(SeqRecord(seq=Seq(None, length=498454), id='tupBel1.scaffold_114895.1-498454', name='', description='', dbxrefs=[]), (331078, 326933), 'I'),
 (SeqRecord(seq=Seq(None, length=100002), id='echTel1.scaffold_288249', name='', description='', dbxrefs=[]), (87661, 95225), 'I'),
 (SeqRecord(seq=Seq(None, length=359464), id='otoGar1.scaffold_334.1-359464', name='', description='', dbxrefs=[]), (178247, 175316), 'I'),
 (SeqRecord(seq=Seq(None, length=174210431), id='ponAbe2.chr6', name='', description='', dbxrefs=[]), (158048983, 158040939), 'I')]

例如，这表明在 ponAbe2.chr6 基因组序列的相反链上，从位置 158040939 到 158048983 插入了一些非比对碱基（“I”）。再次，有关“e” 行的完整定义，请参阅 UCSC 文档。

要以 MAF 格式打印比对，可以使用 Python 的内置 format 函数。

>>> print(format(alignment, "MAF"))
a score=19159.000000
s mm9.chr10                         3014644   45 + 129993255 CCTGTACC---CTTTGGTGAGAATTTTTGTTTCAGTGTTAAAAGTTTG
s hg18.chr6                        15870786   46 - 170899992 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT
i hg18.chr6                        I 9085 C 0
s panTro2.chr6                     16389355   46 - 173908612 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT
q panTro2.chr6                                               99999999999999999999999-9999999999999999999-9999
i panTro2.chr6                     I 9106 C 0
s calJac1.Contig6394                   6182   46 +    133105 CCTATACCTTTCTTTCATGAGAA-TTTTGTTTGAATCCTAAAC-TTTT
i calJac1.Contig6394               N 0 C 0
s loxAfr1.scaffold_75566               1167   34 -     10574 ------------TTTGGTTAGAA-TTATGCTTTAATTCAAAAC-TTCC
q loxAfr1.scaffold_75566                                     ------------99999699899-9999999999999869998-9997
i loxAfr1.scaffold_75566           N 0 C 0
e tupBel1.scaffold_114895.1-498454   167376 4145 -    498454 I
e echTel1.scaffold_288249             87661 7564 +    100002 I
e otoGar1.scaffold_334.1-359464      181217 2931 -    359464 I
e ponAbe2.chr6                     16161448 8044 - 174210431 I

要写入完整的 MAF 文件，请使用 Bio.Align.write(alignments, "myfilename.maf", fmt="maf")，其中 myfilename.maf 是输出 MAF 文件的名称。或者，可以使用（文本）流进行输出。文件头信息将从 alignments 对象的 metadata 属性获取。如果要从头开始创建比对，可以使用 Alignments（复数）类来创建一个类似列表的 alignments 对象（请参见第 Alignments 类节）并为其指定一个 metadata 属性。

bigMaf

bigMaf 文件是一个 bigBed 文件，具有 BED3+1 格式，包含 3 个必需的 BED 字段以及一个自定义字段，该字段将 MAF 比对块作为字符串存储，从而创建 MAF 文件的索引二进制版本（请参见第多重比对格式 (MAF) 节）。相关 AutoSql 文件 bigMaf.as 由 UCSC 提供。要创建 bigMaf 文件，可以使用 UCSC 的 mafToBigMaf 和 bedToBigBed 程序，或者可以使用 Biopython 通过使用 fmt="bigmaf" 调用 Bio.Align.write 函数。虽然这两种方法应该生成相同的 bigMaf 文件，但 UCSC 工具的速度要快得多，并且可能更可靠，因为它是黄金标准。由于 bigMaf 文件是 bigBed 文件，因此它们自带内置索引，使您能够快速搜索参考基因组的特定区域。

Biopython 分发版 Tests/MAF 子目录中的 ucsc_test.bb 文件是 bigMaf 文件的示例。该文件等效于 MAF 文件 ucsc_test.maf（请参见第多重比对格式 (MAF) 节）。要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("ucsc_test.bb", "bigmaf")
>>> len(alignments)
3
>>> print(alignments.declaration)  
table bedMaf
"Bed3 with MAF block"
(
   string  chrom;         "Reference sequence chromosome or scaffold"
   uint    chromStart;    "Start position in chromosome"
   uint    chromEnd;      "End position in chromosome"
   lstring mafBlock;      "MAF block"
)

声明包含由 UCSC 的 bigMaf.as AutoSql 文件定义的列规范。

bigMaf 文件不存储 MAF 文件中找到的标题信息，但它确实定义了一个参考基因组。对应的 SeqRecord 存储在 alignments 对象的 targets 属性中。

>>> alignments.reference
'hg16'
>>> alignments.targets  
[SeqRecord(seq=Seq(None, length=158545518), id='hg16.chr7', ...)]

通过遍历 alignments，我们为 bigMaf 文件中的每个比对块获得一个 Alignment 对象。

>>> alignment = next(alignments)
>>> alignment.score
23262.0
>>> {seq.id: len(seq) for seq in alignment.sequences}  
{'hg16.chr7': 158545518,
 'panTro1.chr6': 161576975,
 'baboon': 4622798,
 'mm4.chr6': 151104725,
 'rn3.chr4': 187371129}
>>> print(alignment.coordinates)
[[27578828 27578829 27578831 27578831 27578850 27578850 27578866]
 [28741140 28741141 28741143 28741143 28741162 28741162 28741178]
 [  116834   116835   116837   116837   116856   116856   116872]
 [53215344 53215344 53215346 53215347 53215366 53215366 53215382]
 [81344243 81344243 81344245 81344245 81344264 81344267 81344283]]
>>> print(alignment)
hg16.chr7  27578828 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG 27578866
panTro1.c  28741140 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG 28741178
baboon       116834 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG   116872
mm4.chr6   53215344 -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG 53215382
rn3.chr4   81344243 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG 81344283

>>> print(format(alignment, "phylip"))
5 42
hg16.chr7 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
panTro1.chAAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
baboon    AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
mm4.chr6  -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
rn3.chr4  -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG

“i”、“e” 和 “q” 行中的信息存储方式与相应的 MAF 文件相同（请参见第多重比对格式 (MAF) 节）。

>>> from Bio import Align
>>> alignments = Align.parse("ucsc_mm9_chr10.bb", "bigmaf")
>>> for i in range(10):
...     alignment = next(alignments)
...
>>> alignment.score
19159.0
>>> print(alignment)
mm9.chr10   3014644 CCTGTACC---CTTTGGTGAGAATTTTTGTTTCAGTGTTAAAAGTTTG   3014689
hg18.chr6 155029206 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT 155029160
panTro2.c 157519257 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT 157519211
calJac1.C      6182 CCTATACCTTTCTTTCATGAGAA-TTTTGTTTGAATCCTAAAC-TTTT      6228
loxAfr1.s      9407 ------------TTTGGTTAGAA-TTATGCTTTAATTCAAAAC-TTCC      9373

>>> print(format(alignment, "MAF"))
a score=19159.000000
s mm9.chr10                         3014644   45 + 129993255 CCTGTACC---CTTTGGTGAGAATTTTTGTTTCAGTGTTAAAAGTTTG
s hg18.chr6                        15870786   46 - 170899992 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT
i hg18.chr6                        I 9085 C 0
s panTro2.chr6                     16389355   46 - 173908612 CCTATACCTTTCTTTTATGAGAA-TTTTGTTTTAATCCTAAAC-TTTT
q panTro2.chr6                                               99999999999999999999999-9999999999999999999-9999
i panTro2.chr6                     I 9106 C 0
s calJac1.Contig6394                   6182   46 +    133105 CCTATACCTTTCTTTCATGAGAA-TTTTGTTTGAATCCTAAAC-TTTT
i calJac1.Contig6394               N 0 C 0
s loxAfr1.scaffold_75566               1167   34 -     10574 ------------TTTGGTTAGAA-TTATGCTTTAATTCAAAAC-TTCC
q loxAfr1.scaffold_75566                                     ------------99999699899-9999999999999869998-9997
i loxAfr1.scaffold_75566           N 0 C 0
e tupBel1.scaffold_114895.1-498454   167376 4145 -    498454 I
e echTel1.scaffold_288249             87661 7564 +    100002 I
e otoGar1.scaffold_334.1-359464      181217 2931 -    359464 I
e ponAbe2.chr6                     16161448 8044 - 174210431 I


>>> alignment.sequences[1].annotations
{'leftStatus': 'I', 'leftCount': 9085, 'rightStatus': 'C', 'rightCount': 0}
>>> alignment.sequences[2].annotations["quality"]
'9999999999999999999999999999999999999999999999'
>>> alignment.sequences[4].annotations["quality"]
'9999969989999999999999998699989997'
>>> alignment.annotations["empty"]  
[(SeqRecord(seq=Seq(None, length=498454), id='tupBel1.scaffold_114895.1-498454', name='', description='', dbxrefs=[]), (331078, 326933), 'I'),
 (SeqRecord(seq=Seq(None, length=100002), id='echTel1.scaffold_288249', name='', description='', dbxrefs=[]), (87661, 95225), 'I'),
 (SeqRecord(seq=Seq(None, length=359464), id='otoGar1.scaffold_334.1-359464', name='', description='', dbxrefs=[]), (178247, 175316), 'I'),
 (SeqRecord(seq=Seq(None, length=174210431), id='ponAbe2.chr6', name='', description='', dbxrefs=[]), (158048983, 158040939), 'I')]

要写入完整的 bigMaf 文件，请使用 Bio.Align.write(alignments, "myfilename.bb", fmt="bigMaf")，其中 myfilename.bb 是输出 bigMaf 文件的名称。或者，可以使用（二进制）流进行输出。如果要从头开始创建比对，可以使用 Alignments（复数）类来创建一个类似列表的 alignments 对象（请参见第 Alignments 类节）并为其指定一个 targets 属性。后者必须是参考物种染色体 SeqRecord 对象的列表，其顺序与它们在比对中出现的顺序相同。或者，您可以在调用 Bio.Align.write 时使用 targets 关键字参数。每个 SeqRecord 的 id 必须采用 reference.chromosome 的形式，其中 reference 指的是参考物种。 Bio.Align.write 具有额外的关键字参数 compress（默认值为 True），用于指定是否应使用 zlib 压缩数据。其他可选参数是 blockSize（默认值为 256）和 itemsPerSlot（默认值为 512）。请参阅 UCSC 的 bedToBigBed 程序的文档以了解这些参数的描述。

由于 bigMaf 文件是 bigBed 文件的特例，因此可以使用 alignments 对象上的 search 方法查找与参考物种特定区域的比对。例如，我们可以查找染色体 10 上位置 3018000 和 3019000 之间的区域。

>>> selected_alignments = alignments.search("mm9.chr10", 3018000, 3019000)
>>> for alignment in selected_alignments:
...     start, end = alignment.coordinates[0, 0], alignment.coordinates[0, -1]
...     print(start, end)
...
3017743 3018161
3018161 3018230
3018230 3018359
3018359 3018482
3018482 3018644
3018644 3018822
3018822 3018932
3018932 3019271

染色体名称可以为 None 以包括所有染色体，起始位置和结束位置可以为 None 以分别从位置 0 开始搜索或继续搜索到染色体末尾。请注意，我们只能对参考物种的基因组位置进行搜索。

通过在创建 bigMaf 文件时使用 compress=False 和 itemsPerSlot=1，可以加快 bigMaf 文件的搜索速度。

UCSC 链文件格式

链文件描述两个核苷酸序列之间的成对比对，允许两个序列中出现间隙。链文件中只存储每个比对子序列的长度和间隙长度；序列本身不存储。链文件通常用于存储两个基因组组装版本之间的比对，允许将一个基因组组装版本上的比对提升到另一个基因组组装版本。这是一个链文件示例（在 Biopython 分发版的 Tests/Blat 子目录中作为 psl_34_001.chain 提供）。

chain 16 chr4 191154276 + 61646095 61646111 hg18_dna 33 + 11 27 1
16
chain 33 chr1 249250621 + 10271783 10271816 hg18_dna 33 + 0 33 2
33
chain 17 chr2 243199373 + 53575980 53575997 hg18_dna 33 - 8 25 3
17
chain 35 chr9 141213431 + 85737865 85737906 hg19_dna 50 + 9 50 4
41
chain 41 chr8 146364022 + 95160479 95160520 hg19_dna 50 + 8 49 5
41
chain 30 chr22 51304566 + 42144400 42144436 hg19_dna 50 + 11 47 6
36
chain 41 chr2 243199373 + 183925984 183926028 hg19_dna 50 + 1 49 7
6       0       4
38
chain 31 chr19 59128983 + 35483340 35483510 hg19_dna 50 + 10 46 8
25      134     0
11
chain 39 chr18 78077248 + 23891310 23891349 hg19_dna 50 + 10 49 9
39
...

该文件是通过对 PSL 文件 psl_34_001.psl 运行 UCSC 的 pslToChain 程序生成的。根据链文件格式规范，每个链块之后应该有一行空行，但一些工具（包括 pslToChain）显然没有遵循此规则。

要解析此文件，请使用

>>> from Bio import Align
>>> alignments = Align.parse("psl_34_001.chain", "chain")

遍历比对以获得每个链的一个 Alignment 对象。

>>> for alignment in alignments:
...     print(alignment.target.id, alignment.query.id)  
...
chr4 hg18_dna
chr1 hg18_dna
chr2 hg18_dna
chr9 hg19_dna
chr8 hg19_dna
chr22 hg19_dna
chr2 hg19_dna
...
chr1 hg19_dna

从起始位置开始遍历，直到到达第七个比对。

>>> alignments = iter(alignments)
>>> for i in range(7):
...     alignment = next(alignments)
...

检查比对分数和链 ID（每个链块的标题行中的第一个和最后一个数字），以确认我们获得了第七个比对。

>>> alignment.score
41.0
>>> alignment.annotations["id"]
'7'

我们可以以链文件格式打印比对。比对坐标与链块中的信息一致，包含一个 6 个核苷酸的比对部分、一个 4 个核苷酸的间隙以及一个 38 个核苷酸的比对部分。

>>> print(format(alignment, "chain"))  
chain 41 chr2 243199373 + 183925984 183926028 hg19_dna 50 + 1 49 7
6   0   4
38


>>> alignment.coordinates
array([[183925984, 183925990, 183925990, 183926028],
       [        1,         7,        11,        49]])
>>> print(alignment)
chr2      183925984 ??????----?????????????????????????????????????? 183926028
                  0 ||||||----||||||||||||||||||||||||||||||||||||||        48
hg19_dna          1 ????????????????????????????????????????????????        49

我们还可以以几种其他比对文件格式打印比对。

>>> print(format(alignment, "BED"))  
chr2    183925984   183926028   hg19_dna    41  +   183925984   183926028   0   2   6,38,   0,6,

>>> print(format(alignment, "PSL"))  
44  0   0   0   1   4   0   0   +   hg19_dna    50  1   49  chr2    243199373   183925984   183926028   2   6,38,   1,11,   183925984,183925990,

>>> print(format(alignment, "exonerate"))
vulgar: hg19_dna 1 49 + chr2 183925984 183926028 + 41 M 6 6 G 4 0 M 38 38

>>> print(alignment.format("exonerate", "cigar"))
cigar: hg19_dna 1 49 + chr2 183925984 183926028 + 41 M 6 I 4 M 38

>>> print(format(alignment, "sam"))  
hg19_dna    0   chr2    183925985   255 1S6M4I38M1S *   0   0   *   *   AS:i:41 id:A:7