多序列比對物件

本章介紹較舊的 MultipleSeqAlignment 類別，以及 Bio.AlignIO 中的解析器，這些解析器會解析序列比對軟體的輸出，產生 MultipleSeqAlignment 物件。多序列比對指的是已對齊在一起的多個序列集合 — 通常是插入間隙字元，並加入前導或後續間隙 — 使得所有序列字串的長度都相同。這種比對可以視為字母矩陣，其中每一列都在內部以 SeqRecord 物件保存。

我們將介紹保存此類資料的 MultipleSeqAlignment 物件，以及用於以各種檔案格式讀寫它們的 Bio.AlignIO 模組（遵循上一章 Bio.SeqIO 模組的設計）。請注意，Bio.SeqIO 和 Bio.AlignIO 都可以讀寫序列比對檔案。適當的選擇在很大程度上取決於您想要如何處理資料。

本章的最後一部分是有關於從 Python 使用常見的多序列比對工具（如 ClustalW 和 MUSCLE），並使用 Biopython 解析結果。

解析或讀取序列比對

我們有兩個函數用於讀取序列比對，Bio.AlignIO.read() 和 Bio.AlignIO.parse()，它們遵循 Bio.SeqIO 中介紹的慣例，分別用於包含一個或多個比對的檔案。

使用 Bio.AlignIO.parse() 將會返回一個 迭代器，它會給出 MultipleSeqAlignment 物件。迭代器通常在 for 迴圈中使用。您會有許多不同比對的情況範例，包括來自 PHYLIP 工具 seqboot 的重新採樣比對，或來自 EMBOSS 工具 water 或 needle 的多個成對比對，或 Bill Pearson 的 FASTA 工具。

但是，在許多情況下，您將處理僅包含單一比對的檔案。在這種情況下，您應該使用 Bio.AlignIO.read() 函數，它會返回一個單一 MultipleSeqAlignment 物件。

這兩個函數都需要兩個強制參數

1. 第一個參數是一個 句柄，用於讀取資料，通常是一個開啟的檔案（請參閱第 句柄到底是什麼？ 節），或是一個檔案名稱。

2. 第二個參數是一個小寫字串，用於指定比對格式。如同在 Bio.SeqIO 中一樣，我們不會嘗試為您猜測檔案格式！請參閱 https://biopython.dev.org.tw/wiki/AlignIO 以取得完整支援格式列表。

還有一個可選的 seq_count 參數，在下面第 模糊比對 節中討論，用於處理可能包含多個比對的模糊檔案格式。

單一比對

作為範例，請考慮以下 PFAM 或斯德哥爾摩檔案格式中註解豐富的蛋白質比對

# STOCKHOLM 1.0
#=GS COATB_BPIKE/30-81  AC P03620.1
#=GS COATB_BPIKE/30-81  DR PDB; 1ifl ; 1-52;
#=GS Q9T0Q8_BPIKE/1-52  AC Q9T0Q8.1
#=GS COATB_BPI22/32-83  AC P15416.1
#=GS COATB_BPM13/24-72  AC P69541.1
#=GS COATB_BPM13/24-72  DR PDB; 2cpb ; 1-49;
#=GS COATB_BPM13/24-72  DR PDB; 2cps ; 1-49;
#=GS COATB_BPZJ2/1-49   AC P03618.1
#=GS Q9T0Q9_BPFD/1-49   AC Q9T0Q9.1
#=GS Q9T0Q9_BPFD/1-49   DR PDB; 1nh4 A; 1-49;
#=GS COATB_BPIF1/22-73  AC P03619.2
#=GS COATB_BPIF1/22-73  DR PDB; 1ifk ; 1-50;
COATB_BPIKE/30-81             AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA
#=GR COATB_BPIKE/30-81  SS    -HHHHHHHHHHHHHH--HHHHHHHH--HHHHHHHHHHHHHHHHHHHHH----
Q9T0Q8_BPIKE/1-52             AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA
COATB_BPI22/32-83             DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA
COATB_BPM13/24-72             AEGDDP...AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
#=GR COATB_BPM13/24-72  SS    ---S-T...CHCHHHHCCCCTCCCTTCHHHHHHHHHHHHHHHHHHHHCTT--
COATB_BPZJ2/1-49              AEGDDP...AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA
Q9T0Q9_BPFD/1-49              AEGDDP...AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
#=GR Q9T0Q9_BPFD/1-49   SS    ------...-HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH--
COATB_BPIF1/22-73             FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA
#=GR COATB_BPIF1/22-73  SS    XX-HHHH--HHHHHH--HHHHHHH--HHHHHHHHHHHHHHHHHHHHHHH---
#=GC SS_cons                  XHHHHHHHHHHHHHHHCHHHHHHHHCHHHHHHHHHHHHHHHHHHHHHHHC--
#=GC seq_cons                 AEssss...AptAhDSLpspAT-hIu.sWshVsslVsAsluIKLFKKFsSKA
//

這是 Phage_Coat_Gp8 (PF05371) PFAM 條目的種子比對，從 PFAM 的過時版本（來自 https://pfam.xfam.org/）下載。我們可以如下載入此檔案（假設它已在目前的工作目錄中儲存為「PF05371_seed.sth」）

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")

此程式碼將會列印比對的摘要

>>> print(alignment)
Alignment with 7 rows and 52 columns
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRL...SKA COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKL...SRA Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRL...SKA COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73

您會在上面的輸出中注意到序列已被截斷。我們可以改為編寫自己的程式碼，將其格式化為我們喜歡的樣子，方法是將各列以 SeqRecord 物件進行迭代

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> print("Alignment length %i" % alignment.get_alignment_length())
Alignment length 52
>>> for record in alignment:
...     print("%s - %s" % (record.seq, record.id))
...
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA - COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA - Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA - COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA - COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA - COATB_BPIF1/22-73

您也可以對比對物件呼叫 Python 的內建 format 函數，以特定的檔案格式顯示它 – 詳細資訊請參閱第 將你的比對物件轉換為格式化的字串 節。

您是否注意到上面原始檔案中的一些序列包括對 PDB 的資料庫交叉參考，以及相關的已知二級結構？試試這個

>>> for record in alignment:
...     if record.dbxrefs:
...         print("%s %s" % (record.id, record.dbxrefs))
...
COATB_BPIKE/30-81 ['PDB; 1ifl ; 1-52;']
COATB_BPM13/24-72 ['PDB; 2cpb ; 1-49;', 'PDB; 2cps ; 1-49;']
Q9T0Q9_BPFD/1-49 ['PDB; 1nh4 A; 1-49;']
COATB_BPIF1/22-73 ['PDB; 1ifk ; 1-50;']

若要查看所有序列註解，請嘗試這個

>>> for record in alignment:
...     print(record)
...

PFAM 在 http://pfam.xfam.org/family/PF05371 提供了一個不錯的網頁介面，它實際上可讓您以其他幾種格式下載此比對。以下是該檔案在 FASTA 檔案格式中的樣子

>COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA
>Q9T0Q8_BPIKE/1-52
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA
>COATB_BPI22/32-83
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA
>COATB_BPM13/24-72
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
>COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA
>Q9T0Q9_BPFD/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA
>COATB_BPIF1/22-73
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA

請注意，該網站應該有一個選項，說明如何將間隙顯示為句號（點）或破折號，我們在上面顯示了破折號。假設您下載並將其儲存為檔案「PF05371_seed.faa」，則可以使用幾乎完全相同的程式碼來載入它

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.faa", "fasta")
>>> print(alignment)

此程式碼中唯一變更的是檔案名稱和格式字串。您會獲得與之前相同的輸出，序列和記錄識別碼都相同。但是，正如您所預期的那樣，如果您檢查每個 SeqRecord，則沒有註解或資料庫交叉參考，因為這些未包含在 FASTA 檔案格式中。

請注意，您可以使用 Bio.AlignIO 將原始的斯德哥爾摩格式檔案轉換為您自己的 FASTA 檔案，而不是使用 Sanger 網站（請參閱下文）。

使用任何受支援的檔案格式，您都可以透過變更格式字串，以完全相同的方式載入比對。例如，PHYLIP 檔案使用「phylip」，NEXUS 檔案使用「nexus」，EMBOSS 工具輸出的比對則使用「emboss」。維基頁面上（https://biopython.dev.org.tw/wiki/AlignIO）和內建的文件中都有完整的列表，Bio.AlignIO

>>> from Bio import AlignIO
>>> help(AlignIO)

多重比對

上一節重點介紹讀取包含單一比對的檔案。然而，一般來說，檔案可以包含多個比對，若要讀取這些檔案，我們必須使用 Bio.AlignIO.parse() 函數。

假設您有一個小型的 PHYLIP 格式比對

    5    6
Alpha     AACAAC
Beta      AACCCC
Gamma     ACCAAC
Delta     CCACCA
Epsilon   CCAAAC

如果您想使用 PHYLIP 工具進行自助法系統發生樹，其中一個步驟是使用工具 bootseq 建立一組重新取樣的比對。這將會產生類似這樣的輸出，為了簡潔起見，已經縮短了

    5     6
Alpha     AAACCA
Beta      AAACCC
Gamma     ACCCCA
Delta     CCCAAC
Epsilon   CCCAAA
    5     6
Alpha     AAACAA
Beta      AAACCC
Gamma     ACCCAA
Delta     CCCACC
Epsilon   CCCAAA
    5     6
Alpha     AAAAAC
Beta      AAACCC
Gamma     AACAAC
Delta     CCCCCA
Epsilon   CCCAAC
...
    5     6
Alpha     AAAACC
Beta      ACCCCC
Gamma     AAAACC
Delta     CCCCAA
Epsilon   CAAACC

如果您想使用 Bio.AlignIO 讀取它，則可以使用

>>> from Bio import AlignIO
>>> alignments = AlignIO.parse("resampled.phy", "phylip")
>>> for alignment in alignments:
...     print(alignment)
...     print()
...

這將會產生以下輸出，為了顯示，再次進行了縮短

Alignment with 5 rows and 6 columns
AAACCA Alpha
AAACCC Beta
ACCCCA Gamma
CCCAAC Delta
CCCAAA Epsilon

Alignment with 5 rows and 6 columns
AAACAA Alpha
AAACCC Beta
ACCCAA Gamma
CCCACC Delta
CCCAAA Epsilon

Alignment with 5 rows and 6 columns
AAAAAC Alpha
AAACCC Beta
AACAAC Gamma
CCCCCA Delta
CCCAAC Epsilon

...

Alignment with 5 rows and 6 columns
AAAACC Alpha
ACCCCC Beta
AAAACC Gamma
CCCCAA Delta
CAAACC Epsilon

如同函數 Bio.SeqIO.parse()，使用 Bio.AlignIO.parse() 會返回一個迭代器。如果您想一次將所有比對都保存在記憶體中，這樣您就可以按任何順序存取它們，則將迭代器轉換為列表

>>> from Bio import AlignIO
>>> alignments = list(AlignIO.parse("resampled.phy", "phylip"))
>>> last_align = alignments[-1]
>>> first_align = alignments[0]

模糊比對

許多比對檔案格式可以明確地儲存多個比對，並且每個比對之間的區隔都很清楚。但是，當使用一般序列檔案格式時，沒有這種區塊結構。最常見的情況是比對已儲存為 FASTA 檔案格式。例如，請考慮以下情況

>Alpha
ACTACGACTAGCTCAG--G
>Beta
ACTACCGCTAGCTCAGAAG
>Gamma
ACTACGGCTAGCACAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Beta
ACTACCGCTAGCTCAGAAG
>Gamma
ACTACGGCTAGCACAGAAG

這可能是一個包含六個序列的單一比對（具有重複的識別碼）。或者，從識別碼來看，這可能是兩個不同的比對，每個比對都有三個序列，它們的長度恰好都相同。

下一個範例呢？

>Alpha
ACTACGACTAGCTCAG--G
>Beta
ACTACCGCTAGCTCAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Gamma
ACTACGGCTAGCACAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Delta
ACTACGGCTAGCACAGAAG

同樣地，這可能是一個包含六個序列的單一比對。但是這次根據識別碼，我們可以猜測這是三個成對的比對，它們碰巧都具有相同的長度。

最後一個範例類似

>Alpha
ACTACGACTAGCTCAG--G
>XXX
ACTACCGCTAGCTCAGAAG
>Alpha
ACTACGACTAGCTCAGG
>YYY
ACTACGGCAAGCACAGG
>Alpha
--ACTACGAC--TAGCTCAGG
>ZZZ
GGACTACGACAATAGCTCAGG

在第三個範例中，由於長度不同，因此不能將其視為包含所有六個記錄的單一比對。但是，它可以是三個成對的比對。

顯然，嘗試在 FASTA 檔案中儲存多個比對並非理想。但是，如果您被迫將這些作為輸入檔案處理，則 Bio.AlignIO 可以應付所有比對都具有相同記錄數量的最常見情況。其中一個範例是成對比對的集合，它們可以由 EMBOSS 工具 needle 和 water 產生 — 儘管在這種情況下，Bio.AlignIO 應該能夠使用「emboss」作為格式字串來理解它們的原生輸出。

若要將這些 FASTA 範例解讀為幾個獨立的比對結果，我們可以搭配選用的 seq_count 參數來使用 Bio.AlignIO.parse()，此參數指定每個比對結果中預期的序列數量（在此範例中，分別為 3、2 和 2）。例如，以第三個範例作為輸入資料：

>>> for alignment in AlignIO.parse(handle, "fasta", seq_count=2):
...     print("Alignment length %i" % alignment.get_alignment_length())
...     for record in alignment:
...         print("%s - %s" % (record.seq, record.id))
...     print()
...

