序列輸入/輸出

在本章中，我們將更詳細地討論 Bio.SeqIO 模組，該模組在第快速入門 – Biopython 可以做什麼？章中簡要介紹過，也在第序列註解物件章中使用過。此模組旨在提供一個簡單的介面，以統一的方式處理各種序列檔案格式。另請參閱 Bio.SeqIO wiki 頁面 (https://biopython.dev.org.tw/wiki/SeqIO) 和內建文件 Bio.Seq

>>> from Bio import SeqIO
>>> help(SeqIO)

「重點」是您必須使用 SeqRecord 物件（請參閱第序列註解物件章），其中包含 Seq 物件（請參閱第序列物件章）以及識別碼和描述等註解。請注意，當處理非常大的 FASTA 或 FASTQ 檔案時，使用所有這些物件的額外負擔可能會使腳本過於緩慢。在這種情況下，請考慮使用底層的 SimpleFastaParser 和 FastqGeneralIterator 解析器，它們只會為每個記錄傳回字串的元組（請參閱第底層 FASTA 和 FASTQ 解析器節）。

解析或讀取序列

主要函數 Bio.SeqIO.parse() 用於以 SeqRecord 物件的形式讀取序列資料。此函數需要兩個引數

第一個引數是從中讀取資料的句柄，或檔案名稱。句柄通常是為了讀取而開啟的檔案，但也可能是命令行程式的輸出，或從網際網路下載的資料（請參閱第從網路解析序列節）。有關句柄的詳細資訊，請參閱第什麼是句柄？節。
第二個引數是指定序列格式的小寫字串 – 我們不會嘗試為您猜測檔案格式！有關支援格式的完整列表，請參閱 https://biopython.dev.org.tw/wiki/SeqIO。

Bio.SeqIO.parse() 函數會傳回一個迭代器，該迭代器會產生 SeqRecord 物件。迭代器通常在如下所示的 for 迴圈中使用。

有時您會發現自己處理的檔案只包含單一記錄。對於這種情況，請使用具有相同引數的 Bio.SeqIO.read() 函數。前提是檔案中只有一條記錄，否則會將其傳回為 SeqRecord 物件。否則會引發例外狀況。

讀取序列檔案

一般而言，Bio.SeqIO.parse() 用於以 SeqRecord 物件的形式讀取序列檔案，通常與如下所示的 for 迴圈一起使用

from Bio import SeqIO

for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

上面的範例重複了第解析序列檔案格式節中的簡介，並將會載入 FASTA 格式檔案 ls_orchid.fasta 中的蘭花 DNA 序列。如果改為想要載入像 ls_orchid.gbk 這樣的 GenBank 格式檔案，您只需要變更檔案名稱和格式字串即可

from Bio import SeqIO

for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

同樣地，如果想要讀取其他檔案格式的檔案，假設 Bio.SeqIO.parse() 支援該格式，您只需要適當地變更格式字串即可，例如 SwissProt 檔案的「swiss」或 EMBL 文字檔案的「embl」。在 wiki 頁面 (https://biopython.dev.org.tw/wiki/SeqIO) 和內建文件 Bio.SeqIO 中有一個完整列表。

在清單理解式（或產生器表達式）中使用 Python 迭代器是另一種非常常見的方式。例如，如果您只想從檔案中提取記錄識別碼的列表，我們可以輕鬆地使用下列清單理解式來完成此操作

>>> from Bio import SeqIO
>>> identifiers = [seq_record.id for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]
>>> identifiers
['Z78533.1', 'Z78532.1', 'Z78531.1', 'Z78530.1', 'Z78529.1', 'Z78527.1', ..., 'Z78439.1']

在第序列解析加上簡單的繪圖節中，有更多範例使用了類似於此的清單理解式中的 SeqIO.parse()（例如，用於繪製序列長度或 GC%）。

迭代處理序列檔案中的記錄

在上面的範例中，我們通常使用 for 迴圈逐一迭代處理所有記錄。您可以將 for 迴圈與所有支援迭代介面的各種 Python 物件（包括清單、元組和字串）一起使用。

Bio.SeqIO 傳回的物件實際上是一個迭代器，該迭代器會傳回 SeqRecord 物件。您可以依序查看每個記錄，但只能查看一次。好處是，當處理大型檔案時，迭代器可以節省記憶體。

除了使用 for 迴圈之外，您也可以在迭代器上使用 next() 函數逐步執行項目，如下所示

from Bio import SeqIO

record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")

first_record = next(record_iterator)
print(first_record.id)
print(first_record.description)

second_record = next(record_iterator)
print(second_record.id)
print(second_record.description)

請注意，如果您嘗試使用 next() 並且沒有更多結果，您將會收到特殊的 StopIteration 例外狀況。

要考慮的一種特殊情況是，當您的序列檔案有多個記錄，但您只想要第一個記錄時。在這種情況下，下列程式碼非常簡潔

from Bio import SeqIO

first_record = next(SeqIO.parse("ls_orchid.gbk", "genbank"))

這裡有一個警告 – 像這樣使用 next() 函數會靜默地忽略檔案中的任何其他記錄。如果您的檔案有一個且只有一個記錄，就像本章稍後的一些線上範例或單一染色體的 GenBank 檔案一樣，請改用新的 Bio.SeqIO.read() 函數。這將會檢查是否沒有額外意外的記錄存在。

取得序列檔案中記錄的列表

在上一節中，我們討論了 Bio.SeqIO.parse() 會產生 SeqRecord 迭代器，並且您會逐一取得記錄。通常您需要能夠以任何順序存取記錄。Python list 資料類型非常適合此目的，而且我們可以使用內建的 Python 函數 list() 將記錄迭代器轉換為 SeqRecord 物件的列表，如下所示

from Bio import SeqIO

records = list(SeqIO.parse("ls_orchid.gbk", "genbank"))

print("Found %i records" % len(records))

print("The last record")
last_record = records[-1]  # using Python's list tricks
print(last_record.id)
print(repr(last_record.seq))
print(len(last_record))

print("The first record")
first_record = records[0]  # remember, Python counts from zero
print(first_record.id)
print(repr(first_record.seq))
print(len(first_record))

得到

Found 94 records
The last record
Z78439.1
Seq('CATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCTGTTTACT...GCC')
592
The first record
Z78533.1
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC')
740

您當然仍然可以使用 for 迴圈和 SeqRecord 物件的列表。使用列表比使用迭代器更靈活（例如，您可以從列表的長度判斷記錄的數量），但確實需要更多的記憶體，因為它會一次將所有記錄保存在記憶體中。

提取資料

SeqRecord 物件及其註解結構在第序列註解物件章中有更完整的說明。作為說明如何儲存註解的範例，我們將查看解析 GenBank 檔案 ls_orchid.gbk 中第一個記錄的輸出。

from Bio import SeqIO

record_iterator = SeqIO.parse("ls_orchid.gbk", "genbank")
first_record = next(record_iterator)
print(first_record)

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

ID: Z78533.1
Name: Z78533
Description: C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA.
Number of features: 5
/sequence_version=1
/source=Cypripedium irapeanum
/taxonomy=['Eukaryota', 'Viridiplantae', 'Streptophyta', ..., 'Cypripedium']
/keywords=['5.8S ribosomal RNA', '5.8S rRNA gene', ..., 'ITS1', 'ITS2']
/references=[...]
/accessions=['Z78533']
/data_file_division=PLN
/date=30-NOV-2006
/organism=Cypripedium irapeanum
/gi=2765658
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC')

這段文字提供了 SeqRecord 大部分註釋資料的人類可讀摘要。在這個範例中，我們將使用 .annotations 屬性，它只是一個 Python 字典。當我們在上方列印紀錄時，會顯示這個註釋字典的內容。您也可以直接將它們列印出來。

print(first_record.annotations)

就像任何 Python 字典一樣，您可以輕鬆取得鍵值（keys）

print(first_record.annotations.keys())

或值（values）

print(first_record.annotations.values())

一般而言，註釋值是字串或字串列表。一個特殊情況是，檔案中的任何參考都會儲存為參考物件。

假設您想從 ls_orchid.gbk GenBank 檔案中提取物種列表。我們想要的資訊，Cypripedium irapeanum，保存在註釋字典的 ‘source’ 和 ‘organism’ 下，我們可以像這樣存取：

>>> print(first_record.annotations["source"])
Cypripedium irapeanum

或

>>> print(first_record.annotations["organism"])
Cypripedium irapeanum

一般來說，‘organism’ 用於科學名稱（拉丁文，例如 Arabidopsis thaliana），而 ‘source’ 通常會是俗名（例如，thale cress）。在這個例子中，如同常見的情況，這兩個欄位是相同的。

現在，讓我們瀏覽所有紀錄，建立一個每個蘭花序列來自的物種列表

from Bio import SeqIO

all_species = []
for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    all_species.append(seq_record.annotations["organism"])
print(all_species)

另一種寫這段程式碼的方式是使用列表推導式

from Bio import SeqIO

all_species = [
    seq_record.annotations["organism"]
    for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")
]
print(all_species)

無論哪種方式，結果都是：

['Cypripedium irapeanum', 'Cypripedium californicum', ..., 'Paphiopedilum barbatum']

太棒了。這非常容易，因為 GenBank 檔案是以標準化的方式註釋的。

現在，假設您想從 FASTA 檔案而不是 GenBank 檔案中提取物種列表。壞消息是您必須編寫一些程式碼，從記錄的描述行中提取您想要的資料 — 如果資訊一開始就在檔案中！我們的範例 FASTA 格式檔案 ls_orchid.fasta 開始如下：

>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG
AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG
...

您可以手動檢查，但對於每個紀錄，物種名稱在描述行中是第二個單字。這表示如果我們將每個紀錄的 .description 以空格分隔，那麼物種名稱就在第一個欄位（欄位零是紀錄識別符）。這表示我們可以這樣做：

>>> from Bio import SeqIO
>>> all_species = []
>>> for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
...     all_species.append(seq_record.description.split()[1])
...
>>> print(all_species)  
['C.irapeanum', 'C.californicum', 'C.fasciculatum', ..., 'P.barbatum']

使用列表推導式的簡潔替代方案是：

>>> from Bio import SeqIO
>>> all_species = [
...     seq_record.description.split()[1]
...     for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta")
... ]
>>> print(all_species)  
['C.irapeanum', 'C.californicum', 'C.fasciculatum', ..., 'P.barbatum']

一般來說，從 FASTA 描述行提取資訊並不是很方便。如果您可以將序列以良好註釋的檔案格式（如 GenBank 或 EMBL）取得，那麼這種註釋資訊會更容易處理。

修改資料

在上一節中，我們示範了如何從 SeqRecord 中提取資料。另一個常見的任務是更改此資料。SeqRecord 的屬性可以直接修改，例如：

>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")
>>> first_record = next(record_iterator)
>>> first_record.id
'gi|2765658|emb|Z78533.1|CIZ78533'
>>> first_record.id = "new_id"
>>> first_record.id
'new_id'

請注意，如果您想要更改寫入檔案時 FASTA 的輸出方式（請參閱寫入序列檔案章節），那麼您應該同時修改 id 和 description 屬性。為確保正確的行為，最好在所需的 description 的開頭包含 id 和一個空格。

>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")
>>> first_record = next(record_iterator)
>>> first_record.id = "new_id"
>>> first_record.description = first_record.id + " " + "desired new description"
>>> print(first_record.format("fasta")[:200])
>new_id desired new description
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAA
CGATCGAGTGAATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGT
GACCCTGATTTGTTGTTGGGCCGCCTCGGGAGCGTCCATGGCGGGT

從壓縮檔案剖析序列

在上一節中，我們研究了如何從檔案中剖析序列資料。除了使用檔案名稱，您可以提供 Bio.SeqIO 一個處理（handle）（請參閱什麼是處理（handle）？章節），在本節中，我們將使用處理來剖析壓縮檔案中的序列。

如您在上方所見，我們可以將 Bio.SeqIO.read() 或 Bio.SeqIO.parse() 與檔案名稱一起使用 — 例如，這個快速範例使用生成器表達式計算多紀錄 GenBank 檔案中序列的總長度：

>>> from Bio import SeqIO
>>> print(sum(len(r) for r in SeqIO.parse("ls_orchid.gbk", "gb")))
67518

在這裡，我們改用檔案處理，並使用 with 陳述式自動關閉處理：

>>> from Bio import SeqIO
>>> with open("ls_orchid.gbk") as handle:
...     print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
...
67518

或者，使用老式的方法手動關閉處理：

>>> from Bio import SeqIO
>>> handle = open("ls_orchid.gbk")
>>> print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
67518
>>> handle.close()

現在，假設我們有一個 gzip 壓縮檔案？這些在 Linux 上非常常用。我們可以使用 Python 的 gzip 模組開啟壓縮檔案以進行讀取 — 這會給我們一個處理物件：

>>> import gzip
>>> from Bio import SeqIO
>>> with gzip.open("ls_orchid.gbk.gz", "rt") as handle:
...     print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
...
67518

類似地，如果我們有一個 bzip2 壓縮檔案：

>>> import bz2
>>> from Bio import SeqIO
>>> with bz2.open("ls_orchid.gbk.bz2", "rt") as handle:
...     print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
...
67518

有一種名為 BGZF（區塊 GNU Zip 格式）的 gzip（GNU Zip）變體，可以像普通的 gzip 檔案一樣進行讀取，但稍後在索引壓縮檔案章節中討論的隨機存取方面具有優勢。

從網路剖析序列

在前面的章節中，我們研究了如何從檔案（使用檔案名稱或處理）和從壓縮檔案（使用處理）剖析序列資料。在這裡，我們將使用 Bio.SeqIO 和另一種處理，網路連線，從網際網路下載和剖析序列。

請注意，僅僅因為您可以一次性下載序列資料並將其剖析為 SeqRecord 物件並不表示這是一個好主意。一般來說，您應該一次下載序列並將其儲存到檔案中以供重複使用。

從網路剖析 GenBank 記錄

EFetch：從 Entrez 下載完整記錄章節更詳細地討論了 Entrez EFetch 介面，但現在讓我們連接到 NCBI 並使用其 GI 號碼從 GenBank 取得一些 Opuntia（仙人掌）序列。

首先，讓我們只提取一筆記錄。如果您不關心註釋和特徵，下載 FASTA 檔案是一個不錯的選擇，因為它們很精簡。現在請記住，當您預期處理只包含一筆記錄時，請使用 Bio.SeqIO.read() 函式：

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"
with Entrez.efetch(
    db="nucleotide", rettype="fasta", retmode="text", id="6273291"
) as handle:
    seq_record = SeqIO.read(handle, "fasta")
print("%s with %i features" % (seq_record.id, len(seq_record.features)))

預期輸出：

gi|6273291|gb|AF191665.1|AF191665 with 0 features

NCBI 也允許您要求以其他格式提供檔案，特別是以 GenBank 檔案的形式。在 2009 年復活節之前，Entrez EFetch API 允許您使用 “genbank” 作為返回類型，但 NCBI 現在堅持使用 “gb”（或蛋白質的 “gp”）的官方返回類型，如 EFetch for Sequence and other Molecular Biology Databases 中所述。因此，在 Biopython 1.50 及更高版本中，我們在 Bio.SeqIO 中支援 “gb” 作為 “genbank” 的別名。

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"
with Entrez.efetch(
    db="nucleotide", rettype="gb", retmode="text", id="6273291"
) as handle:
    seq_record = SeqIO.read(handle, "gb")  # using "gb" as an alias for "genbank"
print("%s with %i features" % (seq_record.id, len(seq_record.features)))

這個範例的預期輸出是：

AF191665.1 with 3 features

請注意，這次我們有三個特徵。

現在讓我們提取多筆記錄。這次處理包含多筆記錄，因此我們必須使用 Bio.SeqIO.parse() 函式：

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"
with Entrez.efetch(
    db="nucleotide", rettype="gb", retmode="text", id="6273291,6273290,6273289"
) as handle:
    for seq_record in SeqIO.parse(handle, "gb"):
        print("%s %s..." % (seq_record.id, seq_record.description[:50]))
        print(
            "Sequence length %i, %i features, from: %s"
            % (
                len(seq_record),
                len(seq_record.features),
                seq_record.annotations["source"],
            )
        )

應該會得到以下輸出：

AF191665.1 Opuntia marenae rpl16 gene; chloroplast gene for c...
Sequence length 902, 3 features, from: chloroplast Opuntia marenae
AF191664.1 Opuntia clavata rpl16 gene; chloroplast gene for c...
Sequence length 899, 3 features, from: chloroplast Grusonia clavata
AF191663.1 Opuntia bradtiana rpl16 gene; chloroplast gene for...
Sequence length 899, 3 features, from: chloroplast Opuntia bradtianaa

有關 Bio.Entrez 模組的更多資訊，請參閱存取 NCBI 的 Entrez 資料庫章節，並務必閱讀有關使用 Entrez 的 NCBI 指南（Entrez 指南章節）。

從網路剖析 SwissProt 序列

現在，讓我們使用處理從 ExPASy 下載 SwissProt 檔案，這在 Swiss-Prot 和 ExPASy 章節中更深入地介紹。如上所述，當您預期處理只包含一筆記錄時，請使用 Bio.SeqIO.read() 函式：

from Bio import ExPASy
from Bio import SeqIO

with ExPASy.get_sprot_raw("O23729") as handle:
    seq_record = SeqIO.read(handle, "swiss")
print(seq_record.id)
print(seq_record.name)
print(seq_record.description)
print(repr(seq_record.seq))
print("Length %i" % len(seq_record))
print(seq_record.annotations["keywords"])

假設您的網路連線正常，您應該會收到：

O23729
CHS3_BROFI
RecName: Full=Chalcone synthase 3; EC=2.3.1.74; AltName: Full=Naringenin-chalcone synthase 3;
Seq('MAPAMEEIRQAQRAEGPAAVLAIGTSTPPNALYQADYPDYYFRITKSEHLTELK...GAE')
Length 394
['Acyltransferase', 'Flavonoid biosynthesis', 'Transferase']

將序列檔案視為字典

對 SeqIO.parse 返回的迭代器進行一次迴圈會耗盡檔案。對於自我索引的檔案（例如 twoBit 格式的檔案），SeqIO.parse 的返回值也可以用作字典，允許隨機存取序列內容。由於在這種情況下是按需進行剖析，因此只要存取序列資料，檔案就必須保持開啟狀態。

>>> from Bio import SeqIO
>>> handle = open("sequence.bigendian.2bit", "rb")
>>> records = SeqIO.parse(handle, "twobit")
>>> records.keys()
dict_keys(['seq11111', 'seq222', 'seq3333', 'seq4', 'seq555', 'seq6'])
>>> records["seq222"]
SeqRecord(seq=Seq('TTGATCGGTGACAAATTTTTTACAAAGAACTGTAGGACTTGCTACTTCTCCCTC...ACA'), id='seq222', name='<unknown name>', description='<unknown description>', dbxrefs=[])
>>> records["seq222"].seq
Seq('TTGATCGGTGACAAATTTTTTACAAAGAACTGTAGGACTTGCTACTTCTCCCTC...ACA')
>>> handle.close()
>>> records["seq222"].seq
Traceback (most recent call last):
...
ValueError: cannot retrieve sequence: file is closed

對於其他檔案格式，Bio.SeqIO 提供了三個相關的函式模組，允許對多序列檔案進行類似字典的隨機存取。在這裡，靈活性和記憶體使用之間需要權衡。總之：

Bio.SeqIO.to_dict() 是最靈活但也是記憶體需求最高的選項（請參閱將序列檔案視為字典 – 記憶體中章節）。這基本上是一個輔助函式，用於建構一個普通的 Python dictionary，其中每個條目都以 SeqRecord 物件的形式保存在記憶體中，允許您修改記錄。
Bio.SeqIO.index() 是一個有用的中間方案，其作用類似於唯讀字典，並根據需求將序列解析為 SeqRecord 物件（請參閱將序列檔案作為字典 – 索引檔案章節）。
Bio.SeqIO.index_db() 的作用也類似於唯讀字典，但它將識別符和檔案偏移量儲存在磁碟上的檔案中（作為 SQLite3 資料庫），這表示它具有非常低的記憶體需求（請參閱將序列檔案作為字典 – 資料庫索引檔案章節），但速度會稍微慢一些。

請參閱概述的討論（討論章節）。

將序列檔案作為字典 – 記憶體中

接下來，我們將使用我們無處不在的蘭花檔案，展示如何對它們進行索引，並像使用 Python dictionary 資料類型（類似 Perl 中的 hash）一樣存取它們。這對於只需要存取檔案中某些元素的適中大小的檔案非常有用，並且可以建立一個快速且簡單的資料庫。對於處理記憶體成為問題的較大檔案，請參閱以下將序列檔案作為字典 – 索引檔案章節。

您可以使用函數 Bio.SeqIO.to_dict() 來建立一個 SeqRecord 字典（在記憶體中）。預設情況下，這會使用每個記錄的識別符（即 .id 屬性）作為鍵。讓我們嘗試使用我們的 GenBank 檔案。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.gbk", "genbank"))

Bio.SeqIO.to_dict() 只有一個必需的參數，即提供 SeqRecord 物件的列表或產生器。在這裡，我們僅使用了 SeqIO.parse 函數的輸出。顧名思義，這會傳回一個 Python 字典。

由於此變數 orchid_dict 是一個普通的 Python 字典，我們可以查看所有可用的鍵。

>>> len(orchid_dict)
94

>>> list(orchid_dict.keys())
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

在 Python 3 下，字典方法（例如「.keys()」和「.values()」）是迭代器，而不是列表。

如果您真的想這樣做，您甚至可以一次查看所有記錄。

>>> list(orchid_dict.values())  # lots of output!

我們可以通過鍵存取單個 SeqRecord 物件，並像平常一樣操作該物件。

>>> seq_record = orchid_dict["Z78475.1"]
>>> print(seq_record.description)
P.supardii 5.8S rRNA gene and ITS1 and ITS2 DNA
>>> seq_record.seq
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...GGT')

因此，建立我們 GenBank 記錄的記憶體中「資料庫」非常容易。接下來，我們將嘗試對 FASTA 檔案執行此操作。

請注意，有先前 Python 經驗的您應該都能「手動」建構這樣的字典。但是，典型的字典建構方法不會很好地處理重複鍵的情況。使用 Bio.SeqIO.to_dict() 會明確檢查是否有重複鍵，如果找到任何重複鍵，則會引發例外。

指定字典鍵

使用與上述相同的程式碼，但改為使用 FASTA 檔案。

from Bio import SeqIO

orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.fasta", "fasta"))
print(orchid_dict.keys())

這次的鍵是

['gi|2765596|emb|Z78471.1|PDZ78471', 'gi|2765646|emb|Z78521.1|CCZ78521', ...
 ..., 'gi|2765613|emb|Z78488.1|PTZ78488', 'gi|2765583|emb|Z78458.1|PHZ78458']

您應該能從我們在簡單 FASTA 解析範例章節中較早解析 FASTA 檔案時識別出這些字串。假設您希望將其他內容作為鍵，例如登錄號。這將我們引導至 SeqIO.to_dict() 的可選參數 key_function，該參數可讓您定義要將哪些內容用作記錄的字典鍵。

首先，您必須編寫自己的函數，以在給定 SeqRecord 物件時傳回您想要的鍵（作為字串）。一般而言，函數的詳細資訊將取決於您正在處理的輸入記錄類型。但是，對於我們的蘭花，我們可以使用「管線」字元（垂直線）分割記錄的識別符，然後傳回第四個條目（第三個欄位）。

def get_accession(record):
    """Given a SeqRecord, return the accession number as a string.

    e.g. "gi|2765613|emb|Z78488.1|PTZ78488" -> "Z78488.1"
    """
    parts = record.id.split("|")
    assert len(parts) == 5 and parts[0] == "gi" and parts[2] == "emb"
    return parts[3]

然後，我們可以將此函數提供給 SeqIO.to_dict() 函數，以在建立字典時使用。

from Bio import SeqIO

orchid_dict = SeqIO.to_dict(
    SeqIO.parse("ls_orchid.fasta", "fasta"), key_function=get_accession
)
print(orchid_dict.keys())

最後，正如所希望的，新的字典鍵。

>>> print(orchid_dict.keys())
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

希望這不會太複雜！

使用 SEGUID 檢查碼為字典建立索引

為了提供另一個使用 SeqRecord 物件字典的範例，我們將使用 SEGUID 檢查碼函數。這是一個相對較新的檢查碼，碰撞應該非常罕見（即兩個不同的序列具有相同的檢查碼），這是對 CRC64 檢查碼的改進。

再次使用蘭花 GenBank 檔案。

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() 函數的 key_function 參數，該參數預期一個將 SeqRecord 轉換為字串的函數。我們無法直接使用 seguid() 函數，因為它預期會被給予一個 Seq 物件（或字串）。但是，我們可以使用 Python 的 lambda 功能來建立一個「一次性」函數，以傳遞給 Bio.SeqIO.to_dict()。

>>> 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)
Z78532.1
>>> print(record.description)
C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA

這應該會檢索檔案中的第二個條目 Z78532.1 記錄。

將序列檔案作為字典 – 索引檔案

正如先前的幾個範例試圖說明的那樣，使用 Bio.SeqIO.to_dict() 非常靈活。但是，由於它將所有內容都保留在記憶體中，因此您可以使用的檔案大小會受到電腦 RAM 的限制。一般而言，這僅適用於小型到中型的檔案。

對於較大的檔案，您應該考慮使用 Bio.SeqIO.index()，其工作方式略有不同。儘管它仍然會傳回類似字典的物件，但這並不會將所有內容都保留在記憶體中。相反，它只會記錄每個記錄在檔案中的位置 – 當您要求特定記錄時，它會根據需求解析該記錄。

舉例來說，我們使用與之前相同的 GenBank 檔案。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.gbk", "genbank")
>>> len(orchid_dict)
94

>>> orchid_dict.keys()
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

>>> seq_record = orchid_dict["Z78475.1"]
>>> print(seq_record.description)
P.supardii 5.8S rRNA gene and ITS1 and ITS2 DNA
>>> seq_record.seq
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...GGT')
>>> orchid_dict.close()

請注意，Bio.SeqIO.index() 不會接受控制代碼，而只會接受檔案名稱。這樣做的原因很充分，但有點技術性。第二個參數是檔案格式（一個小寫字串，與其他 Bio.SeqIO 函數中使用的一樣）。您可以使用許多其他簡單的檔案格式，包括 FASTA 和 FASTQ 檔案（請參閱為 FASTQ 檔案建立索引章節中的範例）。但是，不支援 PHYLIP 或 Clustal 等對齊格式。最後，作為一個可選參數，您可以提供一個鍵函數。

以下是使用 FASTA 檔案的相同範例 – 我們所做的只是變更檔案名稱和格式名稱。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.fasta", "fasta")
>>> len(orchid_dict)
94
>>> orchid_dict.keys()
['gi|2765596|emb|Z78471.1|PDZ78471', 'gi|2765646|emb|Z78521.1|CCZ78521', ...
 ..., 'gi|2765613|emb|Z78488.1|PTZ78488', 'gi|2765583|emb|Z78458.1|PHZ78458']

指定字典鍵

假設您想要使用與之前相同的鍵？與指定字典鍵章節中的 Bio.SeqIO.to_dict() 範例非常相似，您需要編寫一個微小的函數，以將 FASTA 識別符（作為字串）對應到您想要的鍵。

def get_acc(identifier):
    """Given a SeqRecord identifier string, return the accession number as a string.

    e.g. "gi|2765613|emb|Z78488.1|PTZ78488" -> "Z78488.1"
    """
    parts = identifier.split("|")
    assert len(parts) == 5 and parts[0] == "gi" and parts[2] == "emb"
    return parts[3]

然後，我們可以將此函數提供給 Bio.SeqIO.index() 函數，以在建立字典時使用。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.fasta", "fasta", key_function=get_acc)
>>> print(orchid_dict.keys())
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

當您知道如何操作時，就會變得容易？

取得記錄的原始資料

來自 Bio.SeqIO.index() 的類字典物件會將每個條目當作 SeqRecord 物件提供給您。但是，有時能夠直接從檔案中取得原始資料非常有用。為此，請使用 get_raw() 方法，該方法接受單個引數（記錄識別符）並傳回位元組字串（從檔案中提取且未經修改）。

一個令人信服的範例是從大型檔案中提取記錄子集，其中 Bio.SeqIO.write() 尚不支援輸出檔案格式（例如，純文字 SwissProt 檔案格式），或者您需要完全保留文字（例如，來自 Biopython 的 GenBank 或 EMBL 輸出尚未保留每一位註解）。

假設您已從其 FTP 站點 (ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.dat.gz) 以下載整個 UniProt 的純文字 SwissPort 檔案格式，並將其解壓縮為 uniprot_sprot.dat 檔案，而您只想從中提取一些記錄。

>>> from Bio import SeqIO
>>> uniprot = SeqIO.index("uniprot_sprot.dat", "swiss")
>>> with open("selected.dat", "wb") as out_handle:
...     for acc in ["P33487", "P19801", "P13689", "Q8JZQ5", "Q9TRC7"]:
...         out_handle.write(uniprot.get_raw(acc))
...

請注意，從 Python 3 開始，我們必須以二進制模式開啟檔案進行寫入，因為 get_raw() 方法會傳回位元組字串。

在排序序列檔案章節中有一個較長的範例，其中使用 SeqIO.index() 函數對大型序列檔案進行排序（而無需一次將所有內容載入記憶體中）。

將序列檔案作為字典 – 資料庫索引檔案

Biopython 1.57 引入了替代方案 Bio.SeqIO.index_db()，由於它將記錄資訊儲存在磁碟上的檔案中（使用 SQLite3 資料庫），而不是在記憶體中，因此可以處理非常大的檔案。此外，您可以一起為多個檔案建立索引（前提是所有記錄識別符都是唯一的）。

Bio.SeqIO.index() 函數接受三個必需的參數。

索引檔案名稱，我們建議使用以 .idx 結尾的名稱。此索引檔案實際上是一個 SQLite3 資料庫。
要建立索引的序列檔案清單（或單個檔案名稱）。
檔案格式（與 SeqIO 模組的其餘部分使用相同的小寫字串）。

舉例來說，請考慮 NCBI FTP 站點 (ftp://ftp.ncbi.nih.gov/genbank/) 中的 GenBank 平面檔案版本，它們是經過 gzip 壓縮的 GenBank 檔案。

截至 GenBank 版本 \(210\)，共有 \(38\) 個檔案組成了病毒序列，即 gbvrl1.seq, ..., gbvrl38.seq，解壓縮後佔用磁碟上約 8GB 的空間，總共包含近 200 萬條記錄。

如果您對病毒感興趣，可以使用 rsync 命令從命令列輕鬆下載所有病毒檔案，然後使用 gunzip 解壓縮它們。

# For illustration only, see reduced example below
$ rsync -avP "ftp.ncbi.nih.gov::genbank/gbvrl*.seq.gz" .
$ gunzip gbvrl*.seq.gz

除非您很在意病毒，否則為了這個範例下載這麼多資料實在太多了 - 因此，我們只下載前四個區塊（每個壓縮後約 25MB），並將它們解壓縮（總共佔用約 1GB 的空間）。

# Reduced example, download only the first four chunks
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl1.seq.gz
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl2.seq.gz
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl3.seq.gz
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl4.seq.gz
$ gunzip gbvrl*.seq.gz

現在，在 Python 中，將這些 GenBank 檔案編入索引，如下所示：

>>> import glob
>>> from Bio import SeqIO
>>> files = glob.glob("gbvrl*.seq")
>>> print("%i files to index" % len(files))
4
>>> gb_vrl = SeqIO.index_db("gbvrl.idx", files, "genbank")
>>> print("%i sequences indexed" % len(gb_vrl))
272960 sequences indexed

在我的機器上，為完整的病毒 GenBank 檔案集建立索引大約花了十分鐘，而前四個檔案則約花了一分鐘左右。

但是，一旦完成，重複執行此操作將會在極短的時間內重新載入索引檔案 gbvrl.idx。

您可以將索引作為唯讀的 Python 字典使用，而不必擔心序列來自哪個檔案，例如：

>>> print(gb_vrl["AB811634.1"].description)
Equine encephalosis virus NS3 gene, complete cds, isolate: Kimron1.

取得記錄的原始資料

就像上面在取得記錄的原始資料章節中討論的 Bio.SeqIO.index() 函式一樣，這個類似字典的物件也讓您可以取得每個記錄的原始位元組。

>>> print(gb_vrl.get_raw("AB811634.1"))
LOCUS       AB811634                 723 bp    RNA     linear   VRL 17-JUN-2015
DEFINITION  Equine encephalosis virus NS3 gene, complete cds, isolate: Kimron1.
ACCESSION   AB811634
...
//

為壓縮檔建立索引

通常在為序列檔案建立索引時，檔案可能會非常大，因此您可能會想將其壓縮在磁碟上。不幸的是，對於較常見的檔案格式（如 gzip 和 bzip2），高效的隨機存取很困難。在這種情況下，BGZF（Blocked GNU Zip Format）會很有幫助。這是 gzip 的一種變體（可以使用標準 gzip 工具解壓縮），因 BAM 檔案格式、samtools 和 tabix 而廣受歡迎。

若要建立 BGZF 壓縮檔，您可以使用隨附於 samtools 的命令列工具 bgzip。在我們的範例中，我們使用副檔名 *.bgz，以便將它們與一般的 gzipped 檔案（命名為 *.gz）區分開來。您也可以使用 Bio.bgzf 模組在 Python 中讀取和寫入 BGZF 檔案。

Bio.SeqIO.index() 和 Bio.SeqIO.index_db() 都可以與 BGZF 壓縮檔一起使用。例如，如果您從未壓縮的 GenBank 檔案開始：

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.gbk", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.close()

您可以使用以下命令在命令列中壓縮此檔案（同時保留原始檔案）- 但不用擔心，壓縮檔案已經包含在其他範例檔案中。

$ bgzip -c ls_orchid.gbk > ls_orchid.gbk.bgz

您可以使用與之前完全相同的方式使用壓縮檔案。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.gbk.bgz", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.close()

或

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index_db("ls_orchid.gbk.bgz.idx", "ls_orchid.gbk.bgz", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.close()

SeqIO 索引會自動偵測 BGZF 壓縮。請注意，您不能為未壓縮和壓縮的檔案使用相同的索引檔案。

討論

那麼，您應該使用這些方法中的哪一種，為什麼呢？這取決於您嘗試做的事情（以及您正在處理的資料量）。但是，一般來說，選擇 Bio.SeqIO.index() 是一個不錯的起點。如果您正在處理數百萬條記錄、多個檔案或重複的分析，那麼請考慮使用 Bio.SeqIO.index_db()。

選擇 Bio.SeqIO.to_dict() 而不是 Bio.SeqIO.index() 或 Bio.SeqIO.index_db() 的原因，歸根結底是儘管其記憶體需求很高，但仍需要彈性。將 SeqRecord 物件儲存在記憶體中的優點是，可以隨意更改、新增或移除它們。除了消耗大量記憶體的缺點之外，建立索引也可能需要更長的時間，因為必須完整剖析所有記錄。

Bio.SeqIO.index() 和 Bio.SeqIO.index_db() 都只在需要時剖析記錄。在建立索引時，它們會掃描檔案一次，尋找每個記錄的開頭，並盡可能少地工作來擷取識別碼。

選擇 Bio.SeqIO.index() 而不是 Bio.SeqIO.index_db() 的原因包括：

建立索引的速度更快（在簡單的檔案格式中更明顯）。
以 SeqRecord 物件存取的速度略快（但差異只有在剖析簡單的檔案格式時才真正明顯）。
可以使用任何不可變的 Python 物件作為字典鍵（例如字串的元組或凍結的集合），而不僅僅是字串。
如果正在編製索引的序列檔案已變更，則無需擔心索引資料庫過時。

選擇 Bio.SeqIO.index_db() 而不是 Bio.SeqIO.index() 的原因包括：

不受記憶體限制 - 這對於來自第二代定序的檔案已經很重要，其中數千萬個序列很常見，而使用 Bio.SeqIO.index() 可能需要超過 4GB 的 RAM，因此需要 64 位元的 Python 版本。
由於索引會保留在磁碟上，因此可以重複使用。雖然建立索引資料庫檔案需要更長的時間，但如果您有腳本會在未來對相同的資料檔案重新執行，從長遠來看，這可以節省時間。
一起為多個檔案建立索引
get_raw() 方法可能會快很多，因為對於大多數檔案格式，每個記錄的長度和偏移量都會被儲存。

寫入序列檔案

我們已經討論過使用 Bio.SeqIO.parse() 來進行序列輸入（讀取檔案），現在我們將看看 Bio.SeqIO.write()，它是用於序列輸出（寫入檔案）。這是一個採用三個參數的函式：一些 SeqRecord 物件、要寫入的控制代碼或檔案名稱，以及序列格式。

這裡有一個範例，我們首先以比較困難的方式（手動，而不是從檔案載入）建立一些 SeqRecord 物件：

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

rec1 = SeqRecord(
    Seq(
        "MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD"
        "GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK"
        "NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM"
        "SSAC",
    ),
    id="gi|14150838|gb|AAK54648.1|AF376133_1",
    description="chalcone synthase [Cucumis sativus]",
)

rec2 = SeqRecord(
    Seq(
        "YPDYYFRITNREHKAELKEKFQRMCDKSMIKKRYMYLTEEILKENPSMCEYMAPSLDARQ"
        "DMVVVEIPKLGKEAAVKAIKEWGQ",
    ),
    id="gi|13919613|gb|AAK33142.1|",
    description="chalcone synthase [Fragaria vesca subsp. bracteata]",
)

rec3 = SeqRecord(
    Seq(
        "MVTVEEFRRAQCAEGPATVMAIGTATPSNCVDQSTYPDYYFRITNSEHKVELKEKFKRMC"
        "EKSMIKKRYMHLTEEILKENPNICAYMAPSLDARQDIVVVEVPKLGKEAAQKAIKEWGQP"
        "KSKITHLVFCTTSGVDMPGCDYQLTKLLGLRPSVKRFMMYQQGCFAGGTVLRMAKDLAEN"
        "NKGARVLVVCSEITAVTFRGPNDTHLDSLVGQALFGDGAAAVIIGSDPIPEVERPLFELV"
        "SAAQTLLPDSEGAIDGHLREVGLTFHLLKDVPGLISKNIEKSLVEAFQPLGISDWNSLFW"
        "IAHPGGPAILDQVELKLGLKQEKLKATRKVLSNYGNMSSACVLFILDEMRKASAKEGLGT"
        "TGEGLEWGVLFGFGPGLTVETVVLHSVAT",
    ),
    id="gi|13925890|gb|AAK49457.1|",
    description="chalcone synthase [Nicotiana tabacum]",
)

my_records = [rec1, rec2, rec3]

現在我們有一個 SeqRecord 物件列表，我們將它們寫入 FASTA 格式的檔案：

from Bio import SeqIO

SeqIO.write(my_records, "my_example.faa", "fasta")

如果您在您最愛的文字編輯器中開啟此檔案，它應該看起來像這樣：

>gi|14150838|gb|AAK54648.1|AF376133_1 chalcone synthase [Cucumis sativus]
MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD
GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK
NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM
SSAC
>gi|13919613|gb|AAK33142.1| chalcone synthase [Fragaria vesca subsp. bracteata]
YPDYYFRITNREHKAELKEKFQRMCDKSMIKKRYMYLTEEILKENPSMCEYMAPSLDARQ
DMVVVEIPKLGKEAAVKAIKEWGQ
>gi|13925890|gb|AAK49457.1| chalcone synthase [Nicotiana tabacum]
MVTVEEFRRAQCAEGPATVMAIGTATPSNCVDQSTYPDYYFRITNSEHKVELKEKFKRMC
EKSMIKKRYMHLTEEILKENPNICAYMAPSLDARQDIVVVEVPKLGKEAAQKAIKEWGQP
KSKITHLVFCTTSGVDMPGCDYQLTKLLGLRPSVKRFMMYQQGCFAGGTVLRMAKDLAEN
NKGARVLVVCSEITAVTFRGPNDTHLDSLVGQALFGDGAAAVIIGSDPIPEVERPLFELV
SAAQTLLPDSEGAIDGHLREVGLTFHLLKDVPGLISKNIEKSLVEAFQPLGISDWNSLFW
IAHPGGPAILDQVELKLGLKQEKLKATRKVLSNYGNMSSACVLFILDEMRKASAKEGLGT
TGEGLEWGVLFGFGPGLTVETVVLHSVAT

假設您想知道 Bio.SeqIO.write() 函式寫入控制代碼的記錄數量？如果您的記錄在列表中，您可以直接使用 len(my_records)，但是當您的記錄來自產生器/迭代器時，您無法這麼做。Bio.SeqIO.write() 函式會傳回寫入檔案的 SeqRecord 物件數量。

注意 - 如果您告訴 Bio.SeqIO.write() 函式寫入已經存在的檔案，則會覆寫舊檔案，而不會發出任何警告。

往返

有些人希望他們的剖析器是「可往返的」，這意味著如果您讀取檔案並將其寫回，它會保持不變。這需要剖析器擷取足夠的資訊才能完全重現原始檔案。Bio.SeqIO 的目標不是這樣做。

舉一個簡單的例子，允許 FASTA 檔案中序列資料的任何換行。從剖析以下兩個僅行分隔符號不同的範例中，會給出相同的 SeqRecord：

>YAL068C-7235.2170 Putative promoter sequence
TACGAGAATAATTTCTCATCATCCAGCTTTAACACAAAATTCGCACAGTTTTCGTTAAGA
GAACTTAACATTTTCTTATGACGTAAATGAAGTTTATATATAAATTTCCTTTTTATTGGA

>YAL068C-7235.2170 Putative promoter sequence
TACGAGAATAATTTCTCATCATCCAGCTTTAACACAAAATTCGCA
CAGTTTTCGTTAAGAGAACTTAACATTTTCTTATGACGTAAATGA
AGTTTATATATAAATTTCCTTTTTATTGGA

若要建立可往返的 FASTA 剖析器，您需要追蹤序列換行的位置，而這些額外資訊通常毫無意義。相反地，Biopython 在輸出上使用預設的 \(60\) 個字元的換行。空白字元的相同問題也適用於許多其他檔案格式。在某些情況下，另一個問題是 Biopython 尚未保留每最後一點註釋（例如 GenBank 和 EMBL）。

有時，保留原始佈局（及其可能有的任何怪癖）非常重要。請參閱關於 Bio.SeqIO.index() 類似字典物件的 get_raw() 方法的取得記錄的原始資料章節，以取得一種可能的解決方案。

在序列檔案格式之間轉換

在先前的範例中，我們使用 SeqRecord 物件列表作為 Bio.SeqIO.write() 函式的輸入，但它也會接受像我們從 Bio.SeqIO.parse() 取得的 SeqRecord 迭代器 - 這讓我們可以透過組合這兩個函式來進行檔案轉換。

在這個範例中，我們將讀取 GenBank 格式的檔案 ls_orchid.gbk 並以 FASTA 格式寫出：

from Bio import SeqIO

records = SeqIO.parse("ls_orchid.gbk", "genbank")
count = SeqIO.write(records, "my_example.fasta", "fasta")
print("Converted %i records" % count)

不過，這有點複雜。因此，由於檔案轉換是一項很常見的工作，因此有一個輔助函式可讓您用下列內容替換它：

from Bio import SeqIO

count = SeqIO.convert("ls_orchid.gbk", "genbank", "my_example.fasta", "fasta")
print("Converted %i records" % count)

Bio.SeqIO.convert() 函式會接受控制代碼或檔案名稱。請注意 - 如果輸出檔案已經存在，它將會覆寫它！若要了解更多資訊，請參閱內建的說明：

>>> from Bio import SeqIO
>>> help(SeqIO.convert)

原則上，只要變更檔案名稱和格式名稱，此程式碼即可用於在 Biopython 中提供的任何檔案格式之間進行轉換。但是，寫入某些格式需要其他檔案格式不包含的資訊（例如品質分數）。例如，雖然您可以將 FASTQ 檔案轉換為 FASTA 檔案，但您無法反向執行此操作。另請參閱章節轉換 FASTQ 檔案和將 FASTA 和 QUAL 檔案轉換為 FASTQ 檔案，這些章節會探討不同 FASTQ 格式之間的相互轉換。

最後，作為使用 Bio.SeqIO.convert() 函式的一個額外誘因（除了程式碼會更簡短之外），以這種方式執行也可能會更快！原因在於轉換函式可以利用幾個檔案格式特定的最佳化和技巧。

將序列檔案轉換為其反向互補序列

假設您有一個核苷酸序列檔案，並且想要將其轉換為包含其反向互補序列的檔案。這次需要一些工作才能將我們從輸入檔案中取得的 SeqRecord 物件轉換成適合儲存到輸出檔案中的格式。

首先，我們將使用 Bio.SeqIO.parse() 從檔案載入一些核苷酸序列，然後使用 Seq 物件內建的 .reverse_complement() 方法列印出它們的反向互補序列（請參閱核苷酸序列與（反向）互補序列章節）。

>>> from Bio import SeqIO
>>> for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
...     print(record.id)
...     print(record.seq.reverse_complement())
...

現在，如果我們想將這些反向互補序列儲存到檔案中，我們需要建立 SeqRecord 物件。我們可以利用 SeqRecord 物件內建的 .reverse_complement() 方法（請參閱反向互補 SeqRecord 物件章節），但我們必須決定如何命名我們的新紀錄。

這是一個很好的地方來展示列表推導式的強大之處，它可以在記憶體中建立一個列表

>>> from Bio import SeqIO
>>> records = [
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
... ]
>>> len(records)
94

現在，列表推導式有一個很棒的小技巧，您可以新增一個條件語句

>>> records = [
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
...     if len(rec) < 700
... ]
>>> len(records)
18

這將建立一個反向互補記錄的記憶體列表，其中序列長度小於 700 個鹼基對。然而，我們可以使用產生器表達式來做完全相同的事情 - 但優點是它不會一次在記憶體中建立所有記錄的列表

>>> records = (
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
...     if len(rec) < 700
... )

作為一個完整的例子

>>> from Bio import SeqIO
>>> records = (
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
...     if len(rec) < 700
... )
>>> SeqIO.write(records, "rev_comp.fasta", "fasta")
18

在翻譯 FASTA 檔案的 CDS 條目章節中有一個相關的例子，將 FASTA 檔案中的每個記錄從核苷酸翻譯成胺基酸。

將您的 SeqRecord 物件轉換為格式化的字串

假設您並不是真的想將您的記錄寫入檔案或處理 - 而是想要一個字串，其中包含特定檔案格式的記錄。Bio.SeqIO 介面基於句柄，但 Python 有一個有用的內建模組，可提供基於字串的句柄。

舉例說明您可能會如何使用它，讓我們從我們的蘭花 GenBank 檔案中載入一堆 SeqRecord 物件，並建立一個包含 FASTA 格式記錄的字串

from Bio import SeqIO
from io import StringIO

records = SeqIO.parse("ls_orchid.gbk", "genbank")
out_handle = StringIO()
SeqIO.write(records, out_handle, "fasta")
fasta_data = out_handle.getvalue()
print(fasta_data)

第一次看到時，這並非完全直觀！好的一面是，對於您希望字串包含特定檔案格式的單個記錄的特殊情況，請使用 SeqRecord 類別的 format() 方法（請參閱format 方法章節）。

請注意，儘管我們不鼓勵這樣做，但您可以使用 format() 方法寫入檔案，例如像這樣

from Bio import SeqIO

with open("ls_orchid_long.tab", "w") as out_handle:
    for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
        if len(record) > 100:
            out_handle.write(record.format("tab"))

雖然這種程式碼風格適用於簡單的循序檔案格式，例如 FASTA 或這裡使用的簡單的 Tab 分隔格式，但它不適用於更複雜或交錯的檔案格式。這就是為什麼我們仍然建議使用 Bio.SeqIO.write()，如下例所示

from Bio import SeqIO

records = (rec for rec in SeqIO.parse("ls_orchid.gbk", "genbank") if len(rec) > 100)
SeqIO.write(records, "ls_orchid.tab", "tab")

單次呼叫 SeqIO.write(...) 也比多次呼叫 SeqRecord.format(...) 方法快得多。

底層 FASTA 和 FASTQ 解析器

當處理大型高通量 FASTA 或 FASTQ 定序檔案，且速度至關重要時，使用底層的 SimpleFastaParser 或 FastqGeneralIterator 通常比 Bio.SeqIO.parse 更實用。如本章引言中所述，檔案格式中性的 Bio.SeqIO 介面即使對於像 FASTA 這樣的簡單格式，也具有建立許多物件的開銷。

在解析 FASTA 檔案時，內部 Bio.SeqIO.parse() 會呼叫帶有檔案句柄的底層 SimpleFastaParser。您可以直接使用它 - 它會遍歷檔案句柄，將每個記錄作為兩個字串的元組傳回，標題行（> 字元之後的所有內容）和序列（作為純字串）

>>> from Bio.SeqIO.FastaIO import SimpleFastaParser
>>> count = 0
>>> total_len = 0
>>> with open("ls_orchid.fasta") as in_handle:
...     for title, seq in SimpleFastaParser(in_handle):
...         count += 1
...         total_len += len(seq)
...
>>> print("%i records with total sequence length %i" % (count, total_len))
94 records with total sequence length 67518

只要您不關心換行（對於短讀高通量資料來說，您可能不需要），從這些字串輸出 FASTA 格式也非常快

...
out_handle.write(">%s\n%s\n" % (title, seq))
...

同樣地，在解析 FASTQ 檔案時，內部 Bio.SeqIO.parse() 會呼叫帶有檔案句柄的底層 FastqGeneralIterator。如果您不需要將品質分數轉換為整數，或者可以使用它們作為 ASCII 字串，這就很理想

>>> from Bio.SeqIO.QualityIO import FastqGeneralIterator
>>> count = 0
>>> total_len = 0
>>> with open("example.fastq") as in_handle:
...     for title, seq, qual in FastqGeneralIterator(in_handle):
...         count += 1
...         total_len += len(seq)
...
>>> print("%i records with total sequence length %i" % (count, total_len))
3 records with total sequence length 75

在 Cookbook（Cookbook – 使用它的酷炫方法章節）中有更多這方面的範例，包括如何使用此程式碼片段從字串有效率地輸出 FASTQ

...
out_handle.write("@%s\n%s\n+\n%s\n" % (title, seq, qual))
...