計算分子生物學的 Python 工具

Biopython 版本 1.84
© 2024。保留所有權利。
在 GitHub 上編輯此頁面

SeqIO 簡介

此頁面描述 Bio.SeqIO，這是 BioPython 1.43 及更高版本的標準序列輸入/輸出介面。 有關實作細節，請參閱 SeqIO 開發頁面。

Python 初學者可能會發現 Peter 的入門 Biopython 工作坊 很有用，該工作坊從使用 SeqIO 處理序列檔案開始。

教學（PDF）中有一整章關於 Bio.SeqIO 的內容，雖然有一些重疊，但除了這個 WIKI 頁面之外，也值得閱讀。 還有 API 文件（您可以在線上閱讀，或從 Python 中使用 help 命令閱讀）。

目標

Bio.SeqIO 提供了一個簡單統一的介面，用於輸入和輸出各種序列檔案格式（包括多序列比對），但只處理作為 SeqRecord 物件的序列。 有一個姊妹介面 Bio.AlignIO，用於直接處理序列比對檔案作為 Alignment 物件。

這個設計部分靈感來自於 BioPerl 的 SeqIO 的簡潔性。 從長遠來看，我們希望與 BioPerl 令人印象深刻的受支援的 序列檔案格式和多重比對格式列表相匹配。

請注意，在 Biopython 中包含 Bio.SeqIO（和 Bio.AlignIO）確實會導致在處理某些檔案格式時出現一些重複或選擇。 例如，Bio.Nexus 也會從 Nexus 檔案讀取序列 - 但 Bio.Nexus 還可以做更多的事情，例如讀取 Nexus 檔案中的任何系統發生樹。

我的願景是，對於處理序列資料，您應該嘗試將 Bio.SeqIO 作為您的首選。 除非您有一些非常具體的要求，否則我希望這應該足夠了。

Peter

檔案格式

下表列出了 Bio.SeqIO 可以讀取、寫入和索引的檔案格式，以及首次支援此功能的 Biopython 版本（或 git 表示此功能在我們最新的開發程式碼中受支援）。 格式名稱是一個簡單的小寫字串。 在可能的情況下，我們使用與 BioPerl 的 SeqIO 和 EMBOSS 相同的名稱。