會得到：

Alignment length 19
ACTACGACTAGCTCAG--G - Alpha
ACTACCGCTAGCTCAGAAG - XXX

Alignment length 17
ACTACGACTAGCTCAGG - Alpha
ACTACGGCAAGCACAGG - YYY

Alignment length 21
--ACTACGAC--TAGCTCAGG - Alpha
GGACTACGACAATAGCTCAGG - ZZZ

若使用 Bio.AlignIO.read() 或 Bio.AlignIO.parse() 而不帶 seq_count 參數，則會得到一個包含前兩個範例中所有六個記錄的單一比對結果。至於第三個範例，則會引發例外狀況，因為長度不同，導致無法將其轉換為單一比對結果。

如果檔案格式本身具有區塊結構，允許 Bio.AlignIO 直接判斷每個比對結果中的序列數量，則不需要 seq_count 參數。如果提供了此參數，且與檔案內容不符，則會引發錯誤。

請注意，這個選用的 seq_count 參數假設檔案中的每個比對結果都具有相同的序列數量。理論上，您可能會遇到更奇怪的情況，例如，一個 FASTA 檔案包含數個比對結果，每個比對結果都有不同的序列數量 — 雖然我很樂意聽到這類真實世界的範例。假設您無法取得更理想的檔案格式資料，則沒有直接的方式可使用 Bio.AlignIO 處理這種情況。在這種情況下，您可以考慮使用 Bio.SeqIO 讀取序列本身，並將它們分批組合在一起，以建立適當的比對結果。

寫入比對結果

我們已討論過使用 Bio.AlignIO.read() 和 Bio.AlignIO.parse() 進行比對輸入（讀取檔案），現在我們將探討用於比對輸出（寫入檔案）的 Bio.AlignIO.write()。這是一個接受三個參數的函數：一些 MultipleSeqAlignment 物件（為了回溯相容性，也可以使用過時的 Alignment 物件）、要寫入的控制代碼或檔案名稱，以及序列格式。

以下是一個範例，其中我們首先以較辛苦的方式（手動，而不是從檔案載入）建立一些 MultipleSeqAlignment 物件。請注意，我們建立了一些 SeqRecord 物件，以從中建構比對結果。

>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> from Bio.Align import MultipleSeqAlignment
>>> align1 = MultipleSeqAlignment(
...     [
...         SeqRecord(Seq("ACTGCTAGCTAG"), id="Alpha"),
...         SeqRecord(Seq("ACT-CTAGCTAG"), id="Beta"),
...         SeqRecord(Seq("ACTGCTAGDTAG"), id="Gamma"),
...     ]
... )
>>> align2 = MultipleSeqAlignment(
...     [
...         SeqRecord(Seq("GTCAGC-AG"), id="Delta"),
...         SeqRecord(Seq("GACAGCTAG"), id="Epsilon"),
...         SeqRecord(Seq("GTCAGCTAG"), id="Zeta"),
...     ]
... )
>>> align3 = MultipleSeqAlignment(
...     [
...         SeqRecord(Seq("ACTAGTACAGCTG"), id="Eta"),
...         SeqRecord(Seq("ACTAGTACAGCT-"), id="Theta"),
...         SeqRecord(Seq("-CTACTACAGGTG"), id="Iota"),
...     ]
... )
>>> my_alignments = [align1, align2, align3]

現在我們有 Alignment 物件的清單，我們將其寫入 PHYLIP 格式的檔案中：

>>> from Bio import AlignIO
>>> AlignIO.write(my_alignments, "my_example.phy", "phylip")

如果您在喜愛的文字編輯器中開啟此檔案，它應該如下所示：

 3 12
Alpha      ACTGCTAGCT AG
Beta       ACT-CTAGCT AG
Gamma      ACTGCTAGDT AG
 3 9
Delta      GTCAGC-AG
Epislon    GACAGCTAG
Zeta       GTCAGCTAG
 3 13
Eta        ACTAGTACAG CTG
Theta      ACTAGTACAG CT-
Iota       -CTACTACAG GTG

比較常見的做法是載入現有的比對結果，並儲存它，可能在進行一些簡單的操作之後，例如移除特定的列或欄。

假設您想知道 Bio.AlignIO.write() 函數將多少比對結果寫入控制代碼中？如果您的比對結果像上述範例一樣位於清單中，則您可以直接使用 len(my_alignments)，但是當您的記錄來自產生器/迭代器時，您無法這麼做。因此，Bio.AlignIO.write() 函數會傳回寫入檔案的比對結果數量。

注意 - 如果您告知 Bio.AlignIO.write() 函數寫入已存在的檔案，則舊檔案將在沒有任何警告的情況下被覆寫。

在序列比對檔案格式之間轉換

使用 Bio.AlignIO 在序列比對檔案格式之間轉換的方式，與使用 Bio.SeqIO 在序列檔案格式之間轉換的方式相同（請參閱 在序列檔案格式之間轉換 章節）。我們通常使用 Bio.AlignIO.parse() 載入比對結果，然後使用 Bio.AlignIO.write() 儲存它們 — 或直接使用 Bio.AlignIO.convert() 輔助函數。

在此範例中，我們將載入先前使用的 PFAM/Stockholm 格式檔案，並將其儲存為 Clustal W 格式的檔案：

>>> from Bio import AlignIO
>>> count = AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.aln", "clustal")
>>> print("Converted %i alignments" % count)
Converted 1 alignments

或者，使用 Bio.AlignIO.parse() 和 Bio.AlignIO.write()：

>>> from Bio import AlignIO
>>> alignments = AlignIO.parse("PF05371_seed.sth", "stockholm")
>>> count = AlignIO.write(alignments, "PF05371_seed.aln", "clustal")
>>> print("Converted %i alignments" % count)
Converted 1 alignments

Bio.AlignIO.write() 函數預期會收到多個比對物件。在上面的範例中，我們傳遞給它由 Bio.AlignIO.parse() 傳回的比對迭代器。

在這種情況下，我們知道檔案中只有一個比對結果，所以我們可以使用 Bio.AlignIO.read() 來取代，但是請注意，我們必須將此比對結果作為單一元素清單傳遞給 Bio.AlignIO.write()：

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> AlignIO.write([alignment], "PF05371_seed.aln", "clustal")

無論哪種方式，您都應該得到相同的新的 Clustal W 格式檔案「PF05371_seed.aln」，其內容如下：

CLUSTAL X (1.81) multiple sequence alignment


COATB_BPIKE/30-81                   AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSS
Q9T0Q8_BPIKE/1-52                   AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVS
COATB_BPI22/32-83                   DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSS
COATB_BPM13/24-72                   AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTS
COATB_BPZJ2/1-49                    AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFAS
Q9T0Q9_BPFD/1-49                    AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTS
COATB_BPIF1/22-73                   FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVS

