成对序列比对

成对序列比对是通过优化两个序列之间的相似度得分来将两个序列彼此比对的过程。 Bio.Align 模块包含 PairwiseAligner 类，用于使用 Needleman-Wunsch、Smith-Waterman、Gotoh（三态）和 Waterman-Smith-Beyer 全局和局部成对比对算法进行全局和局部比对，并具有许多选项来更改比对参数。我们参考 Durbin 等人 [Durbin1998] 以获取有关序列比对算法的深入信息。

基本用法

要生成成对比对，首先创建一个 PairwiseAligner 对象

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()

PairwiseAligner 对象 aligner（参见第成对比对器对象节）存储用于成对比对的比对参数。这些属性可以在对象的构造函数中设置

>>> aligner = Align.PairwiseAligner(match_score=1.0)

或者在创建对象后设置

>>> aligner.match_score = 1.0

使用 aligner.score 方法计算两个序列之间的比对得分

>>> target = "GAACT"
>>> query = "GAT"
>>> score = aligner.score(target, query)
>>> score
3.0

aligner.align 方法返回 Alignment 对象，每个对象代表两个序列之间的一个比对

>>> alignments = aligner.align(target, query)
>>> alignment = alignments[0]
>>> alignment  
<Alignment object (2 rows x 5 columns) at ...>

遍历 Alignment 对象并打印它们以查看比对

>>> for alignment in alignments:
...     print(alignment)
...
target            0 GAACT 5
                  0 ||--| 5
query             0 GA--T 3

target            0 GAACT 5
                  0 |-|-| 5
query             0 G-A-T 3

每个比对都存储比对得分

>>> alignment.score
3.0

以及对已比对序列的指针

>>> alignment.target
'GAACT'
>>> alignment.query
'GAT'

在内部，比对以序列坐标的形式存储

>>> alignment = alignments[0]
>>> alignment.coordinates
array([[0, 2, 4, 5],
       [0, 2, 2, 3]])

这里，两行分别指目标序列和查询序列。这些坐标表明比对由以下三个块组成

target[0:2] 比对到 query[0:2]；
target[2:4] 比对到一个间隙，因为 query[2:2] 是一个空字符串；
target[4:5] 比对到 query[2:3]。

成对比对的已比对序列数量始终为 2

>>> len(alignment)
2

比对长度定义为比对中列的数量（如打印所示）。这等于匹配次数、不匹配次数以及目标和查询中间隙的总长度之和

>>> alignment.length
5

aligned 属性返回已比对子序列的起始和结束索引，对于第一个比对返回两个长度为 2 的元组

>>> alignment.aligned
array([[[0, 2],
        [4, 5]],

       [[0, 2],
        [2, 3]]])

而对于备选比对，则返回两个长度为 3 的元组

>>> alignment = alignments[1]
>>> print(alignment)
target            0 GAACT 5
                  0 |-|-| 5
query             0 G-A-T 3

>>> alignment.aligned
array([[[0, 1],
        [2, 3],
        [4, 5]],

       [[0, 1],
        [1, 2],
        [2, 3]]])

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

>>> aligner.mode = "global"
>>> aligner.mismatch_score = -10
>>> alignments = aligner.align("AAACAAA", "AAAGAAA")
>>> len(alignments)
2
>>> print(alignments[0])
target            0 AAAC-AAA 7
                  0 |||--||| 8
query             0 AAA-GAAA 7

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

       [[0, 3],
        [4, 7]]])
>>> print(alignments[1])
target            0 AAA-CAAA 7
                  0 |||--||| 8
query             0 AAAG-AAA 7

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

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

map 方法可以应用于成对比对 alignment1 以找到 alignment2 的查询与 alignment1 的目标的成对比对，其中 alignment2 的目标与 alignment1 的查询相同。一个典型的例子是， alignment1 是染色体和转录本之间的成对比对， alignment2 是转录本和序列（例如 RNA-seq 阅读）之间的成对比对，我们希望找到序列与染色体的比对

>>> aligner.mode = "local"
>>> aligner.open_gap_score = -1
>>> aligner.extend_gap_score = 0
>>> chromosome = "AAAAAAAACCCCCCCAAAAAAAAAAAGGGGGGAAAAAAAA"
>>> transcript = "CCCCCCCGGGGGG"
>>> alignments1 = aligner.align(chromosome, transcript)
>>> len(alignments1)
1
>>> alignment1 = alignments1[0]
>>> print(alignment1)
target            8 CCCCCCCAAAAAAAAAAAGGGGGG 32
                  0 |||||||-----------|||||| 24
query             0 CCCCCCC-----------GGGGGG 13

>>> sequence = "CCCCGGGG"
>>> alignments2 = aligner.align(transcript, sequence)
>>> len(alignments2)
1
>>> alignment2 = alignments2[0]
>>> print(alignment2)
target            3 CCCCGGGG 11
                  0 ||||||||  8
query             0 CCCCGGGG  8

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