格式名稱 讀取 寫入 索引 註釋
abi 1.58 否 不適用 讀取 ABI「Sanger」毛細管序列追蹤檔案，包括鹼基讀取的 PHRED 品質分數。 這允許 ABI 轉換為 FASTQ。 請注意，每個 ABI 檔案只包含一個序列（因此索引檔案沒有意義）。
abi-trim 1.71 否 不適用 與「abi」相同，但使用 Mott 演算法進行品質修剪。
ace 1.47 否 1.52 從 ACE 組裝檔案讀取重疊群序列。 內部使用 Bio.Sequencing.Ace
cif-atom 1.73 否 否 使用 Bio.PDB.MMCIFParser 根據原子座標確定結構中出現的（部分）蛋白質序列。
cif-seqres 1.73 否 否 讀取巨分子結晶資訊檔案 (mmCIF) 檔案，以確定由 _pdbx_poly_seq_scheme 記錄定義的完整蛋白質序列。
clustal 1.43 1.43 否 Clustal X 和 Clustal W 的比對格式。
embl 1.43 1.54 1.52 EMBL 平面檔案格式。 內部使用 Bio.GenBank。
fasta 1.43 1.43 1.52 這指的是 Bill Pearson 的 FASTA 工具引入的輸入 FASTA 檔案格式，其中每個記錄都以「>」行開始。
fasta-2line 1.71 1.71 否 FASTA 格式變體，沒有換行且每個記錄恰好有兩行。
fastq-sanger 或 fastq 1.50 1.50 1.52 FASTQ 檔案有點像 FASTA 檔案，但還包括定序品質。 在 Biopython 中，「fastq」（或別名「fastq-sanger」）指的是 Sanger 樣式的 FASTQ 檔案，該檔案使用 33 的 ASCII 偏移量對 PHRED 品質進行編碼。 另請參閱早期 Solexa/Illumina 管道中使用的不相容的「fastq-solexa」和「fastq-illumina」變體，Illumina 管道 1.8 產生 Sanger FASTQ。
fastq-solexa 1.50 1.50 1.52 在 Biopython 中，「fastq-solexa」指的是原始的 Solexa/Illumina 樣式的 FASTQ 檔案，該檔案使用 64 的 ASCII 偏移量對 Solexa 品質進行編碼。 另請參閱我們稱之為「fastq-illumina」的格式。
fastq-illumina 1.51 1.51 1.52 在 Biopython 中，「fastq-illumina」指的是早期的 Solexa/Illumina 樣式的 FASTQ 檔案（來自管道版本 1.3 到 1.7），該檔案使用 64 的 ASCII 偏移量對 PHRED 品質進行編碼。 對於品質良好的讀取，PHRED 和 Solexa 分數大致相等，因此「fastq-solexa」和「fastq-illumina」變體幾乎等效。
gck 1.75 否 否 Gene Construction Kit 使用的原生格式。
genbank 或 gb 1.43 1.48 / 1.51 1.52 GenBank 或 GenPept 平面檔案格式。 內部使用 Bio.GenBank 進行解析。 Biopython 1.48 到 1.50 寫入僅具有最少註釋的基本 GenBank 檔案，而 1.51 及更高版本也會寫入特徵表。
ig 1.47 否 1.52 這指的是 IntelliGenetics 檔案格式，顯然與 MASE 比對格式相同。
imgt 1.56 1.56 1.56 這指的是 EMBL 純文字檔案格式的 IMGT 變體。
nexus 1.43 1.48 否 NEXUS 多重比對格式，也稱為 PAUP 格式。 內部使用 Bio.Nexus。
pdb-seqres 1.61 否 否 讀取蛋白質資料庫 (PDB) 檔案以確定標頭中顯示的完整蛋白質序列（不依賴 Bio.PDB 和 NumPy）。
pdb-atom 1.61 否 否 使用 Bio.PDB 根據檔案的原子座標部分確定結構中出現的（部分）蛋白質序列（需要 NumPy）。
phd 1.46 1.52 1.52 PHD 檔案是 PHRED 的輸出，由 PHRAP 和 CONSED 用於輸入。 內部使用 Bio.Sequencing.Phd。
phylip 1.43 1.43 否 PHYLIP 檔案。 將名稱截斷為 10 個字元。
pir 1.48 1.71 1.52 由國家生物醫學研究基金會 (NBRF) 為 蛋白質資訊資源 (PIR) 資料庫引入的「類似 FASTA」格式，現在是 UniProt 的一部分。
seqxml 1.58 1.58 否 簡單的 序列 XML 檔案格式。
sff 1.54 1.54 1.54 Roche 454 和 IonTorrent/IonProton 定序機產生的標準流程圖格式 (SFF) 二進位檔案。
sff-trim 1.54 否 1.54 套用檔案中列出的修剪的標準流程圖格式。
snapgene 1.75 否 否 SnapGene 使用的原生格式。
stockholm 1.43 1.43 否 Stockholm 比對格式也稱為 PFAM 格式。
swiss 1.43 否 1.52 Swiss-Prot 也稱為 UniProt 格式。 內部使用 Bio.SwissProt。 另請參閱 UniProt XML 格式。
tab 1.48 1.48 1.52 簡單的兩欄索引標籤分隔序列檔案，其中每行都包含記錄的識別碼和序列。 例如，Aligent 的 eArray 軟體在以最小的索引標籤分隔文字檔案儲存微陣列探針時使用此格式。
qual 1.50 1.50 1.52 Qual 檔案有點像 FASTA 檔案，但不是序列，而是記錄以空格分隔的整數定序值作為 PHRED 品質分數。 一對匹配的 FASTA 和 QUAL 檔案通常用作單個 FASTQ 檔案的替代方案。
uniprot-xml 1.56 否 1.56 UniProt XML 格式，是純文字 Swiss-Prot 格式的後繼格式。
xdna 1.75 1.75 否 Christian Marck 的 DNA Strider 和 Serial Cloner 使用的原生格式。
         

使用 Bio.SeqIO，您可以將序列比對檔案格式視為任何其他序列檔案，但新的 Bio.AlignIO 模組旨在直接處理此類比對檔案。 您也可以將任何檔案格式的一組 SeqRecord 物件轉換為比對 - 前提是它們的長度都相同。 請注意，當使用 Bio.SeqIO 將序列寫入比對檔案格式時，所有（帶間隙）序列的長度應相同。

序列輸入