COATB_BPIKE/30-81                   KA
Q9T0Q8_BPIKE/1-52                   RA
COATB_BPI22/32-83                   KA
COATB_BPM13/24-72                   KA
COATB_BPZJ2/1-49                    KA
Q9T0Q9_BPFD/1-49                    KA
COATB_BPIF1/22-73                   RA

或者，您可以建立一個 PHYLIP 格式的檔案，我們將其命名為「PF05371_seed.phy」：

>>> from Bio import AlignIO
>>> AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.phy", "phylip")

這次的輸出結果如下：

 7 52
COATB_BPIK AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIRLFKKFSS
Q9T0Q8_BPI AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIKLFKKFVS
COATB_BPI2 DGTSTATSYA TEAMNSLKTQ ATDLIDQTWP VVTSVAVAGL AIRLFKKFSS
COATB_BPM1 AEGDDP---A KAAFNSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPZJ AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFAS
Q9T0Q9_BPF AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPIF FAADDATSQA KAAFDSLTAQ ATEMSGYAWA LVVLVVGATV GIKLFKKFVS

           KA
           RA
           KA
           KA
           KA
           KA
           RA

原始 PHYLIP 比對檔案格式的一大缺點是，序列識別碼會嚴格截斷為十個字元。在此範例中，正如您所見，產生的名稱仍然是唯一的，但它們不太容易閱讀。因此，原始 PHYLIP 格式的更寬鬆變體現在已廣泛使用。

>>> from Bio import AlignIO
>>> AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.phy", "phylip-relaxed")

這次的輸出結果如下，使用較長的縮排，以允許完整顯示所有識別碼：

 7 52
COATB_BPIKE/30-81  AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIRLFKKFSS
Q9T0Q8_BPIKE/1-52  AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIKLFKKFVS
COATB_BPI22/32-83  DGTSTATSYA TEAMNSLKTQ ATDLIDQTWP VVTSVAVAGL AIRLFKKFSS
COATB_BPM13/24-72  AEGDDP---A KAAFNSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPZJ2/1-49   AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFAS
Q9T0Q9_BPFD/1-49   AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
COATB_BPIF1/22-73  FAADDATSQA KAAFDSLTAQ ATEMSGYAWA LVVLVVGATV GIKLFKKFVS

                   KA
                   RA
                   KA
                   KA
                   KA
                   KA
                   RA

如果您必須使用原始的嚴格 PHYLIP 格式，則可能需要以某種方式壓縮識別碼，或指定您自己的名稱或編號系統。以下程式碼片段會在儲存輸出之前操控記錄識別碼：

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> name_mapping = {}
>>> for i, record in enumerate(alignment):
...     name_mapping[i] = record.id
...     record.id = "seq%i" % i
...
>>> print(name_mapping)
{0: 'COATB_BPIKE/30-81', 1: 'Q9T0Q8_BPIKE/1-52', 2: 'COATB_BPI22/32-83', 3: 'COATB_BPM13/24-72', 4: 'COATB_BPZJ2/1-49', 5: 'Q9T0Q9_BPFD/1-49', 6: 'COATB_BPIF1/22-73'}
>>> AlignIO.write([alignment], "PF05371_seed.phy", "phylip")

此程式碼使用 Python 字典來記錄從新序列系統到原始識別碼的簡單對應：

{
    0: "COATB_BPIKE/30-81",
    1: "Q9T0Q8_BPIKE/1-52",
    2: "COATB_BPI22/32-83",
    # ...
}

以下是新的（嚴格）PHYLIP 格式輸出：

 7 52
seq0       AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIRLFKKFSS
seq1       AEPNAATNYA TEAMDSLKTQ AIDLISQTWP VVTTVVVAGL VIKLFKKFVS
seq2       DGTSTATSYA TEAMNSLKTQ ATDLIDQTWP VVTSVAVAGL AIRLFKKFSS
seq3       AEGDDP---A KAAFNSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
seq4       AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFAS
seq5       AEGDDP---A KAAFDSLQAS ATEYIGYAWA MVVVIVGATI GIKLFKKFTS
seq6       FAADDATSQA KAAFDSLTAQ ATEMSGYAWA LVVLVVGATV GIKLFKKFVS

           KA
           RA
           KA
           KA
           KA
           KA
           RA

一般來說，由於識別碼的限制，使用嚴格 PHYLIP 檔案格式不應該是您的首選。另一方面，使用 PFAM/Stockholm 格式可讓您記錄許多額外的註釋。

將您的比對物件作為格式化的字串

Bio.AlignIO 介面是以控制代碼為基礎，這表示如果您想要將比對結果以特定檔案格式放入字串中，則需要進行一些額外的工作（請參閱下方）。不過，您可能會偏好對比對物件呼叫 Python 內建的 format 函數。這會採用一個輸出格式規格作為單一參數，這是一個由 Bio.AlignIO 支援作為輸出格式的小寫字串。例如：

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> print(format(alignment, "clustal"))
CLUSTAL X (1.81) multiple sequence alignment


COATB_BPIKE/30-81                   AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSS
Q9T0Q8_BPIKE/1-52                   AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVS
COATB_BPI22/32-83                   DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSS
...

如果沒有輸出格式規格，format 會傳回與 str 相同的輸出。

如 格式方法 章節所述，SeqRecord 物件具有類似的方法，使用 Bio.SeqIO 支援的輸出格式。

在內部，format 正在以 StringIO 控制代碼呼叫 Bio.AlignIO.write()。如果您使用的是較舊版本的 Biopython，則可以在自己的程式碼中執行此操作：

>>> from io import StringIO
>>> from Bio import AlignIO
>>> alignments = AlignIO.parse("PF05371_seed.sth", "stockholm")
>>> out_handle = StringIO()
>>> AlignIO.write(alignments, out_handle, "clustal")
1
>>> clustal_data = out_handle.getvalue()
>>> print(clustal_data)
CLUSTAL X (1.81) multiple sequence alignment


COATB_BPIKE/30-81                   AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSS
Q9T0Q8_BPIKE/1-52                   AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVS
COATB_BPI22/32-83                   DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSS
COATB_BPM13/24-72                   AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTS
...

操控比對結果

現在我們已介紹完載入和儲存比對結果，我們將探討您還可以使用它們做些什麼。

切片比對結果

首先，在某些方面，比對物件的作用就像 Python list 的 SeqRecord 物件（列）。記住此模型，希望 len() 的動作（列數）和迭代（每一列作為 SeqRecord）都有意義：

>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> print("Number of rows: %i" % len(alignment))
Number of rows: 7
>>> for record in alignment:
...     print("%s - %s" % (record.seq, record.id))
...
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA - COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA - Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA - COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA - COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA - COATB_BPIF1/22-73

您也可以使用類似清單的 append 和 extend 方法將更多列新增至比對結果（作為 SeqRecord 物件）。記住清單的隱喻，簡單的比對結果切片也應該有意義 - 它會選取一些列，並傳回另一個比對物件：

