成對序列比對

成對序列比對是將兩個序列彼此比對，並最佳化它們之間的相似性分數的過程。Bio.Align 模組包含 PairwiseAligner 類別，用於使用 Needleman-Wunsch、Smith-Waterman、Gotoh（三態）和 Waterman-Smith-Beyer 全域和局部成對比對演算法進行全域和局部比對，並提供許多選項來變更比對參數。關於序列比對演算法的深入資訊，我們參考 Durbin et al. [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]]])

在成對比對 alignment1 上可以使用 map 方法，以尋找 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] 對於使用 megablast 進行的核苷酸比對，會使用 \(+1\) 的匹配分數和 \(-2\) 的不匹配分數，間隙罰分為 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] 對於使用 blastp 進行的蛋白質比對，會使用 BLOSUM62 取代矩陣。此取代矩陣可從 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 物件的以下十二個屬性，來精細控制間隙計分方案

開啟分數

延伸分數

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

這些屬性允許序列內部和任一端有不同的間隙分數，如本範例所示

目標

查詢

分數

A

查詢左側開啟間隙分數

C

查詢左側延伸間隙分數

C

查詢左側延伸間隙分數

G

G

匹配分數

G

T

不匹配分數

G

查詢內部開啟間隙分數

A

查詢內部延伸間隙分數

A

查詢內部延伸間隙分數

T

T

匹配分數

A

A

匹配分數

G

查詢內部開啟間隙分數

C

C

匹配分數

C

目標內部開啟間隙分數

C

目標內部延伸間隙分數

C

C

匹配分數

T

G

不匹配分數

C

C

匹配分數

C

目標內部開啟間隙分數

A

A

匹配分數

T

目標右側開啟間隙分數

A

目標右側延伸間隙分數

A

目標右側延伸間隙分數

為了方便起見，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）。雖然它們看起來類似於 Alignment 物件的 tuple 或 list，但它們的不同之處在於每個 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 陣列的子類別，因此可以這樣使用。如果出現在數學運算中的 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 陣列的子類別，因此您可以使用非常相似的方式對 Array 物件應用數學運算。在此，我們透過計算來自大腸桿菌和枯草桿菌的 16S 核醣體 RNA 基因序列比對的評分矩陣來說明這一點。首先，我們建立一個 PairwiseAligner 物件（請參閱第 成對序列比對 節），並使用 blastn 使用的預設分數初始化它

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

接下來，我們讀取大腸桿菌和枯草桿菌的 16S 核醣體 RNA 基因序列（在 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 種胺基酸、3 種含糊胺基酸（B：天冬醯胺或天冬胺酸，Z：麩醯胺或麩胺酸，以及 X：代表任何胺基酸）和以星號表示的終止密碼子組成。

>>> m.alphabet
'ARNDCQEGHILKMFPSTWYVBZX*'

要取得可用替換矩陣的完整列表，請使用不帶引數的 load：

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

請注意，Schneider *et al.* [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 *et al.* [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 對於 CTG 的偏好高於 CTC，儘管這三個密碼子都編碼亮胺酸。

>>> 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 類別相同（請參閱第 成對比對物件 節）。每當比對中出現 -2、-1、+1 或 +2 的讀碼框位移時，frameshift_minus_two_score、frameshift_minus_one_score、frameshift_plus_one_score 和 frameshift_plus_two_score 參數指定的值會加到比對分數中。預設情況下，讀碼框位移分數設定為 -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 et al. (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')

執行麥克唐納-克雷特曼檢定

麥克唐納-克雷特曼檢定透過比較物種內的同義取代和非同義取代與物種間的同義取代和非同義取代，來評估適應性演化的程度，以查看它們是否來自相同的演化過程。該檢定需要從同一物種的不同個體中取樣的基因序列。在以下範例中，我們將使用果蠅的 Adh 基因。這些資料包括來自黑腹果蠅的 11 個個體、來自擬果蠅的 4 個個體和來自屋久果蠅的 12 個個體。蛋白質比對資料和核苷酸序列可在 Biopython 發行版的 Tests/codonalign 目錄中找到，檔案分別為 adh.aln 和 drosophila.fasta。Bio.Align.analysis 中的函數 mktest 實作了麥克唐納-克雷特曼檢定。

>>> 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 提供。