主要函數是 Bio.SeqIO.parse()，它接受檔案處理常式（或檔案名稱）和格式名稱，並傳回 SeqRecord 迭代器。 這允許您執行以下操作，例如

from Bio import SeqIO

for record in SeqIO.parse("example.fasta", "fasta"):
    print(record.id)

或使用處理常式

from Bio import SeqIO

with open("example.fasta") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

在上面的範例中，我們使用內建的 python 函數 open 開啟檔案。 with-語句確保在讀取檔案後正確關閉檔案。 如果您只使用檔案名稱，則這一切都應該會自動發生。

請注意，您必須明確指定檔案格式，這與 BioPerl 的 SeqIO 不同，後者可以嘗試使用副檔名和/或檔案內容來猜測。 請參閱明確勝於隱含（Python 之禪）。

如果您有不同類型的檔案，例如 Clustalw 比對檔案（例如 opuntia.aln，其中包含七個序列），唯一的區別是您指定 "clustal" 而不是 "fasta"

from Bio import SeqIO

with open("opuntia.aln") as handle:
    for record in SeqIO.parse(handle, "clustal"):
        print(record.id)

當您只需要依序取得檔案中找到的記錄時，迭代器非常有用。 對於某些任務，您可能需要以任何順序隨機存取記錄。 在這種情況下，請使用內建的 python list() 函數將迭代器轉換為列表

from Bio import SeqIO

records = list(SeqIO.parse("example.fasta", "fasta"))
print(records[0].id)  # first record
print(records[-1].id)  # last record

另一個常見的任務是依據某個識別符來索引您的記錄。對於小型檔案，我們有一個函數 Bio.SeqIO.to_dict()，可以將 SeqRecord 迭代器（或列表）轉換為字典（在記憶體中）。

from Bio import SeqIO

record_dict = SeqIO.to_dict(SeqIO.parse("example.fasta", "fasta"))
print(record_dict["gi:12345678"])  # use any record ID

函數 Bio.SeqIO.to_dict() 預設會使用記錄 ID 作為字典的鍵，但您可以使用其可選參數 key_function 指定您想要的任何映射。

對於較大的檔案，不可能將所有內容都保存在記憶體中，因此 Bio.SeqIO.to_dict 不適用。Biopython 1.52 或更新的版本包含了 Bio.SeqIO.index 函數來處理這種情況，但您也可以考慮使用 BioSQL。

from Bio import SeqIO

record_dict = SeqIO.index("example.fasta", "fasta")
print(record_dict["gi:12345678"])  # use any record ID

Biopython 1.45 引入了另一個函數 Bio.SeqIO.read()，它和 Bio.SeqIO.parse() 一樣，需要一個句柄和格式。它適用於當句柄僅包含一筆記錄的情況，並將其作為單個 SeqRecord 物件返回。如果沒有記錄，或有多於一筆的記錄，則會引發例外。

from Bio import SeqIO

record = SeqIO.read("single.fasta", "fasta")

對於您只想取得第一筆記錄（並樂於忽略後續記錄）的相關情況，您可以使用內建的 Python 函數 next。

from Bio import SeqIO

first_record = next(SeqIO.parse("example.fasta", "fasta"))

序列輸出

要將記錄寫入檔案，請使用函數 Bio.SeqIO.write()，它接受 SeqRecord 迭代器（或列表）、輸出句柄（或檔案名稱）和格式字串。

from Bio import SeqIO

sequences = ...  # add code here
with open("example.fasta", "w") as output_handle:
    SeqIO.write(sequences, output_handle, "fasta")

或者

from Bio import SeqIO

sequences = ...  # add code here
SeqIO.write(sequences, "example.fasta", "fasta")

在後續關於檔案格式轉換的章節中有更多範例。

請注意，如果您要寫入對齊檔案格式，則所有序列的長度必須相同。

如果您將序列作為 SeqRecord 迭代器提供，那麼對於像 Fasta 或 GenBank 這樣的循序檔案格式，可以逐個寫入記錄。因為一次只建立一筆記錄，所以只需要很少的記憶體。請參閱下方過濾一組記錄的範例。

另一方面，對於像 Clustal 這樣交錯或非循序的檔案格式，Bio.SeqIO.write() 函數將被迫自動將迭代器轉換為列表。這將破壞使用產生器/迭代器方法所帶來的任何潛在記憶體節省。

檔案格式轉換