>>> print(alignment)
Alignment with 7 rows and 52 columns
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRL...SKA COATB_BPIKE/30-81
AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKL...SRA Q9T0Q8_BPIKE/1-52
DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRL...SKA COATB_BPI22/32-83
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73
>>> print(alignment[3:7])
Alignment with 4 rows and 52 columns
AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49
AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49
FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73

如果您想要按欄選取，該怎麼辦？那些使用過 NumPy 矩陣或陣列物件的人不會對此感到驚訝 - 您會使用雙索引。

>>> print(alignment[2, 6])
T

使用兩個整數索引會提取單一字母，這是此操作的簡寫形式：

>>> print(alignment[2].seq[6])
T

您可以像這樣提取單一欄作為字串：

>>> print(alignment[:, 6])
TTT---T

您也可以選取一系列欄。例如，若要挑選出我們先前提取的相同三列，但只取它們的前六欄：

>>> print(alignment[3:6, :6])
Alignment with 3 rows and 6 columns
AEGDDP COATB_BPM13/24-72
AEGDDP COATB_BPZJ2/1-49
AEGDDP Q9T0Q9_BPFD/1-49

將第一個索引保留為 : 表示取所有列：

>>> print(alignment[:, :6])
Alignment with 7 rows and 6 columns
AEPNAA COATB_BPIKE/30-81
AEPNAA Q9T0Q8_BPIKE/1-52
DGTSTA COATB_BPI22/32-83
AEGDDP COATB_BPM13/24-72
AEGDDP COATB_BPZJ2/1-49
AEGDDP Q9T0Q9_BPFD/1-49
FAADDA COATB_BPIF1/22-73

這將引導我們找到一種移除區段的巧妙方式。請注意第 7、8 和 9 欄，這七個序列中有三個是間隙：

>>> print(alignment[:, 6:9])
Alignment with 7 rows and 3 columns
TNY COATB_BPIKE/30-81
TNY Q9T0Q8_BPIKE/1-52
TSY COATB_BPI22/32-83
--- COATB_BPM13/24-72
--- COATB_BPZJ2/1-49
--- Q9T0Q9_BPFD/1-49
TSQ COATB_BPIF1/22-73

同樣地，您可以切片以取得第九欄之後的所有內容：

>>> print(alignment[:, 9:])
Alignment with 7 rows and 43 columns
ATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA COATB_BPIKE/30-81
ATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA Q9T0Q8_BPIKE/1-52
ATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA COATB_BPI22/32-83
AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA COATB_BPM13/24-72
AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA COATB_BPZJ2/1-49
AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA Q9T0Q9_BPFD/1-49
AKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA COATB_BPIF1/22-73

現在，有趣的是，比對物件的加法是按欄進行的。這讓您可以透過這種方式移除一組欄：

>>> edited = alignment[:, :6] + alignment[:, 9:]
>>> print(edited)
Alignment with 7 rows and 49 columns
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA COATB_BPIKE/30-81
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA Q9T0Q8_BPIKE/1-52
DGTSTAATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA COATB_BPI22/32-83
AEGDDPAKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA COATB_BPM13/24-72
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA COATB_BPZJ2/1-49
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA Q9T0Q9_BPFD/1-49
FAADDAAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA COATB_BPIF1/22-73

比對加法的另一個常見用途是將數個不同基因的比對結果合併為元比對結果。不過請注意 — 識別碼需要比對（請參閱 新增 SeqRecord 物件 章節，瞭解新增 SeqRecord 物件的方式）。您可能會發現先依 ID 以字母順序排序比對列會有所幫助：

>>> edited.sort()
>>> print(edited)
Alignment with 7 rows and 49 columns
DGTSTAATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA COATB_BPI22/32-83
FAADDAAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA COATB_BPIF1/22-73
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA COATB_BPIKE/30-81
AEGDDPAKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA COATB_BPM13/24-72
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA COATB_BPZJ2/1-49
AEPNAAATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA Q9T0Q8_BPIKE/1-52
AEGDDPAKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA Q9T0Q9_BPFD/1-49

請注意，只有當兩個比對結果的列數相同時，才能將它們相加。

將比對結果作為陣列

根據您正在執行的操作，將比對物件轉換為字母陣列可能會更有用 — 而您可以使用 NumPy 執行此操作：