>>> format(mapped_alignment, "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'

比对映射不依赖于序列内容。如果我们删除序列内容，则在 PSL 格式中会找到相同的比对（尽管我们显然失去了打印序列比对的能力）

>>> from Bio.Seq import Seq
>>> alignment1.target = Seq(None, len(alignment1.target))
>>> alignment1.query = Seq(None, len(alignment1.query))
>>> alignment2.target = Seq(None, len(alignment2.target))
>>> alignment2.query = Seq(None, len(alignment2.query))
>>> mapped_alignment = alignment1.map(alignment2)
>>> format(mapped_alignment, "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'

默认情况下，执行全局成对比对，这将在 target 和 query 的整个长度上找到最佳比对。相反，局部比对将找到具有最高比对得分的 target 和 query 的子序列。通过将 aligner.mode 设置为 "local" 可以生成局部比对

>>> aligner.mode = "local"
>>> target = "AGAACTC"
>>> query = "GAACT"
>>> score = aligner.score(target, query)
>>> score
5.0
>>> alignments = aligner.align(target, query)
>>> for alignment in alignments:
...     print(alignment)
...
target            1 GAACT 6
                  0 ||||| 5
query             0 GAACT 5

请注意，如果可以向比对添加得分 0 的片段，则最佳局部比对的定义存在一些歧义。我们遵循 Waterman 和 Eggert 的建议 [Waterman1987] 并禁止此类扩展。

成对比对器对象

PairwiseAligner 对象存储用于成对比对的所有比对参数。要查看所有参数值的概述，请使用

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner(match_score=1.0, mode="local")
>>> print(aligner)
Pairwise sequence aligner with parameters
  wildcard: None
  match_score: 1.000000
  mismatch_score: 0.000000
  target_internal_open_gap_score: 0.000000
  target_internal_extend_gap_score: 0.000000
  target_left_open_gap_score: 0.000000
  target_left_extend_gap_score: 0.000000
  target_right_open_gap_score: 0.000000
  target_right_extend_gap_score: 0.000000
  query_internal_open_gap_score: 0.000000
  query_internal_extend_gap_score: 0.000000
  query_left_open_gap_score: 0.000000
  query_left_extend_gap_score: 0.000000
  query_right_open_gap_score: 0.000000
  query_right_extend_gap_score: 0.000000
  mode: local

有关这些参数的定义，请参见下面的第替换分数节、第仿射间隙分数节和第一般间隙分数节。属性 mode（如上面第基本用法节中所述）可以设置为 "global" 或 "local"，分别指定全局或局部成对比对。

根据间隙评分参数（参见第仿射间隙分数节和第一般间隙分数节）和模式， PairwiseAligner 对象会自动选择适当的算法用于成对序列比对。要验证选择的算法，请使用

>>> aligner.algorithm
'Smith-Waterman'

此属性为只读属性。

PairwiseAligner 对象还存储比对期间使用的精度 \(\epsilon\)。 \(\epsilon\) 的值存储在属性 aligner.epsilon 中，默认值为 \(10^{-6}\)

>>> aligner.epsilon
1e-06

如果两个分数之间的绝对差小于 \(\epsilon\)，则这两个分数将被视为彼此相等，用于比对目的。

替换分数

替换分数定义了当两个字母（核苷酸或氨基酸）相互比对时添加到总得分中的值。 PairwiseAligner 要使用的替换分数可以通过两种方式指定

通过为相同的字母指定匹配分数，为不匹配的字母指定不匹配分数。核苷酸序列比对通常基于匹配分数和不匹配分数。例如，默认情况下，BLAST [Altschul1990] 使用 \(+1\) 的匹配分数和 \(-2\) 的不匹配分数进行 megablast 的核苷酸比对，间隙罚分为 2.5（有关间隙分数的更多信息，请参见第仿射间隙分数节）。可以通过设置 PairwiseAligner 对象的 match 和 mismatch 属性来指定匹配分数和不匹配分数
```
>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> aligner.match_score
1.0
>>> aligner.mismatch_score
0.0
>>> score = aligner.score("ACGT", "ACAT")
>>> print(score)
3.0
>>> aligner.match_score = 1.0
>>> aligner.mismatch_score = -2.0
>>> aligner.gap_score = -2.5
>>> score = aligner.score("ACGT", "ACAT")
>>> print(score)
1.0
```
使用匹配分数和不匹配分数时，可以指定一个通配符字符（默认情况下为 None）用于未知字母。这些字母在比对中将获得零分，无论匹配分数或不匹配分数的值如何
```
>>> aligner.wildcard = "?"
>>> score = aligner.score("ACGT", "AC?T")
>>> print(score)
3.0
```

或者，可以使用 PairwiseAligner 对象的 substitution_matrix 属性来指定替换矩阵。这使你能够对不同的匹配字母和不匹配字母对应用不同的分数。这通常用于氨基酸序列比对。例如，默认情况下，BLAST [Altschul1990] 使用 BLOSUM62 替换矩阵进行 blastp 的蛋白质比对。此替换矩阵可从 Biopython 获取

>>> from Bio.Align import substitution_matrices
>>> substitution_matrices.load()  
['BENNER22', 'BENNER6', 'BENNER74', 'BLASTN', 'BLASTP', 'BLOSUM45', 'BLOSUM50', 'BLOSUM62', ..., 'TRANS']
>>> matrix = substitution_matrices.load("BLOSUM62")
>>> print(matrix)  
#  Matrix made by matblas from blosum62.iij
...
     A    R    N    D    C    Q ...
A  4.0 -1.0 -2.0 -2.0  0.0 -1.0 ...
R -1.0  5.0  0.0 -2.0 -3.0  1.0 ...
N -2.0  0.0  6.0  1.0 -3.0  0.0 ...
D -2.0 -2.0  1.0  6.0 -3.0  0.0 ...
C  0.0 -3.0 -3.0 -3.0  9.0 -3.0 ...
Q -1.0  1.0  0.0  0.0 -3.0  5.0 ...
...
>>> aligner.substitution_matrix = matrix
>>> score = aligner.score("ACDQ", "ACDQ")
>>> score
24.0
>>> score = aligner.score("ACDQ", "ACNQ")
>>> score
19.0

使用替换矩阵时， X 不被解释为未知字符。相反，将使用替换矩阵提供的分数

>>> matrix["D", "X"]
-1.0
>>> score = aligner.score("ACDQ", "ACXQ")
>>> score
17.0

默认情况下，aligner.substitution_matrix 为 None。如果 aligner.substitution_matrix 不为 None，则属性 aligner.match_score 和 aligner.mismatch_score 将被忽略。将 aligner.match_score 或 aligner.mismatch_score 设置为有效值将重置 aligner.substitution_matrix 为 None。

仿射间隙评分

仿射间隙评分由一个打开间隙的得分和一个扩展现有间隙的得分定义。

\(\textrm{间隙评分} = \textrm{打开间隙评分} + (n-1) \times \textrm{扩展间隙评分}\),

其中 \(n\) 是间隙的长度。Biopython 的成对序列比对器允许通过指定 PairwiseAligner 对象的以下 12 个属性来对间隙评分方案进行精细控制。

打开评分	扩展评分
`query_left_open_gap_score`	`query_left_extend_gap_score`
`query_internal_open_gap_score`	`query_internal_extend_gap_score`
`query_right_open_gap_score`	`query_right_extend_gap_score`
`target_left_open_gap_score`	`target_left_extend_gap_score`
`target_internal_open_gap_score`	`target_internal_extend_gap_score`
`target_right_open_gap_score`	`target_right_extend_gap_score`

这些属性允许对内部间隙以及序列两端使用不同的间隙评分，如以下示例所示。

target	query	score
A		query 左端打开间隙评分
C		query 左端扩展间隙评分
C		query 左端扩展间隙评分
G	G	匹配评分
G	T	不匹配评分
G		query 内部打开间隙评分
A		query 内部扩展间隙评分
A		query 内部扩展间隙评分
T	T	匹配评分
A	A	匹配评分
G		query 内部打开间隙评分
C	C	匹配评分
	C	target 内部打开间隙评分
	C	target 内部扩展间隙评分
C	C	匹配评分
T	G	不匹配评分
C	C	匹配评分
	C	target 内部打开间隙评分
A	A	匹配评分
	T	target 右端打开间隙评分
	A	target 右端扩展间隙评分
	A	target 右端扩展间隙评分

为了方便起见，PairwiseAligner 对象具有其他属性，这些属性共同引用这些值中的多个值，如表成对比对器对象的元属性。中所示（按层次结构）。

表 1 成对比对器对象的元属性。
元属性	它映射到的属性
`gap_score`	`target_gap_score`, `query_gap_score`
`open_gap_score`	`target_open_gap_score`, `query_open_gap_score`
`extend_gap_score`	`target_extend_gap_score`, `query_extend_gap_score`
`internal_gap_score`	`target_internal_gap_score`, `query_internal_gap_score`
`internal_open_gap_score`	`target_internal_open_gap_score`, `query_internal_open_gap_score`
`internal_extend_gap_score`	`target_internal_extend_gap_score`, `query_internal_extend_gap_score`
`end_gap_score`	`target_end_gap_score`, `query_end_gap_score`
`end_open_gap_score`	`target_end_open_gap_score`, `query_end_open_gap_score`
`end_extend_gap_score`	`target_end_extend_gap_score`, `query_end_extend_gap_score`
`left_gap_score`	`target_left_gap_score`, `query_left_gap_score`
`right_gap_score`	`target_right_gap_score`, `query_right_gap_score`
`left_open_gap_score`	`target_left_open_gap_score`, `query_left_open_gap_score`
`left_extend_gap_score`	`target_left_extend_gap_score`, `query_left_extend_gap_score`
`right_open_gap_score`	`target_right_open_gap_score`, `query_right_open_gap_score`
`right_extend_gap_score`	`target_right_extend_gap_score`, `query_right_extend_gap_score`
`target_open_gap_score`	`target_internal_open_gap_score`, `target_left_open_gap_score`, `target_right_open_gap_score`
`target_extend_gap_score`	`target_internal_extend_gap_score`, `target_left_extend_gap_score`, `target_right_extend_gap_score`
`target_gap_score`	`target_open_gap_score`, `target_extend_gap_score`
`query_open_gap_score`	`query_internal_open_gap_score`, `query_left_open_gap_score`, `query_right_open_gap_score`
`query_extend_gap_score`	`query_internal_extend_gap_score`, `query_left_extend_gap_score`, `query_right_extend_gap_score`
`query_gap_score`	`query_open_gap_score`, `query_extend_gap_score`
`target_internal_gap_score`	`target_internal_open_gap_score`, `target_internal_extend_gap_score`
`target_end_gap_score`	`target_end_open_gap_score`, `target_end_extend_gap_score`
`target_end_open_gap_score`	`target_left_open_gap_score`, `target_right_open_gap_score`
`target_end_extend_gap_score`	`target_left_extend_gap_score`, `target_right_extend_gap_score`
`target_left_gap_score`	`target_left_open_gap_score`, `target_left_extend_gap_score`
`target_right_gap_score`	`target_right_open_gap_score`, `target_right_extend_gap_score`
`query_end_gap_score`	`query_end_open_gap_score`, `query_end_extend_gap_score`
`query_end_open_gap_score`	`query_left_open_gap_score`, `query_right_open_gap_score`
`query_end_extend_gap_score`	`query_left_extend_gap_score`, `query_right_extend_gap_score`
`query_internal_gap_score`	`query_internal_open_gap_score`, `query_internal_extend_gap_score`
`query_left_gap_score`	`query_left_open_gap_score`, `query_left_extend_gap_score`
`query_right_gap_score`	`query_right_open_gap_score`, `query_right_extend_gap_score`

通用间隙评分

为了对间隙评分进行更精细的控制，您可以指定间隙评分函数。例如，以下间隙评分函数不允许在查询序列中的两个核苷酸之后出现间隙。

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> def my_gap_score_function(start, length):
...     if start == 2:
...         return -1000
...     else:
...         return -1 * length
...
>>> aligner.query_gap_score = my_gap_score_function
>>> alignments = aligner.align("AACTT", "AATT")
>>> for alignment in alignments:
...     print(alignment)
...
target            0 AACTT 5
                  0 -|.|| 5
query             0 -AATT 4

target            0 AACTT 5
                  0 |-.|| 5
query             0 A-ATT 4

target            0 AACTT 5
                  0 ||.-| 5
query             0 AAT-T 4

target            0 AACTT 5
                  0 ||.|- 5
query             0 AATT- 4

使用预定义的替换矩阵和间隙评分

默认情况下，PairwiseAligner 对象的初始化使用 +1.0 的匹配评分、0.0 的不匹配评分以及所有间隙评分等于 0.0。虽然这在简单评分方案方面具有优势，但通常不会提供最佳性能。相反，您可以使用参数 scoring 在初始化 PairwiseAligner 对象时选择预定义的评分方案。目前，提供的评分方案包括 blastn 和 megablast，适用于核苷酸比对，以及 blastp，适用于蛋白质比对。选择这些评分方案将分别将 PairwiseAligner 对象初始化为 BLASTN、MegaBLAST 和 BLASTP 使用的默认评分参数。

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner(scoring="blastn")
>>> print(aligner)  
Pairwise sequence aligner with parameters
  substitution_matrix: <Array object at ...>
  target_internal_open_gap_score: -7.000000
  target_internal_extend_gap_score: -2.000000
  target_left_open_gap_score: -7.000000
  target_left_extend_gap_score: -2.000000
  target_right_open_gap_score: -7.000000
  target_right_extend_gap_score: -2.000000
  query_internal_open_gap_score: -7.000000
  query_internal_extend_gap_score: -2.000000
  query_left_open_gap_score: -7.000000
  query_left_extend_gap_score: -2.000000
  query_right_open_gap_score: -7.000000
  query_right_extend_gap_score: -2.000000
  mode: global

>>> print(aligner.substitution_matrix[:, :])
     A    T    G    C    S    W    R    Y    K    M    B    V    H    D    N
A  2.0 -3.0 -3.0 -3.0 -3.0 -1.0 -1.0 -3.0 -3.0 -1.0 -3.0 -1.0 -1.0 -1.0 -2.0
T -3.0  2.0 -3.0 -3.0 -3.0 -1.0 -3.0 -1.0 -1.0 -3.0 -1.0 -3.0 -1.0 -1.0 -2.0
G -3.0 -3.0  2.0 -3.0 -1.0 -3.0 -1.0 -3.0 -1.0 -3.0 -1.0 -1.0 -3.0 -1.0 -2.0
C -3.0 -3.0 -3.0  2.0 -1.0 -3.0 -3.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -3.0 -2.0
S -3.0 -3.0 -1.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
W -1.0 -1.0 -3.0 -3.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
R -1.0 -3.0 -1.0 -3.0 -1.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
Y -3.0 -1.0 -3.0 -1.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
K -3.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -2.0
M -1.0 -3.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
B -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
V -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
H -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
D -1.0 -1.0 -1.0 -3.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0
N -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0

迭代比对

由 aligner.align 返回的 alignments 是一种不可变可迭代对象（类似于 range）。虽然它们看起来类似于 tuple 或 list 中的 Alignment 对象，但它们的不同之处在于每个 Alignment 对象都是在需要时动态创建的。选择这种方法是因为比对的数量可能非常大，特别是对于较差的比对（请参阅第示例节，了解示例）。

您可以对 alignments 执行以下操作。

len(alignments) 返回存储的比对数量。此函数返回速度很快，即使比对数量非常大。如果比对数量非常大（通常，大于 9,223,372,036,854,775,807，这是 64 位机器上可以存储为 long int 的最大整数），len(alignments) 将引发 OverflowError。大量的比对表明比对质量较低。
```
>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> alignments = aligner.align("AAA", "AA")
>>> len(alignments)
3
```

您可以通过索引提取特定的比对。

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> alignments = aligner.align("AAA", "AA")
>>> print(alignments[2])
target            0 AAA 3
                  0 -|| 3
query             0 -AA 2

>>> print(alignments[0])
target            0 AAA 3
                  0 ||- 3
query             0 AA- 2

您可以迭代比对，例如以下示例所示。
```
>>> for alignment in alignments:
...     print(alignment)
...
```
可以将 alignments 迭代器转换为 list 或 tuple
```
>>> alignments = list(alignments)
```
在尝试调用 list(alignments) 将所有比对结果保存为列表之前，最好通过调用 len(alignments) 检查比对结果的数量。
比对分数（在 alignments 中每个比对结果的值都相同）存储为一个属性。这使您可以在继续提取单个比对结果之前检查比对分数。
```
>>> print(alignments.score)
2.0
```

比对到反向链

默认情况下，成对比对器将查询序列的正向链比对到目标序列的正向链。要计算 query 与 target 的反向链的比对分数，请使用 strand="-"

>>> from Bio import Align
>>> from Bio.Seq import reverse_complement
>>> target = "AAAACCC"
>>> query = "AACC"
>>> aligner = Align.PairwiseAligner()
>>> aligner.mismatch_score = -1
>>> aligner.internal_gap_score = -1
>>> aligner.score(target, query)  # strand is "+" by default
4.0
>>> aligner.score(target, reverse_complement(query), strand="-")
4.0
>>> aligner.score(target, query, strand="-")
0.0
>>> aligner.score(target, reverse_complement(query))
0.0

可以通过在调用 aligner.align 时指定 strand="-" 来获取针对反向链的比对结果。

>>> alignments = aligner.align(target, query)
>>> len(alignments)
1
>>> print(alignments[0])
target            0 AAAACCC 7
                  0 --||||- 7
query             0 --AACC- 4

>>> print(alignments[0].format("bed"))  
target   2   6   query   4   +   2   6   0   1   4,   0,

>>> alignments = aligner.align(target, reverse_complement(query), strand="-")
>>> len(alignments)
1
>>> print(alignments[0])
target            0 AAAACCC 7
                  0 --||||- 7
query             4 --AACC- 0

>>> print(alignments[0].format("bed"))  
target   2   6   query   4   -   2   6   0   1   4,   0,

>>> alignments = aligner.align(target, query, strand="-")
>>> len(alignments)
2
>>> print(alignments[0])
target            0 AAAACCC----  7
                  0 ----------- 11
query             4 -------GGTT  0

>>> print(alignments[1])
target            0 ----AAAACCC  7
                  0 ----------- 11
query             4 GGTT-------  0

请注意，如果左右间隙分数不同，将 query 比对到 target 的反向链的分数可能与将 query 的反向互补序列比对到 target 的正向链的分数不同。

>>> aligner.left_gap_score = -0.5
>>> aligner.right_gap_score = -0.2
>>> aligner.score(target, query)
2.8
>>> alignments = aligner.align(target, query)
>>> len(alignments)
1
>>> print(alignments[0])
target            0 AAAACCC 7
                  0 --||||- 7
query             0 --AACC- 4

>>> aligner.score(target, reverse_complement(query), strand="-")
3.1
>>> alignments = aligner.align(target, reverse_complement(query), strand="-")
>>> len(alignments)
1
>>> print(alignments[0])
target            0 AAAACCC 7
                  0 --||||- 7
query             4 --AACC- 0

替换矩阵

替换矩阵 [Durbin1998] 提供了用于对两个不同残基相互替代的可能性进行分类的评分项。这在进行序列比较时至关重要。Biopython 提供了大量常见的替换矩阵，包括著名的 PAM 和 BLOSUM 矩阵系列，并提供了创建自己的替换矩阵的功能。

数组对象

您可以将替换矩阵视为二维数组，其中索引是字母（核苷酸或氨基酸）而不是整数。 Bio.Align.substitution_matrices 中的 Array 类是 numpy 数组的子类，支持通过整数和特定字符串进行索引。 Array 实例可以是一维数组或方形二维数组。例如，一维 Array 对象可用于存储 DNA 序列的核苷酸频率，而二维 Array 对象可用于表示序列比对的评分矩阵。

要创建一个一维 Array，只需要指定允许字母的字母表。

>>> from Bio.Align.substitution_matrices import Array
>>> counts = Array("ACGT")
>>> print(counts)
A 0.0
C 0.0
G 0.0
T 0.0

允许的字母存储在 alphabet 属性中。

>>> counts.alphabet
'ACGT'

此属性是只读的；修改底层的 _alphabet 属性可能会导致意外结果。元素可以通过字母和整数索引进行访问。

>>> counts["C"] = -3
>>> counts[2] = 7
>>> print(counts)
A  0.0
C -3.0
G  7.0
T  0.0

>>> counts[1]
-3.0

使用字母表中不存在的字母或超出范围的索引会导致 IndexError。

>>> counts["U"]
Traceback (most recent call last):
    ...
IndexError: 'U'
>>> counts["X"] = 6
Traceback (most recent call last):
    ...
IndexError: 'X'
>>> counts[7]
Traceback (most recent call last):
    ...
IndexError: index 7 is out of bounds for axis 0 with size 4

可以通过指定 dims=2 创建二维 Array。

>>> from Bio.Align.substitution_matrices import Array
>>> counts = Array("ACGT", dims=2)
>>> print(counts)
    A   C   G   T
A 0.0 0.0 0.0 0.0
C 0.0 0.0 0.0 0.0
G 0.0 0.0 0.0 0.0
T 0.0 0.0 0.0 0.0

同样，字母和整数都可以用于索引，并且指定字母表中不存在的字母会导致 IndexError。

>>> counts["A", "C"] = 12.0
>>> counts[2, 1] = 5.0
>>> counts[3, "T"] = -2
>>> print(counts)
    A    C   G    T
A 0.0 12.0 0.0  0.0
C 0.0  0.0 0.0  0.0
G 0.0  5.0 0.0  0.0
T 0.0  0.0 0.0 -2.0

>>> counts["X", 1]
Traceback (most recent call last):
    ...
IndexError: 'X'
>>> counts["A", 5]
Traceback (most recent call last):
    ...
IndexError: index 5 is out of bounds for axis 1 with size 4

从二维数组中选择一行或一列将返回一维 Array。

>>> counts = Array("ACGT", dims=2)
>>> counts["A", "C"] = 12.0
>>> counts[2, 1] = 5.0
>>> counts[3, "T"] = -2

>>> counts["G"]
Array([0., 5., 0., 0.],
      alphabet='ACGT')
>>> counts[:, "C"]
Array([12.,  0.,  5.,  0.],
      alphabet='ACGT')

Array 对象因此可以用作数组和字典。它们可以转换为纯 numpy 数组或纯字典对象。

>>> import numpy as np
>>> x = Array("ACGT")
>>> x["C"] = 5

>>> x
Array([0., 5., 0., 0.],
      alphabet='ACGT')
>>> a = np.array(x)  # create a plain numpy array
>>> a
array([0., 5., 0., 0.])
>>> d = dict(x)  # create a plain dictionary
>>> d
{'A': 0.0, 'C': 5.0, 'G': 0.0, 'T': 0.0}

虽然 Array 的字母表通常是一个字符串，但您也可以使用 (不可变) 对象的元组。例如，这用于密码子替换矩阵（如后面显示的 substitution_matrices.load("SCHNEIDER") 示例），其中键不是单个核苷酸或氨基酸，而是三核苷酸密码子。

虽然 Array 的 alphabet 属性是不可变的，但您可以通过从字母表中选择您感兴趣的字母来创建一个新的 Array 对象。例如，

>>> a = Array("ABCD", dims=2, data=np.arange(16).reshape(4, 4))
>>> print(a)
     A    B    C    D
A  0.0  1.0  2.0  3.0
B  4.0  5.0  6.0  7.0
C  8.0  9.0 10.0 11.0
D 12.0 13.0 14.0 15.0

>>> b = a.select("CAD")
>>> print(b)
     C    A    D
C 10.0  8.0 11.0
A  2.0  0.0  3.0
D 14.0 12.0 15.0

请注意，这也允许您重新排序字母表。

字母表中找不到的字母的数据被设置为零。

>>> c = a.select("DEC")
>>> print(c)
     D   E    C
D 15.0 0.0 14.0
E  0.0 0.0  0.0
C 11.0 0.0 10.0

由于 Array 类是 numpy 数组的子类，因此可以将其用作 numpy 数组。如果出现在数学运算中的 Array 对象具有不同的字母表，则会触发 ValueError，例如

>>> from Bio.Align.substitution_matrices import Array
>>> d = Array("ACGT")
>>> r = Array("ACGU")
>>> d + r
Traceback (most recent call last):
    ...
ValueError: alphabets are inconsistent

从成对序列比对中计算替换矩阵

由于 Array 是 numpy 数组的子类，因此您可以像对待 numpy 数组一样对 Array 对象进行数学运算。在这里，我们将通过从 *大肠杆菌* 和 *枯草芽孢杆菌* 的 16S 核糖体 RNA 基因序列的比对中计算评分矩阵来说明这一点。首先，我们创建一个 PairwiseAligner 对象（参见第成对序列比对章），并使用 blastn 使用的默认分数对其进行初始化。

>>> from Bio.Align import PairwiseAligner
>>> aligner = PairwiseAligner(scoring="blastn")
>>> aligner.mode = "local"

接下来，我们读取 *大肠杆菌* 和 *枯草芽孢杆菌* 的 16S 核糖体 RNA 基因序列（在 Biopython 发行版的 Tests/Align/ecoli.fa 和 Tests/Align/bsubtilis.fa 中提供），并将它们相互比对。

>>> from Bio import SeqIO
>>> sequence1 = SeqIO.read("ecoli.fa", "fasta")
>>> sequence2 = SeqIO.read("bsubtilis.fa", "fasta")
>>> alignments = aligner.align(sequence1, sequence2)

生成的比对结果数量非常多。

>>> len(alignments)
1990656

但是，由于它们之间只存在微不足道的差异，我们任意选择了第一个比对结果，并计算了每个替换的次数。

>>> alignment = alignments[0]
>>> substitutions = alignment.substitutions
>>> print(substitutions)
      A     C     G     T
A 307.0  19.0  34.0  19.0
C  15.0 280.0  25.0  29.0
G  34.0  24.0 401.0  20.0
T  24.0  36.0  20.0 228.0

我们对总数进行归一化，以找到每个替换的概率，并创建一个观察频率的对称矩阵。

>>> observed_frequencies = substitutions / substitutions.sum()
>>> observed_frequencies = (observed_frequencies + observed_frequencies.transpose()) / 2.0
>>> print(format(observed_frequencies, "%.4f"))
       A      C      G      T
A 0.2026 0.0112 0.0224 0.0142
C 0.0112 0.1848 0.0162 0.0215
G 0.0224 0.0162 0.2647 0.0132
T 0.0142 0.0215 0.0132 0.1505

背景概率是在每个序列中单独找到 A、C、G 或 T 核苷酸的概率。这可以计算为每行或每列的总和。

>>> background = observed_frequencies.sum(0)
>>> print(format(background, "%.4f"))
A 0.2505
C 0.2337
G 0.3165
T 0.1993

随机预期替换的次数仅仅是背景分布的自乘积。

>>> expected_frequencies = background[:, None].dot(background[None, :])
>>> print(format(expected_frequencies, "%.4f"))
       A      C      G      T
A 0.0627 0.0585 0.0793 0.0499
C 0.0585 0.0546 0.0740 0.0466
G 0.0793 0.0740 0.1002 0.0631
T 0.0499 0.0466 0.0631 0.0397

然后，可以将评分矩阵计算为观察概率与预期概率的几率比的对数。

>>> oddsratios = observed_frequencies / expected_frequencies
>>> import numpy as np
>>> scoring_matrix = np.log2(oddsratios)
>>> print(scoring_matrix)
     A    C    G    T
A  1.7 -2.4 -1.8 -1.8
C -2.4  1.8 -2.2 -1.1
G -1.8 -2.2  1.4 -2.3
T -1.8 -1.1 -2.3  1.9

该矩阵可以用作成对比对器的替换矩阵（参见第成对序列比对章）。

>>> aligner.substitution_matrix = scoring_matrix

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

在此示例中，我们将首先从 Clustalw 文件 protein.aln（也在线这里提供）中读取蛋白质序列比对。

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

第 ClustalW 节包含有关执行此操作的更多信息。

比对的 substitutions 属性存储了不同残基相互替代的次数。

>>> substitutions = alignment.substitutions

为了使示例更易读，我们将仅选择具有极性带电侧链的氨基酸。

>>> substitutions = substitutions.select("DEHKR")
>>> print(substitutions)
       D      E      H      K      R
D 2360.0  270.0   15.0    1.0   48.0
E  241.0 3305.0   15.0   45.0    2.0
H    0.0   18.0 1235.0    8.0    0.0
K    0.0    9.0   24.0 3218.0  130.0
R    2.0    2.0   17.0  103.0 2079.0

矩阵中其他氨基酸的行和列被删除。

接下来，我们对矩阵进行归一化并使其对称。

>>> observed_frequencies = substitutions / substitutions.sum()
>>> observed_frequencies = (observed_frequencies + observed_frequencies.transpose()) / 2.0
>>> print(format(observed_frequencies, "%.4f"))
       D      E      H      K      R
D 0.1795 0.0194 0.0006 0.0000 0.0019
E 0.0194 0.2514 0.0013 0.0021 0.0002
H 0.0006 0.0013 0.0939 0.0012 0.0006
K 0.0000 0.0021 0.0012 0.2448 0.0089
R 0.0019 0.0002 0.0006 0.0089 0.1581

对行或列求和给出每个残基的相对出现频率。

>>> background = observed_frequencies.sum(0)
>>> print(format(background, "%.4f"))
D 0.2015
E 0.2743
H 0.0976
K 0.2569
R 0.1697

>>> sum(background) == 1.0
True

然后，残基对的预期频率为

>>> expected_frequencies = background[:, None].dot(background[None, :])
>>> print(format(expected_frequencies, "%.4f"))
       D      E      H      K      R
D 0.0406 0.0553 0.0197 0.0518 0.0342
E 0.0553 0.0752 0.0268 0.0705 0.0465
H 0.0197 0.0268 0.0095 0.0251 0.0166
K 0.0518 0.0705 0.0251 0.0660 0.0436
R 0.0342 0.0465 0.0166 0.0436 0.0288

这里，background[:, None] 创建了一个二维数组，该数组包含一个包含 expected_frequencies 值的列，而 rxpected_frequencies[None, :] 则是将这些值作为一行的一个二维数组。对它们的点积（内积）创建一个预期频率矩阵，其中每个条目包含两个相乘的 expected_frequencies 值。例如，expected_frequencies['D', 'E'] 等于 residue_frequencies['D'] * residue_frequencies['E']。

现在，我们可以通过将观察频率除以预期频率并取对数来计算对数几率矩阵。

>>> import numpy as np
>>> scoring_matrix = np.log2(observed_frequencies / expected_frequencies)
>>> print(scoring_matrix)
      D    E    H     K    R
D   2.1 -1.5 -5.1 -10.4 -4.2
E  -1.5  1.7 -4.4  -5.1 -8.3
H  -5.1 -4.4  3.3  -4.4 -4.7
K -10.4 -5.1 -4.4   1.9 -2.3
R  -4.2 -8.3 -4.7  -2.3  2.5

此矩阵可以用作执行比对时的替换矩阵。例如，

>>> from Bio.Align import PairwiseAligner
>>> aligner = PairwiseAligner()
>>> aligner.substitution_matrix = scoring_matrix
>>> aligner.gap_score = -3.0
>>> alignments = aligner.align("DEHEK", "DHHKK")
>>> print(alignments[0])
target            0 DEHEK 5
                  0 |.|.| 5
query             0 DHHKK 5

>>> print("%.2f" % alignments.score)
-2.18
>>> score = (
...     scoring_matrix["D", "D"]
...     + scoring_matrix["E", "H"]
...     + scoring_matrix["H", "H"]
...     + scoring_matrix["E", "K"]
...     + scoring_matrix["K", "K"]
... )
>>> print("%.2f" % score)
-2.18

（有关详细信息，请参阅第成对序列比对章）。

从文件中读取 `Array` 对象

Bio.Align.substitution_matrices 包含一个解析器，用于从文件中读取一维和二维 Array 对象。一维数组用简单的两列格式表示，第一列包含键，第二列包含相应的值。例如，文件 hg38.chrom.sizes（从 UCSC 获取），在 Biopython 发行版的 Tests/Align 子目录中提供，包含人类基因组组装 hg38 中每条染色体的核苷酸大小。

chr1    248956422
chr2    242193529
chr3    198295559
chr4    190214555
...
chrUn_KI270385v1    990
chrUn_KI270423v1    981
chrUn_KI270392v1    971
chrUn_KI270394v1    970

要解析此文件，请使用

>>> from Bio.Align import substitution_matrices
>>> with open("hg38.chrom.sizes") as handle:
...     table = substitution_matrices.read(handle)
...
>>> print(table)  
chr1 248956422.0
chr2 242193529.0
chr3 198295559.0
chr4 190214555.0
...
chrUn_KI270423v1       981.0
chrUn_KI270392v1       971.0
chrUn_KI270394v1       970.0

使用 dtype=int 将值读取为整数。

>>> with open("hg38.chrom.sizes") as handle:
...     table = substitution_matrices.read(handle, int)
...
>>> print(table)  
chr1 248956422
chr2 242193529
chr3 198295559
chr4 190214555
...
chrUn_KI270423v1       981
chrUn_KI270392v1       971
chrUn_KI270394v1       970

对于二维数组，我们遵循 NCBI 提供的替换矩阵的文件格式。例如，BLOSUM62 矩阵是 NCBI 的蛋白质-蛋白质 BLAST [Altschul1990] 程序 blastp 的默认替换矩阵，存储如下。

#  Matrix made by matblas from blosum62.iij
#  * column uses minimum score
#  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
#  Blocks Database = /data/blocks_5.0/blocks.dat
#  Cluster Percentage: >= 62
#  Entropy =   0.6979, Expected =  -0.5209
   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4
...

此文件包含在 Biopython 发行版下的 Bio/Align/substitution_matrices/data 中。要解析此文件，请使用

>>> from Bio.Align import substitution_matrices
>>> with open("BLOSUM62") as handle:
...     matrix = substitution_matrices.read(handle)
...
>>> print(matrix.alphabet)
ARNDCQEGHILKMFPSTWYVBZX*
>>> print(matrix["A", "D"])
-2.0

以 # 开头的标题行存储在属性 header 中。

>>> matrix.header[0]
'Matrix made by matblas from blosum62.iij'

现在，我们可以将此矩阵用作比对器对象上的替换矩阵。

>>> from Bio.Align import PairwiseAligner
>>> aligner = PairwiseAligner()
>>> aligner.substitution_matrix = matrix

要保存一个 Array 对象，请先创建一个字符串

>>> text = str(matrix)
>>> print(text)  
#  Matrix made by matblas from blosum62.iij
#  * column uses minimum score
#  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
#  Blocks Database = /data/blocks_5.0/blocks.dat
#  Cluster Percentage: >= 62
#  Entropy =   0.6979, Expected =  -0.5209
     A    R    N    D    C    Q    E    G    H    I    L    K    M    F    P    S ...
A  4.0 -1.0 -2.0 -2.0  0.0 -1.0 -1.0  0.0 -2.0 -1.0 -1.0 -1.0 -1.0 -2.0 -1.0  1.0 ...
R -1.0  5.0  0.0 -2.0 -3.0  1.0  0.0 -2.0  0.0 -3.0 -2.0  2.0 -1.0 -3.0 -2.0 -1.0 ...
N -2.0  0.0  6.0  1.0 -3.0  0.0  0.0  0.0  1.0 -3.0 -3.0  0.0 -2.0 -3.0 -2.0  1.0 ...
D -2.0 -2.0  1.0  6.0 -3.0  0.0  2.0 -1.0 -1.0 -3.0 -4.0 -1.0 -3.0 -3.0 -1.0  0.0 ...
C  0.0 -3.0 -3.0 -3.0  9.0 -3.0 -4.0 -3.0 -3.0 -1.0 -1.0 -3.0 -1.0 -2.0 -3.0 -1.0 ...
...

并将 text 写入文件。

加载预定义的替换矩阵

Biopython 包含文献中定义的大量替换矩阵，包括 BLOSUM（块替换矩阵）[Henikoff1992] 和 PAM（点接受突变）矩阵 [Dayhoff1978]。这些矩阵以平面文件的形式存在于 Bio/Align/substitution_matrices/data 目录中，可以使用 substitution_matrices 子模块中的 load 函数加载到 Python 中。例如，BLOSUM62 矩阵可以通过运行以下代码加载：

>>> from Bio.Align import substitution_matrices
>>> m = substitution_matrices.load("BLOSUM62")

该替换矩阵的字母表包含遗传密码中使用的 20 个氨基酸，三个模糊的氨基酸 B（天冬酰胺或天冬氨酸）、Z（谷氨酰胺或谷氨酸）和 X（代表任何氨基酸），以及用星号表示的终止密码子

>>> m.alphabet
'ARNDCQEGHILKMFPSTWYVBZX*'

要获取可用替换矩阵的完整列表，请在不带参数的情况下使用 load

>>> substitution_matrices.load()  
['BENNER22', 'BENNER6', 'BENNER74', 'BLASTN', 'BLASTP', 'BLOSUM45', 'BLOSUM50', ..., 'TRANS']

请注意，Schneider 等人提供的替换矩阵 [Schneider2005] 使用包含三个核苷酸密码子的字母表

>>> m = substitution_matrices.load("SCHNEIDER")
>>> m.alphabet  
('AAA', 'AAC', 'AAG', 'AAT', 'ACA', 'ACC', 'ACG', 'ACT', ..., 'TTG', 'TTT')

示例

假设您想对上面相同的两个血红蛋白序列（HBA_HUMAN、HBB_HUMAN）进行全局成对比对，这两个序列存储在 alpha.faa 和 beta.faa 中

>>> from Bio import Align
>>> from Bio import SeqIO
>>> seq1 = SeqIO.read("alpha.faa", "fasta")
>>> seq2 = SeqIO.read("beta.faa", "fasta")
>>> aligner = Align.PairwiseAligner()
>>> score = aligner.score(seq1.seq, seq2.seq)
>>> print(score)
72.0

显示对齐得分 72.0。要查看单个对齐，请执行以下操作：

>>> alignments = aligner.align(seq1.seq, seq2.seq)

在这个例子中，最佳对齐的总数是巨大的（超过 \(4 \times 10^{37}\)），调用 len(alignments) 将引发 OverflowError

>>> len(alignments)
Traceback (most recent call last):
...
OverflowError: number of optimal alignments is larger than 9223372036854775807

让我们看看第一个对齐

>>> alignment = alignments[0]

对齐对象存储对齐得分以及对齐本身

>>> print(alignment.score)
72.0
>>> print(alignment)
target            0 MV-LS-PAD--KTN--VK-AA-WGKV-----GAHAGEYGAEALE-RMFLSF----P-TTK
                  0 ||-|--|----|----|--|--||||-----|---||--|--|--|--|------|-|--
query             0 MVHL-TP--EEK--SAV-TA-LWGKVNVDEVG---GE--A--L-GR--L--LVVYPWT--

target           41 TY--FPHF----DLSHGS---AQVK-G------HGKKV--A--DA-LTNAVAHV-DDMPN
                 60 ----|--|----|||------|-|--|------|||||--|--|--|--|--|--|---|
query            39 --QRF--FESFGDLS---TPDA-V-MGNPKVKAHGKKVLGAFSD-GL--A--H-LD---N

target           79 ALS----A-LSD-LHAH--KLR-VDPV-NFK-LLSHC---LLVT--LAAHLPA----EFT
                120 -|-----|-||--||----||--|||--||--||------|-|---||-|-------|||
query            81 -L-KGTFATLS-ELH--CDKL-HVDP-ENF-RLL---GNVL-V-CVLA-H---HFGKEFT

target          119 PA-VH-ASLDKFLAS---VSTV------LTS--KYR- 142
                180 |--|--|------|----|--|------|----||-- 217
query           124 P-PV-QA------A-YQKV--VAGVANAL--AHKY-H 147

通过惩罚间隙通常可以获得更好的对齐：更高的间隙打开成本和更低的间隙扩展成本。对于氨基酸序列匹配分数，通常编码在像 PAM 或 BLOSUM 这样的矩阵中。因此，我们可以通过使用 BLOSUM62 矩阵以及 10 的间隙打开惩罚和 0.5 的间隙扩展惩罚来获得我们示例中更有意义的对齐

>>> from Bio import Align
>>> from Bio import SeqIO
>>> from Bio.Align import substitution_matrices
>>> seq1 = SeqIO.read("alpha.faa", "fasta")
>>> seq2 = SeqIO.read("beta.faa", "fasta")
>>> aligner = Align.PairwiseAligner()
>>> aligner.open_gap_score = -10
>>> aligner.extend_gap_score = -0.5
>>> aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
>>> score = aligner.score(seq1.seq, seq2.seq)
>>> print(score)
292.5
>>> alignments = aligner.align(seq1.seq, seq2.seq)
>>> len(alignments)
2
>>> print(alignments[0].score)
292.5
>>> print(alignments[0])
target            0 MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF-DLS-----HGS
                  0 ||-|.|..|..|.|.||||--...|.|.|||.|.....|.|...|..|-|||-----.|.
query             0 MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGN

target           53 AQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAH
                 60 ..||.|||||..|.....||.|........||.||..||.|||.||.||...|...||.|
query            58 PKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHH

target          113 LPAEFTPAVHASLDKFLASVSTVLTSKYR 142
                120 ...||||.|.|...|..|.|...|..||. 149
query           118 FGKEFTPPVQAAYQKVVAGVANALAHKYH 147

此对齐与我们之前使用相同的序列和相同的参数使用 EMBOSS needle 获得的对齐具有相同的得分。

要执行局部对齐，请将 aligner.mode 设置为 'local'

>>> aligner.mode = "local"
>>> aligner.open_gap_score = -10
>>> aligner.extend_gap_score = -1
>>> alignments = aligner.align("LSPADKTNVKAA", "PEEKSAV")
>>> print(len(alignments))
1
>>> alignment = alignments[0]
>>> print(alignment)
target            2 PADKTNV 9
                  0 |..|..| 7
query             0 PEEKSAV 7

>>> print(alignment.score)
16.0

广义成对比对

在大多数情况下，PairwiseAligner 用于执行由单字母核苷酸或氨基酸组成的序列（字符串或 Seq 对象）的比对。更一般地说，PairwiseAligner 也可以应用于任意对象的列表或元组。本节将描述此类广义成对比对的一些示例。

使用替换矩阵和字母表的广义成对比对

Schneider 等人 [Schneider2005] 创建了一个用于对齐三个核苷酸密码子的替换矩阵（有关更多信息，请参见第下节替换矩阵）。此替换矩阵与包含所有三个字母密码子的字母表相关联

>>> from Bio.Align import substitution_matrices
>>> m = substitution_matrices.load("SCHNEIDER")
>>> m.alphabet  
('AAA', 'AAC', 'AAG', 'AAT', 'ACA', 'ACC', 'ACG', 'ACT', ..., 'TTG', 'TTT')

我们可以使用此矩阵将密码子序列相互比对

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> aligner.substitution_matrix = m
>>> aligner.gap_score = -1.0
>>> s1 = ("AAT", "CTG", "TTT", "TTT")
>>> s2 = ("AAT", "TTA", "TTT")
>>> alignments = aligner.align(s1, s2)
>>> len(alignments)
2
>>> print(alignments[0])
AAT CTG TTT TTT
||| ... ||| ---
AAT TTA TTT ---

>>> print(alignments[1])
AAT CTG TTT TTT
||| ... --- |||
AAT TTA --- TTT

请注意，在这个例子中，将 TTT 与 TTA 比对

AAT CTG TTT TTT
||| --- ... |||
AAT --- TTA TTT

将获得一个低得多的得分

>>> print(m["CTG", "TTA"])
7.6
>>> print(m["TTT", "TTA"])
-0.3

可能是因为 CTG 和 TTA 都编码亮氨酸，而 TTT 编码苯丙氨酸。三个字母密码子替换矩阵还揭示了代表相同氨基酸的密码子之间的偏好。例如，TTA 比 CTC 更偏好 CTG，尽管这三个密码子都编码亮氨酸

>>> s1 = ("AAT", "CTG", "CTC", "TTT")
>>> s2 = ("AAT", "TTA", "TTT")
>>> alignments = aligner.align(s1, s2)
>>> len(alignments)
1
>>> print(alignments[0])
AAT CTG CTC TTT
||| ... --- |||
AAT TTA --- TTT

>>> print(m["CTC", "TTA"])
6.5

使用匹配/不匹配分数和字母表的广义成对比对

使用三个字母的氨基酸符号，上面的序列转换为

>>> s1 = ("Asn", "Leu", "Leu", "Phe")
>>> s2 = ("Asn", "Leu", "Phe")

我们可以通过使用三个字母的氨基酸字母表将这些序列直接相互比对

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> aligner.alphabet = ['Ala', 'Arg', 'Asn', 'Asp', 'Cys',
...                     'Gln', 'Glu', 'Gly', 'His', 'Ile',
...                     'Leu', 'Lys', 'Met', 'Phe', 'Pro',
...                     'Ser', 'Thr', 'Trp', 'Tyr', 'Val']  # fmt: skip
...

我们使用 +6/-1 匹配和不匹配分数作为 BLOSUM62 矩阵的近似值，并将这些序列相互比对

>>> aligner.match = +6
>>> aligner.mismatch = -1
>>> alignments = aligner.align(s1, s2)
>>> print(len(alignments))
2
>>> print(alignments[0])
Asn Leu Leu Phe
||| ||| --- |||
Asn Leu --- Phe

>>> print(alignments[1])
Asn Leu Leu Phe
||| --- ||| |||
Asn --- Leu Phe

>>> print(alignments.score)
18.0

使用匹配/不匹配分数和整数序列的广义成对比对

在内部，执行比对的第一步是用由每个序列中每个字母在与比对器关联的字母表中的索引组成的整数数组替换两个序列。可以通过直接传递整数数组来绕过此步骤

>>> import numpy as np
>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> s1 = np.array([2, 10, 10, 13], np.int32)
>>> s2 = np.array([2, 10, 13], np.int32)
>>> aligner.match = +6
>>> aligner.mismatch = -1
>>> alignments = aligner.align(s1, s2)
>>> print(len(alignments))
2
>>> print(alignments[0])
2 10 10 13
| || -- ||
2 10 -- 13

>>> print(alignments[1])
2 10 10 13
| -- || ||
2 -- 10 13

>>> print(alignments.score)
18.0

请注意，索引应包含 32 位整数，如本示例中 numpy.int32 指定的那样。

可以通过定义通配符字符并使用相应的 Unicode 代码点号作为索引来再次包含未知字母

>>> aligner.wildcard = "?"
>>> ord(aligner.wildcard)
63
>>> s2 = np.array([2, 63, 13], np.int32)
>>> aligner.gap_score = -3
>>> alignments = aligner.align(s1, s2)
>>> print(len(alignments))
2
>>> print(alignments[0])
2 10 10 13
| .. -- ||
2 63 -- 13

>>> print(alignments[1])
2 10 10 13
| -- .. ||
2 -- 63 13

>>> print(alignments.score)
9.0

使用替换矩阵和整数序列的广义成对比对

整数序列也可以使用替换矩阵进行比对，在这种情况下，替换矩阵是一个没有关联字母表的 numpy 方阵。在这种情况下，所有索引值必须是非负的，并且小于替换矩阵的大小

>>> from Bio import Align
>>> import numpy as np
>>> aligner = Align.PairwiseAligner()
>>> m = np.eye(5)
>>> m[0, 1:] = m[1:, 0] = -2
>>> m[2, 2] = 3
>>> print(m)
[[ 1. -2. -2. -2. -2.]
 [-2.  1.  0.  0.  0.]
 [-2.  0.  3.  0.  0.]
 [-2.  0.  0.  1.  0.]
 [-2.  0.  0.  0.  1.]]
>>> aligner.substitution_matrix = m
>>> aligner.gap_score = -1
>>> s1 = np.array([0, 2, 3, 4], np.int32)
>>> s2 = np.array([0, 3, 2, 1], np.int32)
>>> alignments = aligner.align(s1, s2)
>>> print(len(alignments))
2
>>> print(alignments[0])
0 - 2 3 4
| - | . -
0 3 2 1 -

>>> print(alignments[1])
0 - 2 3 4
| - | - .
0 3 2 - 1

>>> print(alignments.score)
2.0

密码子比对

Bio.Align 模块中的 CodonAligner 类实现了一个专门的比对器，用于将核苷酸序列比对到它编码的氨基酸序列。如果在翻译过程中发生移码，此类比对将非常重要。

将核苷酸序列比对到氨基酸序列

要将核苷酸序列比对到氨基酸序列，请首先创建一个 CodonAligner 对象

>>> from Bio import Align
>>> aligner = Align.CodonAligner()

CodonAligner 对象 aligner 存储将用于比对的对齐参数

>>> print(aligner)
Codon aligner with parameters
  wildcard: 'X'
  match_score: 1.0
  mismatch_score: 0.0
  frameshift_minus_two_score: -3.0
  frameshift_minus_one_score: -3.0
  frameshift_plus_one_score: -3.0
  frameshift_plus_two_score: -3.0

参数 wildcard、match_score 和 mismatch_score 的定义方式与上面描述的 PairwiseAligner 类相同（请参见第成对比对器对象节）。参数 frameshift_minus_two_score、frameshift_minus_one_score、frameshift_plus_one_score 和 frameshift_plus_two_score 指定的值在比对中发生 -2、-1、+1 或 +2 移码时会添加到比对得分中。默认情况下，移码得分设置为 -3.0。与 PairwiseAligner 类（表成对比对器对象的元属性。）类似，CodonAligner 类定义了其他属性，这些属性共同引用了许多这些值，如表 CodonAligner 对象的元属性。所示。

表 2 CodonAligner 对象的元属性。
元属性	它映射到的属性
`frameshift_minus_score`	`frameshift_minus_two_score`、`frameshift_minus_one_score`
`frameshift_plus_score`	`frameshift_plus_two_score`、`frameshift_plus_one_score`
`frameshift_two_score`	`frameshift_minus_two_score`、`frameshift_plus_two_score`
`frameshift_one_score`	`frameshift_minus_one_score`、`frameshift_plus_one_score`
`frameshift_score`	`frameshift_minus_two_score`、`frameshift_minus_one_score`、`frameshift_plus_one_score`、`frameshift_plus_two_score`

现在让我们考虑两个核苷酸序列以及它们编码的氨基酸序列

>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> nuc1 = Seq("TCAGGGACTGCGAGAACCAAGCTACTGCTGCTGCTGGCTGCGCTCTGCGCCGCAGGTGGGGCGCTGGAG")
>>> rna1 = SeqRecord(nuc1, id="rna1")
>>> nuc2 = Seq("TCAGGGACTTCGAGAACCAAGCGCTCCTGCTGCTGGCTGCGCTCGGCGCCGCAGGTGGAGCACTGGAG")
>>> rna2 = SeqRecord(nuc2, id="rna2")
>>> aa1 = Seq("SGTARTKLLLLLAALCAAGGALE")
>>> aa2 = Seq("SGTSRTKRLLLLAALGAAGGALE")
>>> pro1 = SeqRecord(aa1, id="pro1")
>>> pro2 = SeqRecord(aa2, id="pro2")

虽然两个蛋白质序列都包含 23 个氨基酸，但第一个核苷酸序列包含 \(3 \times 23 = 69\) 个核苷酸，而第二个核苷酸序列只包含 68 个核苷酸

>>> len(pro1)
23
>>> len(pro2)
23
>>> len(rna1)
69
>>> len(rna2)
68

这是由于在翻译第二个核苷酸序列期间发生了 -1 移码事件。使用 CodonAligner.align 将 rna1 比对到 pro1，并将 rna2 比对到 pro2，返回 Alignment 对象的迭代器

>>> alignments1 = aligner.align(pro1, rna1)
>>> len(alignments1)
1
>>> alignment1 = next(alignments1)
>>> print(alignment1)
pro1              0 S  G  T  A  R  T  K  L  L  L  L  L  A  A  L  C  A  A  G  G
rna1              0 TCAGGGACTGCGAGAACCAAGCTACTGCTGCTGCTGGCTGCGCTCTGCGCCGCAGGTGGG

pro1             20 A  L  E   23
rna1             60 GCGCTGGAG 69

>>> alignment1.coordinates
array([[ 0, 23],
       [ 0, 69]])
>>> alignment1[0]
'SGTARTKLLLLLAALCAAGGALE'
>>> alignment1[1]
'TCAGGGACTGCGAGAACCAAGCTACTGCTGCTGCTGGCTGCGCTCTGCGCCGCAGGTGGGGCGCTGGAG'
>>> alignments2 = aligner.align(pro2, rna2)
>>> len(alignments2)
1
>>> alignment2 = next(alignments2)
>>> print(alignment2)
pro2              0 S  G  T  S  R  T  K  R   8
rna2              0 TCAGGGACTTCGAGAACCAAGCGC 24

pro2              8 L  L  L  L  A  A  L  G  A  A  G  G  A  L  E   23
rna2             23 CTCCTGCTGCTGGCTGCGCTCGGCGCCGCAGGTGGAGCACTGGAG 68

>>> alignment2[0]
'SGTSRTKRLLLLAALGAAGGALE'
>>> alignment2[1]
'TCAGGGACTTCGAGAACCAAGCGCCTCCTGCTGCTGGCTGCGCTCGGCGCCGCAGGTGGAGCACTGGAG'
>>> alignment2.coordinates
array([[ 0,  8,  8, 23],
       [ 0, 24, 23, 68]])

虽然 alignment1 是 69 个核苷酸与 23 个氨基酸的连续比对，但在 alignment2 中，我们发现在 24 个核苷酸之后发生了 -1 移码。由于 alignment2[1] 包含应用 -1 移码后的核苷酸序列，因此它比 nuc2 长一个核苷酸，可以直接翻译，从而产生氨基酸序列 aa2

>>> from Bio.Seq import translate
>>> len(nuc2)
68
>>> len(alignment2[1])
69
>>> translate(alignment2[1])
'SGTSRTKRLLLLAALGAAGGALE'
>>> _ == aa2
True

对齐得分存储为 alignments1 和 alignments2 迭代器以及单个比对 alignment1 和 alignment2 上的属性

>>> alignments1.score
23.0
>>> alignment1.score
23.0
>>> alignments2.score
20.0
>>> alignment2.score
20.0

其中，rna1-pro1 比对的得分等于比对的氨基酸数量，而 rna2-pro2 比对的得分由于移码惩罚而低 3 分。要计算比对得分而不实际进行比对，可以使用 score 方法。

>>> score = aligner.score(pro1, rna1)
>>> print(score)
23.0
>>> score = aligner.score(pro2, rna2)
>>> print(score)
20.0

生成密码子序列的多序列比对

假设我们还有第三个相关的氨基酸序列及其对应的核苷酸序列

>>> aa3 = Seq("MGTALLLLLAALCAAGGALE")
>>> pro3 = SeqRecord(aa3, id="pro3")
>>> nuc3 = Seq("ATGGGAACCGCGCTGCTTTTGCTACTGGCCGCGCTCTGCGCCGCAGGTGGGGCCCTGGAG")
>>> rna3 = SeqRecord(nuc3, id="rna3")
>>> nuc3.translate() == aa3
True

如上，我们使用 CodonAligner 将核苷酸序列比对到氨基酸序列

>>> alignments3 = aligner.align(pro3, rna3)
>>> len(alignments3)
1
>>> alignment3 = next(alignments3)
>>> print(alignment3)
pro3              0 M  G  T  A  L  L  L  L  L  A  A  L  C  A  A  G  G  A  L  E
rna3              0 ATGGGAACCGCGCTGCTTTTGCTACTGGCCGCGCTCTGCGCCGCAGGTGGGGCCCTGGAG

pro3             20
rna3             60

可以使用 ClustalW 等工具将三个氨基酸序列相互比对。在这里，我们手动创建比对结果

>>> import numpy as np
>>> from Bio.Align import Alignment
>>> sequences = [pro1, pro2, pro3]
>>> protein_alignment = Alignment(
...     sequences, coordinates=np.array([[0, 4, 7, 23], [0, 4, 7, 23], [0, 4, 4, 20]])
... )
>>> print(protein_alignment)
pro1              0 SGTARTKLLLLLAALCAAGGALE 23
pro2              0 SGTSRTKRLLLLAALGAAGGALE 23
pro3              0 MGTA---LLLLLAALCAAGGALE 20

现在，我们可以使用 mapall 方法对蛋白质比对进行操作，并将核苷酸-蛋白质成对比对作为参数，以获得相应的密码子比对结果

>>> codon_alignment = protein_alignment.mapall([alignment1, alignment2, alignment3])
>>> print(codon_alignment)
rna1              0 TCAGGGACTGCGAGAACCAAGCTA 24
rna2              0 TCAGGGACTTCGAGAACCAAGCGC 24
rna3              0 ATGGGAACCGCG---------CTG 15

rna1             24 CTGCTGCTGCTGGCTGCGCTCTGCGCCGCAGGTGGGGCGCTGGAG 69
rna2             23 CTCCTGCTGCTGGCTGCGCTCGGCGCCGCAGGTGGAGCACTGGAG 68
rna3             15 CTTTTGCTACTGGCCGCGCTCTGCGCCGCAGGTGGGGCCCTGGAG 60

分析密码子比对

计算每个位点的非同义替换数和同义替换数

密码子比对最重要的应用是估计每个位点的非同义替换数 (dN) 和同义替换数 (dS)。这些可以通过 Bio.Align.analysis 中的 calculate_dn_ds 函数计算。此函数采用成对密码子比对和输入作为参数，以及可选参数 method 指定计算方法，codon_table（默认值为标准密码子表），过渡和颠换速率之比 k，以及 cfreq 指定平衡密码子频率。 Biopython 目前支持三种基于计数的方法（NG86、LWL85、YN00）以及最大似然法 (ML) 来估计 dN 和 dS

NG86：Nei 和 Gojobori (1986) [Nei1986]（默认）。使用此方法，您还可以通过参数 k 指定过渡和颠换速率之比，默认值为 1.0。
LWL85：Li 等人 (1985) [Li1985].
YN00：Yang 和 Nielsen (2000) [Yang2000].
ML：Goldman 和 Yang (1994) [Goldman1994]。使用此方法，您还可以通过 cfreq 参数指定平衡密码子频率，选项如下：
- F1x4：统计提供的密码子序列中的核苷酸频率，并使用它计算背景密码子频率；
- F3x4：（默认）分别统计提供的密码子中第一个、第二个和第三个位置的核苷酸频率，并使用它计算背景密码子频率；
- F61：统计提供的密码子序列中密码子的频率，并使用 0.1 的伪计数。

calculate_dN_dS 方法可以应用于成对密码子比对。通常，不同的计算方法会导致 dN 和 dS 的估计值略有不同

>>> from Bio.Align import analysis
>>> pairwise_codon_alignment = codon_alignment[:2]
>>> print(pairwise_codon_alignment)
rna1              0 TCAGGGACTGCGAGAACCAAGCTA 24
                  0 |||||||||.||||||||||||..
rna2              0 TCAGGGACTTCGAGAACCAAGCGC 24

rna1             24 CTGCTGCTGCTGGCTGCGCTCTGCGCCGCAGGTGGGGCGCTGGAG 69
                 24 ||.||||||||||||||||||.|||||||||||||.||.|||||| 69
rna2             23 CTCCTGCTGCTGGCTGCGCTCGGCGCCGCAGGTGGAGCACTGGAG 68

>>> dN, dS = analysis.calculate_dn_ds(pairwise_codon_alignment, method="NG86")
>>> print(dN, dS)  
0.067715... 0.201197...
>>> dN, dS = analysis.calculate_dn_ds(pairwise_codon_alignment, method="LWL85")
>>> print(dN, dS)  
0.068728... 0.207551...

>>> dN, dS = analysis.calculate_dn_ds(pairwise_codon_alignment, method="YN00")
>>> print(dN, dS)  
0.081468... 0.127706...
>>> dN, dS = analysis.calculate_dn_ds(pairwise_codon_alignment, method="ML")
>>> print(dN, dS)  
0.069475... 0.205754...

对于密码子序列的多重比对，您可以计算 dN 和 dS 值的矩阵

>>> dN, dS = analysis.calculate_dn_ds_matrix(codon_alignment, method="NG86")
>>> print(dN)
rna1    0.000000
rna2    0.067715    0.000000
rna3    0.060204    0.145469    0.000000
    rna1    rna2    rna3
>>> print(dS)
rna1    0.000000
rna2    0.201198    0.000000
rna3    0.664268    0.798957    0.000000
    rna1    rna2    rna3

calculate_dn_ds_matrix 返回的 dN 和 dS 对象是 Bio.Phylo.TreeConstruction 中 DistanceMatrix 类的实例。此函数只将 codon_table 作为可选参数。

从这两个序列，您可以使用 Bio.Phylo.TreeConstruction 创建 dN 树和 dS 树

>>> from Bio.Phylo.TreeConstruction import DistanceTreeConstructor
>>> dn_constructor = DistanceTreeConstructor()
>>> ds_constructor = DistanceTreeConstructor()
>>> dn_tree = dn_constructor.upgma(dN)
>>> ds_tree = ds_constructor.upgma(dS)
>>> print(type(dn_tree))
<class 'Bio.Phylo.BaseTree.Tree'>
>>> print(dn_tree)  
Tree(rooted=True)
    Clade(branch_length=0, name='Inner2')
        Clade(branch_length=0.053296..., name='rna2')
        Clade(branch_length=0.023194..., name='Inner1')
            Clade(branch_length=0.0301021..., name='rna3')
            Clade(branch_length=0.0301021..., name='rna1')
>>> print(ds_tree)  
Tree(rooted=True)
    Clade(branch_length=0, name='Inner2')
        Clade(branch_length=0.365806..., name='rna3')
        Clade(branch_length=0.265207..., name='Inner1')
            Clade(branch_length=0.100598..., name='rna2')
            Clade(branch_length=0.100598..., name='rna1')

执行 McDonald-Kreitman 检验

McDonald-Kreitman 检验通过比较种内同义替换和非同义替换与种间同义替换和非同义替换，来评估适应性进化的程度，以查看它们是否来自相同的进化过程。该检验需要从同一物种的不同个体中取样的基因序列。在以下示例中，我们将使用果蝇的 Adh 基因。数据包括来自黑腹果蝇的 11 个个体、来自拟黑腹果蝇的 4 个个体以及来自雅库布果蝇的 12 个个体。蛋白质比对数据和核苷酸序列可在 Biopython 发行版中的 Tests/codonalign 目录中找到，分别为文件 adh.aln 和 drosophila.fasta。 Bio.Align.analysis 中的 mktest 函数实现了 Mcdonald-Kreitman 检验。

>>> from Bio import SeqIO
>>> from Bio import Align
>>> from Bio.Align import CodonAligner
>>> from Bio.Align.analysis import mktest
>>> aligner = CodonAligner()
>>> nucleotide_records = SeqIO.index("drosophila.fasta", "fasta")
>>> for nucleotide_record in nucleotide_records.values():
...     print(nucleotide_record.description)  
...
gi|9097|emb|X57361.1| Drosophila simulans (individual c) ...
gi|9099|emb|X57362.1| Drosophila simulans (individual d) ...
gi|9101|emb|X57363.1| Drosophila simulans (individual e) ...
gi|9103|emb|X57364.1| Drosophila simulans (individual f) ...
gi|9217|emb|X57365.1| Drosophila yakuba (individual a) ...
gi|9219|emb|X57366.1| Drosophila yakuba (individual b) ...
gi|9221|emb|X57367.1| Drosophila yakuba (individual c) ...
gi|9223|emb|X57368.1| Drosophila yakuba (individual d) ...
gi|9225|emb|X57369.1| Drosophila yakuba (individual e) ...
gi|9227|emb|X57370.1| Drosophila yakuba (individual f) ...
gi|9229|emb|X57371.1| Drosophila yakuba (individual g) ...
gi|9231|emb|X57372.1| Drosophila yakuba (individual h) ...
gi|9233|emb|X57373.1| Drosophila yakuba (individual i) ...
gi|9235|emb|X57374.1| Drosophila yakuba (individual j) ...
gi|9237|emb|X57375.1| Drosophila yakuba (individual k) ...
gi|9239|emb|X57376.1| Drosophila yakuba (individual l) ...
gi|156879|gb|M17837.1|DROADHCK D.melanogaster (strain Ja-F) ...
gi|156863|gb|M19547.1|DROADHCC D.melanogaster (strain Af-S) ...
gi|156877|gb|M17836.1|DROADHCJ D.melanogaster (strain Af-F) ...
gi|156875|gb|M17835.1|DROADHCI D.melanogaster (strain Wa-F) ...
gi|156873|gb|M17834.1|DROADHCH D.melanogaster (strain Fr-F) ...
gi|156871|gb|M17833.1|DROADHCG D.melanogaster (strain Fl-F) ...
gi|156869|gb|M17832.1|DROADHCF D.melanogaster (strain Ja-S) ...
gi|156867|gb|M17831.1|DROADHCE D.melanogaster (strain Fl-2S) ...
gi|156865|gb|M17830.1|DROADHCD D.melanogaster (strain Fr-S) ...
gi|156861|gb|M17828.1|DROADHCB D.melanogaster (strain Fl-1S) ...
gi|156859|gb|M17827.1|DROADHCA D.melanogaster (strain Wa-S) ...
>>> protein_alignment = Align.read("adh.aln", "clustal")
>>> len(protein_alignment)
27
>>> print(protein_alignment)  
gi|9217|e         0 MAFTLTNKNVVFVAGLGGIGLDTSKELVKRDLKNLVILDRIENPAAIAELKAINPKVTVT
gi|9219|e         0 MAFTLTNKNVVFVAGLGGIGLDTSKELVKRDLKNLVILDRIENPAAIAELKAINPKVTVT
gi|9221|e         0 MAFTLTNKNVVFVAGLGGIGLDTSKELVKRDLKNLVILDRIENPAAIAELKAINPKVTVT
...
gi|156859         0 MSFTLTNKNVIFVAGLGGIGLDTSKELLKRDLKNLVILDRIENPAAIAELKAINPKVTVT

...

gi|9217|e       240 GTLEAIQWSKHWDSGI 256
gi|9219|e       240 GTLEAIQWSKHWDSGI 256
gi|9221|e       240 GTLEAIQWSKHWDSGI 256
...
gi|156859       240 GTLEAIQWTKHWDSGI 256

>>> codon_alignments = []
>>> for protein_record in protein_alignment.sequences:
...     nucleotide_record = nucleotide_records[protein_record.id]
...     alignments = aligner.align(protein_record, nucleotide_record)
...     assert len(alignments) == 1
...     codon_alignment = next(alignments)
...     codon_alignments.append(codon_alignment)
...
>>> print(codon_alignment)  
gi|156859         0 M  S  F  T  L  T  N  K  N  V  I  F  V  A  G  L  G  G  I  G
gi|156859         0 ATGTCGTTTACTTTGACCAACAAGAACGTGATTTTCGTTGCCGGTCTGGGAGGCATTGGT

gi|156859        20 L  D  T  S  K  E  L  L  K  R  D  L  K  N  L  V  I  L  D  R
gi|156859        60 CTGGACACCAGCAAGGAGCTGCTCAAGCGCGATCTGAAGAACCTGGTGATCCTCGACCGC

...

gi|156859       240 G  T  L  E  A  I  Q  W  T  K  H  W  D  S  G  I   256
gi|156859       720 GGCACCCTGGAGGCCATCCAGTGGACCAAGCACTGGGACTCCGGCATC 768

>>> nucleotide_records.close()  # Close indexed FASTA file
>>> alignment = protein_alignment.mapall(codon_alignments)
>>> print(alignment)  
gi|9217|e         0 ATGGCGTTTACCTTGACCAACAAGAACGTGGTTTTCGTGGCCGGTCTGGGAGGCATTGGT
gi|9219|e         0 ATGGCGTTTACCTTGACCAACAAGAACGTGGTTTTCGTGGCCGGTCTGGGAGGCATTGGT
gi|9221|e         0 ATGGCGTTTACCTTGACCAACAAGAACGTGGTTTTCGTGGCCGGTCTGGGAGGCATTGGT
...
gi|156859         0 ATGTCGTTTACTTTGACCAACAAGAACGTGATTTTCGTTGCCGGTCTGGGAGGCATTGGT

...

gi|9217|e       720 GGCACCCTGGAGGCCATCCAGTGGTCCAAGCACTGGGACTCCGGCATC 768
gi|9219|e       720 GGCACCCTGGAGGCCATCCAGTGGTCCAAGCACTGGGACTCCGGCATC 768
gi|9221|e       720 GGTACCCTGGAGGCCATCCAGTGGTCCAAGCACTGGGACTCCGGCATC 768
...
gi|156859       720 GGCACCCTGGAGGCCATCCAGTGGACCAAGCACTGGGACTCCGGCATC 768

>>> unique_species = ["Drosophila simulans", "Drosophila yakuba", "D.melanogaster"]
>>> species = []
>>> for record in alignment.sequences:
...     description = record.description
...     for s in unique_species:
...         if s in description:
...             break
...     else:
...         raise Exception(f"Failed to find species for {description}")
...     species.append(s)
...
>>> print(species)
['Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila yakuba', 'Drosophila simulans', 'Drosophila simulans', 'Drosophila simulans', 'Drosophila simulans', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster', 'D.melanogaster']
>>> pvalue = mktest(alignment, species)
>>> print(pvalue)  
0.00206457...

除了多重密码子比对之外，mktest 函数还接收比对中每个序列所属的物种作为输入。可以将密码子表作为可选参数 codon_table 提供。