假設您有一個想要轉換為 Fasta 檔案的 GenBank 檔案。例如，讓我們考慮 cor6_6.gb 這個檔案，它包含在 Biopython 單元測試的 GenBank 目錄下。

您可以使用 Bio.SeqIO.parse() 函數像這樣讀取檔案。

from Bio import SeqIO

with open("cor6_6.gb") as input_handle:
    for record in SeqIO.parse(input_handle, "genbank"):
        print(record)

請注意，此檔案包含六筆記錄。現在，不要印出記錄，而是將 SeqRecord 迭代器傳遞給 Bio.SeqIO.write() 函數，以將這個 GenBank 檔案轉換為 Fasta 檔案。

from Bio import SeqIO

with open("cor6_6.gb") as input_handle, open(
    "cor6_6.fasta", "w"
) as output_handle:
    sequences = SeqIO.parse(input_handle, "genbank")
    count = SeqIO.write(sequences, output_handle, "fasta")

print("Converted %i records" % count)

或者，更簡潔地使用 Bio.SeqIO.convert() 函數（在 Biopython 1.52 或更新版本中），只需：

from Bio import SeqIO

count = SeqIO.convert("cor6_6.gb", "genbank", "cor6_6.fasta", "fasta")
print("Converted %i records" % count)

在這個範例中，GenBank 檔案的開頭是這樣的：

LOCUS       ATCOR66M      513 bp    mRNA            PLN       02-MAR-1992
DEFINITION  A.thaliana cor6.6 mRNA.
ACCESSION   X55053
VERSION     X55053.1  GI:16229
...

產生的 Fasta 檔案看起來像這樣：

>X55053.1 A.thaliana cor6.6 mRNA.
AACAAAACACACATCAAAAACGATTTTACAAGAAAAAAATA...
...

請注意，Fasta 檔案只能儲存識別符、描述和序列。

透過更改格式字串，該程式碼可用於在任何支援的檔案格式之間進行轉換。

範例

輸入/輸出範例 - 依序列長度過濾

雖然您可能只想轉換檔案（如上所示），但更實際的範例是以某種方式操作或過濾資料。

例如，讓我們將所有少於 300 個核苷酸的「短」序列儲存到 Fasta 檔案中。

from Bio import SeqIO

short_sequences = []  # Setup an empty list
for record in SeqIO.parse("cor6_6.gb", "genbank"):
    if len(record.seq) < 300:
        # Add this record to our list
        short_sequences.append(record)

print("Found %i short sequences" % len(short_sequences))

SeqIO.write(short_sequences, "short_seqs.fasta", "fasta")

如果您了解列表推導式，那麼您也可以像這樣編寫上面的範例：

from Bio import SeqIO

input_seq_iterator = SeqIO.parse("cor6_6.gb", "genbank")

# Build a list of short sequences:
short_sequences = [record for record in input_seq_iterator if len(record.seq) < 300]

print("Found %i short sequences" % len(short_sequences))

SeqIO.write(short_sequences, "short_seqs.fasta", "fasta")

我不認為這實際上更容易理解，但它的確比較簡短。

然而，如果您正在處理包含數千筆記錄的非常大的檔案，您可以使用產生器表達式來獲益。這可以避免在記憶體中建立整個所需記錄的列表。

from Bio import SeqIO

input_seq_iterator = SeqIO.parse("cor6_6.gb", "genbank")
short_seq_iterator = (record for record in input_seq_iterator if len(record.seq) < 300)

SeqIO.write(short_seq_iterator, "short_seqs.fasta", "fasta")

請記住，對於像 Fasta 或 GenBank 這樣的循序檔案格式，Bio.SeqIO.write() 將接受 SeqRecord 迭代器。上述程式碼的優點是任何時候都只有一筆記錄在記憶體中。

然而，如輸出章節所述，對於像 Clustal 這樣的非循序檔案格式，Bio.SeqIO.write() 會被迫自動將迭代器轉換為列表，因此這個優勢就消失了。

如果這一切都令人困惑，別慌張，只需忽略那些花俏的東西。對於中等大小的資料集，一次在記憶體中（例如在列表中）存放過多記錄可能不會是問題。

使用 SEGUID 檢查碼

在這個範例中，我們將 Bio.SeqIO 與 Bio.SeqUtils.CheckSum 模組（在 Biopython 1.44 或更新版本中）一起使用。首先，我們只需印出 GenBank 檔案 ls_orchid.gbk 中每個序列的檢查碼。