>>> import numpy as np
>>> from Bio import AlignIO
>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
>>> align_array = np.array(alignment)
>>> print("Array shape %i by %i" % align_array.shape)
Array shape 7 by 52
>>> align_array[:, :10]  
array([['A', 'E', 'P', 'N', 'A', 'A', 'T', 'N', 'Y', 'A'],
       ['A', 'E', 'P', 'N', 'A', 'A', 'T', 'N', 'Y', 'A'],
       ['D', 'G', 'T', 'S', 'T', 'A', 'T', 'S', 'Y', 'A'],
       ['A', 'E', 'G', 'D', 'D', 'P', '-', '-', '-', 'A'],
       ['A', 'E', 'G', 'D', 'D', 'P', '-', '-', '-', 'A'],
       ['A', 'E', 'G', 'D', 'D', 'P', '-', '-', '-', 'A'],
       ['F', 'A', 'A', 'D', 'D', 'A', 'T', 'S', 'Q', 'A']],...

請注意，這會將原始的 Biopython 比對物件和 NumPy 陣列以個別物件的形式保留在記憶體中 - 編輯其中一個物件不會更新另一個物件！

計算取代次數

比對的 substitutions 屬性會報告比對中字母彼此替換的頻率。這是透過取得比對中所有成對的列，計算兩個字母彼此對齊的次數，然後將所有配對的次數加總來計算的。例如：

>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> from Bio.Align import MultipleSeqAlignment
>>> msa = MultipleSeqAlignment(
...     [
...         SeqRecord(Seq("ACTCCTA"), id="seq1"),
...         SeqRecord(Seq("AAT-CTA"), id="seq2"),
...         SeqRecord(Seq("CCTACT-"), id="seq3"),
...         SeqRecord(Seq("TCTCCTC"), id="seq4"),
...     ]
... )
>>> print(msa)
Alignment with 4 rows and 7 columns
ACTCCTA seq1
AAT-CTA seq2
CCTACT- seq3
TCTCCTC seq4
>>> substitutions = msa.substitutions
>>> print(substitutions)
    A    C    T
A 2.0  4.5  1.0
C 4.5 10.0  0.5
T 1.0  0.5 12.0

由於配對的順序是任意的，因此計數會平均分配在對角線的上方和下方。例如，A 與 C 對齊的 9 次計數會儲存在位置 ['A', 'C'] 為 4.5，以及位置 ['C', 'A'] 為 4.5。當從這些計數計算取代矩陣時，這種排列方式有助於簡化數學計算，如第 取代矩陣 節中所述。

請注意，msa.substitutions 只包含比對中出現的字母的條目。您可以使用 select 方法來新增遺失字母的條目，例如：

>>> m = substitutions.select("ATCG")
>>> print(m)
    A    T    C   G
A 2.0  1.0  4.5 0.0
T 1.0 12.0  0.5 0.0
C 4.5  0.5 10.0 0.0
G 0.0  0.0  0.0 0.0

這也允許您變更字母順序。

計算摘要資訊

一旦您有了比對結果，您很可能會想要找出有關它的資訊。我們並沒有嘗試將所有可以產生比對資訊的函數都放在比對物件本身中，而是嘗試將功能分離到單獨的類別中，這些類別會對比對執行操作。

準備計算物件的摘要資訊很快就能完成。假設我們有一個名為 alignment 的比對物件，例如使用第 多序列比對物件 節中所述的 Bio.AlignIO.read(...) 讀取。我們只需要做的是建立一個將計算摘要資訊的物件：

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

summary_align 物件非常有用，並且會為您執行以下方便的操作：

1. 計算快速一致性序列 - 請參閱第 計算快速一致性序列 節

2. 取得比對的位置特定得分矩陣 - 請參閱第 位置特定得分矩陣 節

3. 計算比對的資訊含量 - 請參閱第 資訊含量 節

4. 產生比對中取代的相關資訊 - 第 取代矩陣 節詳細說明如何使用此資訊來產生取代矩陣。

計算快速一致性序列

第 計算摘要資訊 節中所述的 SummaryInfo 物件，提供了計算比對快速一致性的功能。假設我們有一個名為 summary_align 的 SummaryInfo 物件，我們可以透過執行以下操作來計算一致性：

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

顧名思義，這是一個非常簡單的一致性計算器，它只會加總一致性中每個點的所有殘基，如果最常見的值高於某個閾值，則會將該常見殘基加入一致性中。如果它沒有達到閾值，則會將模糊字元加入一致性中。傳回的一致性物件是一個 Seq 物件。

您可以透過傳遞可選參數來調整 dumb_consensus 的運作方式：

閾值

這是指定特定殘基在被加入前，在位置上必須有多常見的閾值。預設值為 \(0.7\)（表示 \(70\%\)）。

模糊字元

這是要使用的模糊字元。預設值為 'N'。

或者，您可以使用 alignment 屬性（請參閱第 取得新樣式的 Alignment 物件 節）將多序列比對物件 msa 轉換為新的 Alignment 物件（請參閱第 Alignment 物件 節）：

>>> alignment = msa.alignment

然後，您可以建立一個 Motif 物件（請參閱第 Motif 物件 節）：

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

並取得快速一致性序列：

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

motif.counts.calculate_consensus 方法（請參閱第 取得一致性序列 節）可讓您詳細指定應如何計算一致性序列。例如：

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

位置特定得分矩陣

位置特定得分矩陣 (PSSM) 以不同於一致性的方式摘要比對資訊，並且可能對不同的任務有用。基本上，PSSM 是一個計數矩陣。對於比對中的每個欄位，會計數並加總每個字母表的字母數量。總數會相對於沿著左軸的某些代表性序列顯示。此序列可能是共識序列，但也可能是比對中的任何序列。

例如，對於上面的比對：

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

我們使用以下方法取得一個沿側具有一致性序列的 PSSM：

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

其中我們在計算 PSSM 時會忽略任何 N 模糊殘基。

關於這一點，應該注意兩件事：

1. 為了保持字母表的嚴格性，您只能在 PSSM 的頂部包含比對物件字母表中的字元。間隙不會包含在 PSSM 的頂軸上。

2. 沿著軸的左側顯示的序列不需要是一致性。例如，如果您想要沿著此軸顯示比對中的第二個序列，您需要執行以下操作：

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

上面的命令會傳回一個 PSSM 物件。您可以透過像 your_pssm[序列編號][殘基計數名稱] 這樣的方式下標來存取 PSSM 的任何元素。例如，要取得上述 PSSM 中第二個元素的 'A' 殘基的計數，您需要執行：

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

PSSM 類的結構希望可以方便地存取元素並美觀地列印矩陣。

或者，您可以使用 alignment 屬性（請參閱第 取得新樣式的 Alignment 物件 節）將多序列比對物件 msa 轉換為新的 Alignment 物件（請參閱第 Alignment 物件 節）：

>>> alignment = msa.alignment

然後，您可以建立一個 Motif 物件（請參閱第 Motif 物件 節）：

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

並取得每個位置中每個核苷酸的計數：

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

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

資訊含量

一個可能對演化保守性有用的衡量標準是序列的資訊含量。

可以在 http://www.lecb.ncifcrf.gov/~toms/paper/primer/ 找到針對分子生物學家的資訊理論的有用介紹。就我們的目的而言，我們將研究一致性序列或一致性序列一部分的資訊含量。我們使用以下公式計算多序列比對中特定欄位的資訊含量：

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

其中：

• \(IC_{j}\) – 比對中第 \(j\) 欄的資訊含量。

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

• \(P_{ij}\) – 特定字母 \(i\) 在第 \(j\) 欄的頻率（例如，如果 G 在比對欄中出現 6 次中的 3 次，則此值為 0.5）

• \(Q_{i}\) – 字母 \(i\) 的預期頻率。這是一個可選的參數，是否使用由使用者自行決定。預設情況下，對於蛋白質字母表，它會自動設定為 \(0.05 = 1/20\)，對於核酸字母表，則設定為 \(0.25 = 1/4\)。這是為了在沒有任何先前分配假設的情況下取得資訊含量。當假設先前分配或使用非標準字母表時，您應該提供 \(Q_{i}\) 的值。

好的，現在我們已經了解了在 Biopython 中計算的資訊含量的概念，讓我們看看如何取得比對特定區域的資訊含量。

首先，我們需要使用我們的比對來取得一個比對摘要物件，我們假設它叫做 summary_align（請參閱第 計算摘要資訊 節，以取得如何取得此物件的指示）。一旦我們取得了此物件，計算區域的資訊含量就如同以下操作一樣簡單：

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

現在，info_content 將包含相對於預期頻率的 [2:6] 區域的相對資訊含量。

傳回的值是使用公式中的對數底數 2 作為底數計算的。您可以透過將參數 log_base 作為您想要的底數來修改此值：

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

依預設，在計算該欄的相對資訊欄時，不會考慮欄中頻率為 0 的核苷酸或氨基酸殘基。如果這不是想要的结果，您可以使用 pseudo_count 來代替。

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

在這種情況下，特定字母 \(i\) 在第 \(j\) 欄的觀察頻率 \(P_{ij}\) 計算如下：

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

其中：

• \(k\) – 您作為引數傳遞的假計數。

• \(k\) – 您作為引數傳遞的假計數。

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

好的，現在您已經準備好計算資訊含量了。如果您想嘗試將其應用於一些真實世界的問題，最好深入研究有關資訊含量的文獻，以了解如何使用它。希望您的深入研究不會發現編碼此函數時發生的任何錯誤！

取得新樣式的 Alignment 物件

使用 alignment 屬性，從舊樣式的 MultipleSeqAlignment 物件建立新的 Alignment 物件（請參閱第 Alignment 物件 節）：

>>> type(msa)
<class 'Bio.Align.MultipleSeqAlignment'>
>>> print(msa)
Alignment with 4 rows and 7 columns
ACTCCTA seq1
AAT-CTA seq2
CCTACT- seq3
TCTCCTC seq4
>>> alignment = msa.alignment
>>> type(alignment)
<class 'Bio.Align.Alignment'>
>>> print(alignment)
seq1              0 ACTCCTA 7
seq2              0 AAT-CTA 6
seq3              0 CCTACT- 6
seq4              0 TCTCCTC 7

請注意，alignment 屬性會建立並回傳一個新的 Alignment 物件，該物件與建立 Alignment 物件時儲存在 MultipleSeqAlignment 物件中的資訊一致。在呼叫 alignment 屬性之後對 MultipleSeqAlignment 所做的任何變更，都不會傳播到 Alignment 物件。不過，您當然可以再次呼叫 alignment 屬性，以建立一個與更新後的 MultipleSeqAlignment 物件一致的新 Alignment 物件。

從多序列比對計算取代矩陣

您可以從比對建立自己的取代矩陣。在此範例中，我們將首先從 Clustalw 檔案 protein.aln（也可在線上取得 此處）讀取蛋白質序列比對。

>>> from Bio import AlignIO
>>> filename = "protein.aln"
>>> msa = AlignIO.read(filename, "clustal")

章節 ClustalW 包含更多關於執行此操作的資訊。

比對的 substitutions 屬性會儲存不同殘基互相取代的次數。

>>> observed_frequencies = msa.substitutions

為了使範例更易讀，我們將僅選擇具有極性帶電荷側鏈的胺基酸。

>>> observed_frequencies = observed_frequencies.select("DEHKR")
>>> print(observed_frequencies)
       D      E      H      K      R
D 2360.0  255.5    7.5    0.5   25.0
E  255.5 3305.0   16.5   27.0    2.0
H    7.5   16.5 1235.0   16.0    8.5
K    0.5   27.0   16.0 3218.0  116.5
R   25.0    2.0    8.5  116.5 2079.0

其他胺基酸的列和欄已從矩陣中移除。

接下來，我們對矩陣進行正規化。

>>> import numpy as np
>>> observed_frequencies /= np.sum(observed_frequencies)

對列或欄求和會得出每個殘基出現的相對頻率。

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

>>> sum(residue_frequencies) == 1.0
True

然後，殘基對的預期頻率為

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

在此，residue_frequencies[:, None] 會建立一個二維陣列，其中包含一個單一欄，其值為 residue_frequencies，而 residue_frequencies[None, :] 則是一個二維陣列，這些值為單一列。取它們的點積（內積）會建立一個預期頻率的矩陣，其中每個項目都由兩個 residue_frequencies 值相乘而成。例如，expected_frequencies['D', 'E'] 等於 residue_frequencies['D'] * residue_frequencies['E']。

現在，我們可以透過將觀察到的頻率除以預期的頻率並取對數來計算對數優勢矩陣。

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

此矩陣可用作執行比對時的取代矩陣。例如，

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

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

比對工具

有很多演算法可用於比對序列，包括成對比對和多序列比對。這些計算相對較慢，通常您不會想要在 Python 中編寫此類演算法。對於成對比對，您可以使用 Biopython 的 PairwiseAligner（請參閱章節 成對序列比對），它以 C 語言實作，因此速度很快。或者，您可以透過從 Python 呼叫外部比對程式來執行它。通常你會

1. 準備未比對序列的輸入檔，通常這會是一個 FASTA 檔案，您可以使用 Bio.SeqIO 來建立（請參閱章節 序列輸入/輸出）。

2. 透過使用 Python 的 subprocess 模組執行其命令來執行比對程式。

3. 從工具讀取輸出，即您已比對的序列，通常使用 Bio.AlignIO（請參閱本章稍早）。

在此，我們將展示此工作流程的幾個範例。

ClustalW

ClustalW 是一個流行的多序列比對命令列工具（還有一個稱為 ClustalX 的圖形介面）。在嘗試從 Python 內部使用 ClustalW 之前，您應該先嘗試在命令列上手動執行 ClustalW 工具，以熟悉其他選項。

對於最基本的使用方式，您只需要有一個 FASTA 輸入檔，例如 opuntia.fasta（可在線上取得或在 Biopython 原始碼的 Doc/examples 子目錄中）。這是一個小的 FASTA 檔案，其中包含七個仙人掌 DNA 序列（來自仙人掌科 Opuntia）。預設情況下，ClustalW 將產生一個比對和引導樹檔案，其名稱基於輸入 FASTA 檔案，在此範例中為 opuntia.aln 和 opuntia.dnd，但您可以覆寫此設定或使其明確。

>>> import subprocess
>>> cmd = "clustalw2 -infile=opuntia.fasta"
>>> results = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, text=True)

