Biopython（https://biopython.org/）是一個計算分子生物學模塊的集合，用它可以實現許多生物信息學項目中的基本任務，比如：
解析文件格式，這些信息包括如基因和蛋白質序列、蛋白質結構、PubMed記錄等等；
從資源庫下載文件，這些資源包括NCBI，ExPASy等；
運行（本地或遠程）常用的生物信息學算法，如BLAST，ClustalW等；
運行Biopython實現的算法，進行聚類、機器學習、數據分析和可視化。
我們也不僅可以用Biopython建立一個研究流程，也可以為特定的任務寫新的代碼，而讓Biopython進行更多標準操作，你還可以修改Biopython的開原始碼以更好地適應你的需求。
讓我們開始吧安裝Biopython直接使用pip安裝：
pip install biopython
升級舊版本：
pip install biopython --upgrade
在前面的筆記中我們已經學習了如何用基本的字符串操作來解析序列文件，現在讓我們用Biopython來試試看。Biopython提供了快速處理序列文件的工具，使利用不同文件格式、注釋序列記錄和把他們寫入文件等操作變得非常簡單。下面實例我們會用到4個從Bio導入的模塊：Seq用來創建序列對象；IUPAC用來定義一個序列對象用的字符集（如DNA或蛋白質）；SeqRecord允許創建一個包含ID、注釋、描述等的序列記錄數據；SeqIO用來讀寫格式化的序列文件。
from Bio import Seq
from Bio.Alphabet import IUPAC
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO
dna = open("hemoglobin-gene.txt").read().strip()
dna = Seq.Seq(dna, IUPAC.unambiguous_dna)
mrna = dna.transcribe()
protein = mrna.translate()
protein_record = SeqRecord(protein, id='sp|P69905.2|HBA_HUMAN', description="Hemoglobin subunit alpha, human")
outfile = open("HBA_HUMAN.fasta", "w")
SeqIO.write(protein_record, outfile, "fasta")
outfile.close()
hemoglobin-gene.txt文件內容如下：
ATGGTGCTGTCTCCTGCCGACAAGACCAACGTCAAGGCCGCCTGGGGTAAGGTCGGCGCGCACGCTGGCGAGTATGGTGCGGAGGCCCTGGAGAGGATGTTCCTGTCCTTCCCCACCACCAAGACCTACTTCCCGCACTTCGACCTGAGCCACGGCTCTGCCCAGGTTAAGGGCCACGGCAAGAAGGTGGCCGACGCGCTGACCAACGCCGTGGCGCACGTGGACGACATGCCCAACGCGCTGTCCGCCCTGAGCGACCTGCACGCGCACAAGCTTCGGGTGGACCCGGTCAACTTCAAGCTCCTAAGCCACTGCCTGCTGGTGACCCTGGCCGCCCACCTCCCCGCCGAGTTCACCCCTGCGGTGCACGCCTCCCTGGACAAGTTCCTGGCTTCTGTGAGCACCGTGCTGACCTCCAAATACCGTTAA
HBA_HUMAN.fasta輸出文件如下：
>sp|P69905.2|HBA_HUMAN Hemoglobin subunit alpha, human
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR*
程序首先從文本文件中讀取序列，然後將它轉錄成mRNA序列，翻譯成肽段序列，最後寫入HBA_HUMAN.fasta輸出文件。讓我們逐步分析導入Python會話3中的對象們。
Seq對象Seq.Seq類創建了一個序列對象，用戶在創建Seq對象時可以指定（或不指定）其字符集。
>>> from Bio.Seq import Seq
>>> my_seq = Seq("AGTACACTGGT")
>>> my_seq
Seq('AGTACACTGGT', Alphabet())
>>> my_seq.alphabet
Alphabet()
當我們指定字符集時：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("AGTACACTGGT", IUPAC.unambiguous_dna)
>>> my_seq
Seq('AGTACACTGGT', IUPACUnambiguousDNA())
>>> my_seq.alphabet
IUPACUnambiguousDNA()
當然，這也許是個胺基酸序列：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("AGTACACTGGT", IUPAC.unambiguous_dna)
>>> my_seq
Seq('AGTACACTGGT', IUPACUnambiguousDNA())
>>> my_seq.alphabet
IUPACUnambiguousDNA()
Biopython包含了一個預編譯好的字符集，覆蓋了所有生物序列類型。最頻繁使用的是IUPAC定義的字符集(http://www. chem. qmw. ac. uk/iupac)。如果用戶需要使用字符集，就必須從Bio. Alphabet模塊導入IUPAC模塊。它包括字符集IUPACUnamiguousDNA(基本的ACTG字母)，IUPACAmbiguousDNA(包含二義字母)，ExtendedIUPACDNA(包含修飾的鹼基)，IUPACUnamiguousRNA, IUPACAmbiguousRNA, ExtendedIUPACRNA，IUPACProtein(IUPAC標準胺基酸)和ExtendedIUPACProtein(包括硒代半胱氨酸，X等)。在上面的實例中定義的dna變量是一個以IUPAC.unambiguous_dna字符集為特徵的序列對象。
轉錄和翻譯序列Seq對象可以用transcribe()方法和translate()方法來進行轉錄和翻譯。transcribe()方法只是把所有的T替換成U，同時把字符集設置成RNA。我們也可以用reverse_complement()方法來得到反向互補序列，如：
>>> from Bio import Seq
>>> my_seq = Seq.Seq("AGTACACTGGT")
>>> cdna = my_seq.reverse_complement()
>>> cdna
Seq('ACCAGTGTACT', Alphabet())
在Biopython中，我們可以像處理字符串一樣處理序列對象。例如，可以索引、切片、分割、轉換序列大小寫，計算出現字符個數等等：
>>> from Bio import Seq
>>> my_seq = Seq.Seq("AGTACACTGGT")
>>> my_seq[0]
'A'
>>> my_seq[0:3]
Seq('AGT', Alphabet())
>>> my_seq.split('T')
[Seq('AG', Alphabet()), Seq('ACAC', Alphabet()), Seq('GG', Alphabet()), Seq('', Alphabet())]
>>> my_seq.count('A')
3
>>> my_seq.count('A')/len(my_seq)
0.2727272727272727
需要注意的是將Seq對象分割後返回的是多個Seq對象，這些Seq對象可以用+連接。
除此之外，我們還可以用find()方法搜索序列中的子串，如果沒有找到會返回-1，找到了則返回目標序列的最左端匹配字符的位置。當然，也可以配合Python的re模塊或Biopython的Bio.motif模塊用正則表達式進行搜索。
SeqRecord對象SeqRecord類提供序列及其注釋的容器，在上面的實例中我們將翻譯得到的protein變量轉換成SeqRecord對象。
protein_record = SeqRecord(protein, id='sp|P69905.2|HBA_HUMAN', description="Hemoglobin subunit alpha, human")
SeqRecord 類包括下列屬性:
.seq
– 序列自身（即 Seq 對象）。
.id
– 序列ID。通常類同於accession number。
.name
– 序列名/id 。可以是accession number, 也可是clone名（類似GenBank record中的LOCUS id）。
.description
– 序列描述。
.letter_annotations
– 對照序列的每個字母逐字注釋（per-letter-annotations），以信息名為鍵（keys），信息內容為值（value）所構成的字典。值與序列等長，用Python列表、元組或字符串表示。.letter_annotations可用於質量分數或二級結構信息 (如 Stockholm/PFAM 比對文件)等數據的存儲。
.annotations
– 用於儲存附加信息的字典。信息名為鍵（keys），信息內容為值（value）。用於保存序列的零散信息（如unstructured information）。
.features
– SeqFeature 對象列表，儲存序列的結構化信息（structured information），如：基因位置, 蛋白結構域。
.dbxrefs
– 儲存資料庫交叉引用信息（cross-references）的字符串列表。
SeqIO模塊Biopython的SeqIO模塊提供了多種常用文件格式的解析器。這些解析器從一個輸入文件（從本地或資料庫中）提取信息，而後自動轉換成SeqRecord對象。SeqIO模塊也提供了一種方法把SeqRecord對象寫入到格式化的文件中。
解析文件序列文件的解析有兩種方式：SeqIO.parse()和SeqIO.read()：這兩種方法都有兩個必須的參數和一個可選參數：
輸入文件，它指定從哪裡讀取數據；
數據格式，如fasta或genbank，完整的支持格式列表參見：http:// biopython.org/wiki/SeqIO
參數指定序列數據的字符集（可選）。
SeqIO.parse()和SeqIO.read()的區別在於：SeqIO.parse()返回的是一個迭代器，即從輸入文件中產生幾個SeqRecord對象，你可以用for或while循環來遍歷他們；而當文件中只包含一條記錄時，就像上面的例子，則必須用SeqIO.read()。換言之，SeqIO.parse()能處理輸入文件中任意數目的記錄，而SeqIO.read()只能處理一條記錄的文件，後者會先檢查文件中是否只有一條記錄，否則將會產生錯誤。
解析大文件對於大量記錄，可以使用SeqIO.index()方法，它需要兩個參數：記錄文件和文件格式。SeqIO.index()方法返回一個字典式對象，用它可以訪問所有記錄而不用把他們都讀取到內存中。字典的鍵是記錄的ID，值包含整個記錄，後者也能用屬性來訪問，如id，description等。需要注意的是，這些類字典對象是只讀的，也就是說創建後不能刪除或插入。
寫文件SeqIO.write()方法可將一個或多個SeqRecord對象寫入指定格式的文件中。用這個方法需要三個參數：一個或多個SeqRecord對象，輸出文件以及輸出的格式。
示例用SeqIO模塊來解析一個多序列FASTA文件代碼會輸出這些序列的標識、序列和長度。
from Bio import SeqIO
fasta_file = open("Uniprot.fasta","r")
for seq_record in SeqIO.parse(fasta_file, "fasta"):
print(seq_record.id)
print(repr(seq_record.seq))
print(len(seq_record))
fasta_file.close()
下面的代碼將文件解析並儲存到列表中，輸出為第一條記錄的標識和序列。
from Bio import SeqIO
uniprot_iterator = SeqIO.parse("Uniprot.fasta", "fasta")
records = list(uniprot_iterator)
print(records[0].id)
print(records[0].seq)
此外可以使用字典，鍵是記錄的ID，值包含記錄的信息：
from Bio import SeqIO
uniprot_iterator = SeqIO.parse("Uniprot.fasta", "fasta")
records = SeqIO.index("Uniprot.fasta","fasta")
print(len(records['sp|P03372|ESR1_HUMAN'].seq))
可以使用SeqIO.parse()和SeqIO.write()來轉換序列格式，下面的腳本把一個Genbank文件轉換成FASTA文件：
from Bio import SeqIO
genbank_file = open ("AY810830.gb", "r")
output_file = open("AY810830.fasta", "w")
records = SeqIO.parse(genbank_file, "genbank")
SeqIO.write(records, output_file, "fasta")
output_file.close()
猜你喜歡
生信基礎知識100講
生信菜鳥團-專題學習目錄（5）
還有更多文章，請移步公眾號閱讀
▼ 如果你生信基本技能已經入門，需要提高自己，請關注下面的生信技能樹，看我們是如何完善生信技能，成為一個生信全棧工程師。
▼ 如果你是初學者，請關注下面的生信菜鳥團，了解生信基礎名詞，概念，紮實的打好基礎，爭取早日入門。