from Bio import SeqIO
from Bio.SeqUtils.CheckSum import seguid

for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    print(record.id + "_" + seguid(record.seq))

您應該會得到這樣的輸出：

Z78533.1_JUEoWn6DPhgZ9nAyowsgtoD9TTo
Z78532.1_MN/s0q9zDoCVEEc+k/IFwCNF2pY
...
Z78439.1_H+JfaShya/4yyAj7IbMqgNkxdxQ

現在，讓我們使用檢查碼函數和 Bio.SeqIO.to_dict()，來建立一個使用 SEGUID 作為鍵的 SeqRecord 字典。這裡的技巧是使用 Python lambda 語法來建立一個臨時函數，以取得每個 SeqRecord 的 SEGUID - 我們不能直接使用 seguid() 函數，因為它只適用於 Seq 物件或字串。

from Bio import SeqIO
from Bio.SeqUtils.CheckSum import seguid

seguid_dict = SeqIO.to_dict(
    SeqIO.parse("ls_orchid.gbk", "genbank"), lambda rec: seguid(rec.seq)
)
record = seguid_dict["MN/s0q9zDoCVEEc+k/IFwCNF2pY"]
print(record.id)
print(record.description)

產生這個輸出：

Z78439.1
P.barbatum 5.8S rRNA gene and ITS1 and ITS2 DNA.

隨機子序列

這個腳本將讀取一個包含整個粒線體基因組的 GenBank 檔案（例如，菸草粒線體，Nicotiana tabacum mitochondrion NC_006581），建立 500 筆包含該基因組隨機片段的記錄，並將其儲存為 Fasta 檔案。這些子序列是使用隨機起始點和 200 的固定長度建立的。

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
from random import randint

# There should be one and only one record, the entire genome:
mito_record = SeqIO.read("NC_006581.gbk", "genbank")

mito_frags = []
limit = len(mito_record.seq)
for i in range(0, 500):
    start = randint(0, limit - 200)
    end = start + 200
    mito_frag = mito_record.seq[start:end]
    record = SeqRecord(mito_frag, "fragment_%i" % (i + 1), "", "")
    mito_frags.append(record)

SeqIO.write(mito_frags, "mitofrags.fasta", "fasta")

這應該會產生類似這樣的輸出檔案：

>fragment_1
TGGGCCTCATATTTATCCTATATACCATGTTCGTATGGTGGCGCGATGTTCTACGTGAAT
CCACGTTCGAAGGACATCATACCAAAGTCGTACAATTAGGACCTCGATATGGTTTTATTC
TGTTTATCGTATCGGAGGTTATGTTCTTTTTTGCTCTTTTTCGGGCTTCTTCTCATTCTT
CTTTGGCACCTACGGTAGAG
...
>fragment_500
ACCCAGTGCCGCTACCCACTTCTACTAAGGCTGAGCTTAATAGGAGCAAGAGACTTGGAG
GCAACAACCAGAATGAAATATTATTTAATCGTGGAAATGCCATGTCAGGCGCACCTATCA
GAATCGGAACAGACCAATTACCAGATCCACCTATCATCGCCGGCATAACCATAAAAAAGA
TCATTAAAAAAGCGTGAGCC

寫入字串

有時您不希望將 SeqRecord 物件寫入檔案，而是寫入字串。例如，您可能正在準備要在網頁上顯示的輸出。如果您想將多筆記錄寫入單個字串，請使用 StringIO 來建立基於字串的句柄。教學課程（PDF）在 SeqIO 章節中提供了一個範例。

對於您想要以給定檔案格式將單筆記錄作為字串的特殊情況，Biopython 1.48 新增了一個新的格式方法：

from Bio import SeqIO

mito_record = SeqIO.read("NC_006581.gbk", "genbank")
print(mito_record.format("fasta"))

格式方法將採用 Bio.SeqIO 支援的任何輸出格式，其中檔案格式可用於單筆記錄（例如，"fasta"、"tab" 或 "genbank"）。

請注意，我們不建議您將此用於檔案輸出 - 使用 Bio.SeqIO.write() 更快且更通用。

求助！

如果您在使用 Bio.SeqIO 時遇到問題，請加入討論郵件列表（請參閱郵件列表）。

如果您認為您發現了一個錯誤，請在專案的 GitHub 頁面上回報。