請注意，這裡我們已將可執行檔名稱指定為 clustalw2，表示我們已安裝第二版，其檔案名稱與第一版（預設的 clustalw）不同。幸運的是，兩個版本都支援命令列上的相同引數集（而且實際上，功能應該相同）。

您可能會發現，即使您已安裝 ClustalW，上述命令也無法運作 – 您可能會收到關於「找不到命令」的訊息（尤其是在 Windows 上）。這表示 ClustalW 可執行檔不在您的 PATH 中（一個環境變數，要搜尋的目錄清單）。您可以更新 PATH 設定以包含 ClustalW 工具的複本位置（如何執行此操作取決於您的作業系統），或直接輸入工具的完整路徑。請記住，在 Python 字串中，\n 和 \t 預設會被解譯為換行符號和 Tab – 這就是為什麼我們在開頭放置字母「r」作為原始字串，不會以這種方式轉換。當指定 Windows 樣式檔案名稱時，這通常是很好的做法。

>>> import os
>>> clustalw_exe = r"C:\Program Files\new clustal\clustalw2.exe"
>>> assert os.path.isfile(clustalw_exe), "Clustal W executable missing"
>>> cmd = clustalw_exe + " -infile=opuntia.fasta"
>>> results = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, text=True)

現在，此時了解命令列工具「如何運作」會有所幫助。當您在命令列執行工具時，它通常會將文字輸出直接列印到螢幕上。此文字可以透過兩個「管道」（稱為標準輸出（正常結果）和標準錯誤（用於錯誤訊息和偵錯訊息））擷取或重新導向。還有標準輸入，它是饋送到工具的任何文字。這些名稱縮寫為 stdin、stdout 和 stderr。當工具完成時，它會有一個傳回碼（一個整數），按照慣例，成功時傳回零，而非零傳回碼表示發生錯誤。

在上面的 ClustalW 範例中，在命令列執行時，所有重要的輸出都會直接寫入輸出檔案。當您等待時，通常列印到螢幕上的所有內容都會擷取到 results.stdout 和 results.stderr 中，而傳回碼則儲存在 results.returncode 中。

我們關心的是兩個輸出檔案，即比對和引導樹。我們沒有告訴 ClustalW 要使用哪些檔案名稱，但它預設會根據輸入檔案挑選名稱。在此情況下，輸出應該位於 opuntia.aln 檔案中。您現在應該能夠算出如何使用 Bio.AlignIO 讀取比對。

>>> from Bio import AlignIO
>>> align = AlignIO.read("opuntia.aln", "clustal")
>>> print(align)
Alignment with 7 rows and 906 columns
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273285|gb|AF191659.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273284|gb|AF191658.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273287|gb|AF191661.1|AF191
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273286|gb|AF191660.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273290|gb|AF191664.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273289|gb|AF191663.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273291|gb|AF191665.1|AF191

如果您有興趣（這是本章主要內容的題外話），ClustalW 建立的 opuntia.dnd 檔案只是一個標準的 Newick 樹狀檔案，而 Bio.Phylo 可以剖析這些檔案。

>>> from Bio import Phylo
>>> tree = Phylo.read("opuntia.dnd", "newick")
>>> Phylo.draw_ascii(tree)
                             _______________ gi|6273291|gb|AF191665.1|AF191665
  __________________________|
 |                          |   ______ gi|6273290|gb|AF191664.1|AF191664
 |                          |__|
 |                             |_____ gi|6273289|gb|AF191663.1|AF191663
 |
_|_________________ gi|6273287|gb|AF191661.1|AF191661
 |
 |__________ gi|6273286|gb|AF191660.1|AF191660
 |
 |    __ gi|6273285|gb|AF191659.1|AF191659
 |___|
     | gi|6273284|gb|AF191658.1|AF191658

章節 使用 Bio.Phylo 進行系統發生學分析 更深入地介紹了 Biopython 對系統發生樹的支援。

MUSCLE

MUSCLE 是一個比 ClustalW 更新的多序列比對工具。與之前一樣，我們建議您先嘗試從命令列使用 MUSCLE，然後再嘗試從 Python 執行它。

對於最基本的使用方式，您只需要有一個 FASTA 輸入檔，例如 opuntia.fasta（可在線上取得或在 Biopython 原始碼的 Doc/examples 子目錄中）。然後，您可以告知 MUSCLE 讀取此 FASTA 檔案，並將比對寫入名為 opuntia.txt 的輸出檔案。

>>> import subprocess
>>> cmd = "muscle -align opuntia.fasta -output opuntia.txt"
>>> results = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, text=True)

MUSCLE 將以 FASTA 檔案格式輸出比對（使用帶間隙的序列）。Bio.AlignIO 模組能夠使用 format="fasta" 讀取此比對。

>>> from Bio import AlignIO
>>> align = AlignIO.read("opuntia.txt", "fasta")
>>> print(align)
Alignment with 7 rows and 906 columns
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273289|gb|AF191663.1|AF191663
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273291|gb|AF191665.1|AF191665
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273290|gb|AF191664.1|AF191664
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273287|gb|AF191661.1|AF191661
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273286|gb|AF191660.1|AF191660
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273285|gb|AF191659.1|AF191659
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273284|gb|AF191658.1|AF191658

您也可以設定其他選用參數；請參閱 MUSCLE 的內建說明以取得詳細資料。

EMBOSS needle 和 water

EMBOSS 套件包含了 water 和 needle 工具，分別用於 Smith-Waterman 演算法的局部比對和 Needleman-Wunsch 演算法的全域比對。這些工具共享相同的介面風格，因此在這兩者之間切換非常容易 – 在這裡我們將只使用 needle。

假設您想要在兩個序列之間進行全域成對比對，這兩個序列以 FASTA 格式準備，如下所示：

>HBA_HUMAN
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR

一個在檔案 alpha.faa 中，另一個在檔案 beta.faa 中

>HBB_HUMAN
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPK
VKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFG
KEFTPPVQAAYQKVVAGVANALAHKYH

您可以在 Biopython 原始碼的 Doc/examples/ 目錄下找到這些範例檔案的副本。

使用 needle 將這兩個序列互相比對的命令如下：

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

何不嘗試在命令提示字元下手動執行此操作？您應該會看到它進行了成對比較，並將輸出記錄在 needle.txt 檔案中（使用預設的 EMBOSS 比對檔案格式）。

即使您已安裝 EMBOSS，執行此命令也可能無法運作 – 您可能會收到有關「找不到命令」的訊息（尤其是在 Windows 上）。這可能表示 EMBOSS 工具不在您的 PATH 環境變數中。您可以更新 PATH 設定，或直接使用工具的完整路徑，例如：

C:\EMBOSS\needle.exe -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5

接下來，我們想使用 Python 來為我們執行此命令。如上所述，為了完全控制，我們建議您使用 Python 的內建模組 subprocess。

>>> import sys
>>> import subprocess
>>> cmd = "needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5"
>>> results = subprocess.run(
...     cmd,
...     stdout=subprocess.PIPE,
...     stderr=subprocess.PIPE,
...     text=True,
...     shell=(sys, platform != "win32"),
... )
>>> print(results.stdout)

>>> print(results.stderr)
Needleman-Wunsch global alignment of two sequences

接下來，我們可以像本章前面討論的那樣，使用 Bio.AlignIO 以 emboss 格式載入輸出檔案。

>>> from Bio import AlignIO
>>> align = AlignIO.read("needle.txt", "emboss")
>>> print(align)
Alignment with 2 rows and 149 columns
MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY...KYR HBA_HUMAN
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRF...KYH HBB_HUMAN

在這個範例中，我們指示 EMBOSS 將輸出寫入檔案，但您 *可以* 指示它將輸出寫入 stdout（如果您不想要捨棄臨時輸出檔案，這會很有用 – 使用 outfile=stdout 參數）。

>>> cmd = "needle -outfile=stdout -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5"
>>> child = subprocess.Popen(
...     cmd,
...     stdout=subprocess.PIPE,
...     stderr=subprocess.PIPE,
...     text=True,
...     shell=(sys.platform != "win32"),
... )
>>> align = AlignIO.read(child.stdout, "emboss")
>>> print(align)
Alignment with 2 rows and 149 columns
MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY...KYR HBA_HUMAN
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRF...KYH HBB_HUMAN

類似地，也可以從 stdin 讀取 *其中一個* 輸入（例如，asequence="stdin"）。

這只是 needle 和 water 功能的冰山一角。一個有用的技巧是，第二個檔案可以包含多個序列（比如五個），然後 EMBOSS 將執行五個成對比對。