- Biopython 教程
- Biopython - 首頁
- Biopython - 簡介
- Biopython - 安裝
- 建立簡單應用程式
- Biopython - 序列
- 高階序列操作
- 序列 I/O 操作
- Biopython - 序列比對
- Biopython - BLAST概述
- Biopython - Entrez 資料庫
- Biopython - PDB 模組
- Biopython - 基序物件
- Biopython - BioSQL 模組
- Biopython - 種群遺傳學
- Biopython - 基因組分析
- Biopython - 表型微陣列
- Biopython - 繪圖
- Biopython - 聚類分析
- Biopython - 機器學習
- Biopython - 測試技術
- Biopython 資源
- Biopython - 快速指南
- Biopython - 有用資源
- Biopython - 討論
Biopython - BLAST概述
BLAST 代表 **基本區域性比對搜尋工具**。它查詢生物序列之間相似區域。Biopython 提供 Bio.Blast 模組來處理 NCBI BLAST 操作。您可以在本地連線或透過網際網路連線執行 BLAST。
讓我們在下節中簡要了解這兩種連線 -
透過網際網路執行
Biopython 提供 Bio.Blast.NCBIWWW 模組來呼叫 BLAST 的線上版本。為此,我們需要匯入以下模組 -
>>> from Bio.Blast import NCBIWWW
NCBIWW 模組提供 qblast 函式來查詢 BLAST 線上版本,https://blast.ncbi.nlm.nih.gov/Blast.cgi。qblast 支援線上版本支援的所有引數。
要獲取有關此模組的任何幫助,請使用以下命令並瞭解其功能 -
>>> help(NCBIWWW.qblast) Help on function qblast in module Bio.Blast.NCBIWWW: qblast( program, database, sequence, url_base = 'https://blast.ncbi.nlm.nih.gov/Blast.cgi', auto_format = None, composition_based_statistics = None, db_genetic_code = None, endpoints = None, entrez_query = '(none)', expect = 10.0, filter = None, gapcosts = None, genetic_code = None, hitlist_size = 50, i_thresh = None, layout = None, lcase_mask = None, matrix_name = None, nucl_penalty = None, nucl_reward = None, other_advanced = None, perc_ident = None, phi_pattern = None, query_file = None, query_believe_defline = None, query_from = None, query_to = None, searchsp_eff = None, service = None, threshold = None, ungapped_alignment = None, word_size = None, alignments = 500, alignment_view = None, descriptions = 500, entrez_links_new_window = None, expect_low = None, expect_high = None, format_entrez_query = None, format_object = None, format_type = 'XML', ncbi_gi = None, results_file = None, show_overview = None, megablast = None, template_type = None, template_length = None ) BLAST search using NCBI's QBLAST server or a cloud service provider. Supports all parameters of the qblast API for Put and Get. Please note that BLAST on the cloud supports the NCBI-BLAST Common URL API (http://ncbi.github.io/blast-cloud/dev/api.html). To use this feature, please set url_base to 'http://host.my.cloud.service.provider.com/cgi-bin/blast.cgi' and format_object = 'Alignment'. For more details, please see 8. Biopython – Overview of BLAST https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE = BlastDocs&DOC_TYPE = CloudBlast Some useful parameters: - program blastn, blastp, blastx, tblastn, or tblastx (lower case) - database Which database to search against (e.g. "nr"). - sequence The sequence to search. - ncbi_gi TRUE/FALSE whether to give 'gi' identifier. - descriptions Number of descriptions to show. Def 500. - alignments Number of alignments to show. Def 500. - expect An expect value cutoff. Def 10.0. - matrix_name Specify an alt. matrix (PAM30, PAM70, BLOSUM80, BLOSUM45). - filter "none" turns off filtering. Default no filtering - format_type "HTML", "Text", "ASN.1", or "XML". Def. "XML". - entrez_query Entrez query to limit Blast search - hitlist_size Number of hits to return. Default 50 - megablast TRUE/FALSE whether to use MEga BLAST algorithm (blastn only) - service plain, psi, phi, rpsblast, megablast (lower case) This function does no checking of the validity of the parameters and passes the values to the server as is. More help is available at: https://ncbi.github.io/blast-cloud/dev/api.html
通常,qblast 函式的引數基本上類似於您可以在 BLAST 網頁上設定的不同引數。這使得 qblast 函式易於理解,並減少了使用它的學習曲線。
連線和搜尋
為了理解連線和搜尋 BLAST 線上版本的流程,讓我們對線上 BLAST 伺服器透過 Biopython 進行一次簡單的序列搜尋(在我們本地的序列檔案中可用)。
**步驟 1** - 在 Biopython 目錄中建立一個名為 **blast_example.fasta** 的檔案,並將以下序列資訊作為輸入
Example of a single sequence in FASTA/Pearson format: >sequence A ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc >sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc
**步驟 2** - 匯入 NCBIWWW 模組。
>>> from Bio.Blast import NCBIWWW
**步驟 3** - 使用 python IO 模組開啟序列檔案 **blast_example.fasta**。
>>> sequence_data = open("blast_example.fasta").read()
>>> sequence_data
'Example of a single sequence in FASTA/Pearson format:\n\n\n> sequence
A\nggtaagtcctctagtacaaacacccccaatattgtgatataattaaaatt
atattcatat\ntctgttgccagaaaaaacacttttaggctatattagagccatcttctttg aagcgttgtc\n\n'
**步驟 4** - 現在,呼叫 qblast 函式,將序列資料作為主要引數傳遞。另一個引數表示資料庫 (nt) 和內部程式 (blastn)。
>>> result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
>>> result_handle
<_io.StringIO object at 0x000001EC9FAA4558>
**blast_results** 儲存我們搜尋的結果。它可以儲存到檔案中以供以後使用,也可以解析以獲取詳細資訊。我們將在下一節中學習如何做到這一點。
**步驟 5** - 同樣的功能也可以使用 Seq 物件來完成,而不是使用整個 fasta 檔案,如下所示 -
>>> from Bio import SeqIO
>>> seq_record = next(SeqIO.parse(open('blast_example.fasta'),'fasta'))
>>> seq_record.id
'sequence'
>>> seq_record.seq
Seq('ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatat...gtc',
SingleLetterAlphabet())
現在,呼叫 qblast 函式,將 Seq 物件 record.seq 作為主要引數傳遞。
>>> result_handle = NCBIWWW.qblast("blastn", "nt", seq_record.seq)
>>> print(result_handle)
<_io.StringIO object at 0x000001EC9FAA4558>
BLAST 會自動為您的序列分配一個識別符號。
**步驟 6** - result_handle 物件將包含整個結果,並可以儲存到檔案中以供以後使用。
>>> with open('results.xml', 'w') as save_file:
>>> blast_results = result_handle.read()
>>> save_file.write(blast_results)
我們將在後面的章節中看到如何解析結果檔案。
執行獨立 BLAST
本節說明如何在本地系統上執行 BLAST。如果您在本地系統上執行 BLAST,它可能會更快,並且允許您建立自己的資料庫來搜尋序列。
連線 BLAST
通常,不建議在本地執行 BLAST,因為它體積龐大,需要額外的執行軟體的工作量,並且涉及成本。線上 BLAST 對於基本和高階用途來說已經足夠了。當然,有時您可能需要在本地安裝它。
假設您正在線上進行頻繁搜尋,這可能需要大量時間和高網路流量,並且如果您有專有序列資料或與 IP 相關的 issue,則建議在本地安裝它。
為此,我們需要按照以下步驟操作 -
**步驟 1** - 使用給定的連結下載並安裝最新的 blast 二進位制檔案 - ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
**步驟 2** - 使用以下連結下載並解壓縮最新的必要資料庫 - ftp://ftp.ncbi.nlm.nih.gov/blast/db/
BLAST 軟體在其網站上提供了許多資料庫。讓我們從 blast 資料庫站點下載 alu.n.gz 檔案並將其解壓縮到 alu 資料夾中。此檔案為 FASTA 格式。要在我們的 blast 應用程式中使用此檔案,我們需要首先將檔案從 FASTA 格式轉換為 blast 資料庫格式。BLAST 提供 makeblastdb 應用程式來執行此轉換。
使用以下程式碼片段 -
cd /path/to/alu makeblastdb -in alu.n -parse_seqids -dbtype nucl -out alun
執行上述程式碼將解析輸入檔案 alun 並建立 BLAST 資料庫,例如多個檔案 alun.nsq、alun.nsi 等。現在,我們可以查詢此資料庫以查詢序列。
我們已在本地伺服器上安裝了 BLAST,並且還有一個示例 BLAST 資料庫 **alun** 可以對其進行查詢。
**步驟 3** - 讓我們建立一個示例序列檔案來查詢資料庫。建立一個檔案 search.fsa 並將以下資料放入其中。
>gnl|alu|Z15030_HSAL001056 (Alu-J) AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCT TGAGCCTAGGAGTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAA AGAAAAAAAAAATAGCTCTGCTGGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTG GGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCCACGATCACACCACT GCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA >gnl|alu|D00596_HSAL003180 (Alu-Sx) AGCCAGGTGTGGTGGCTCACGCCTGTAATCCCACCGCTTTGGGAGGCTGAGTCAGATCAC CTGAGGTTAGGAATTTGGGACCAGCCTGGCCAACATGGCGACACCCCAGTCTCTACTAAT AACACAAAAAATTAGCCAGGTGTGCTGGTGCATGTCTGTAATCCCAGCTACTCAGGAGGC TGAGGCATGAGAATTGCTCACGAGGCGGAGGTTGTAGTGAGCTGAGATCGTGGCACTGTA CTCCAGCCTGGCGACAGAGGGAGAACCCATGTCAAAAACAAAAAAAGACACCACCAAAGG TCAAAGCATA >gnl|alu|X55502_HSAL000745 (Alu-J) TGCCTTCCCCATCTGTAATTCTGGCACTTGGGGAGTCCAAGGCAGGATGATCACTTATGC CCAAGGAATTTGAGTACCAAGCCTGGGCAATATAACAAGGCCCTGTTTCTACAAAAACTT TAAACAATTAGCCAGGTGTGGTGGTGCGTGCCTGTGTCCAGCTACTCAGGAAGCTGAGGC AAGAGCTTGAGGCTACAGTGAGCTGTGTTCCACCATGGTGCTCCAGCCTGGGTGACAGGG CAAGACCCTGTCAAAAGAAAGGAAGAAAGAACGGAAGGAAAGAAGGAAAGAAACAAGGAG AG
序列資料是從 alun 檔案中收集的;因此,它與我們的資料庫匹配。
**步驟 4** - BLAST 軟體提供了許多應用程式來搜尋資料庫,我們使用 blastn。**blastn 應用程式至少需要三個引數,db、query 和 out。db** 指的是要搜尋的資料庫;**query** 是要匹配的序列;**out** 是儲存結果的檔案。現在,執行以下命令執行此簡單查詢 -
blastn -db alun -query search.fsa -out results.xml -outfmt 5
執行上述命令將在 **results.xml** 檔案中搜索並給出輸出,如下所示(部分資料) -
<?xml version = "1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN"
"http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput>
<BlastOutput_program>blastn</BlastOutput_program>
<BlastOutput_version>BLASTN 2.7.1+</BlastOutput_version>
<BlastOutput_reference>Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb
Miller (2000), "A greedy algorithm for aligning DNA sequences", J
Comput Biol 2000; 7(1-2):203-14.
</BlastOutput_reference>
<BlastOutput_db>alun</BlastOutput_db>
<BlastOutput_query-ID>Query_1</BlastOutput_query-ID>
<BlastOutput_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</BlastOutput_query-def>
<BlastOutput_query-len>292</BlastOutput_query-len>
<BlastOutput_param>
<Parameters>
<Parameters_expect>10</Parameters_expect>
<Parameters_sc-match>1</Parameters_sc-match>
<Parameters_sc-mismatch>-2</Parameters_sc-mismatch>
<Parameters_gap-open>0</Parameters_gap-open>
<Parameters_gap-extend>0</Parameters_gap-extend>
<Parameters_filter>L;m;</Parameters_filter>
</Parameters>
</BlastOutput_param>
<BlastOutput_iterations>
<Iteration>
<Iteration_iter-num>1</Iteration_iter-num><Iteration_query-ID>Query_1</Iteration_query-ID>
<Iteration_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</Iteration_query-def>
<Iteration_query-len>292</Iteration_query-len>
<Iteration_hits>
<Hit>
<Hit_num>1</Hit_num>
<Hit_id>gnl|alu|Z15030_HSAL001056</Hit_id>
<Hit_def>(Alu-J)</Hit_def>
<Hit_accession>Z15030_HSAL001056</Hit_accession>
<Hit_len>292</Hit_len>
<Hit_hsps>
<Hsp>
<Hsp_num>1</Hsp_num>
<Hsp_bit-score>540.342</Hsp_bit-score>
<Hsp_score>292</Hsp_score>
<Hsp_evalue>4.55414e-156</Hsp_evalue>
<Hsp_query-from>1</Hsp_query-from>
<Hsp_query-to>292</Hsp_query-to>
<Hsp_hit-from>1</Hsp_hit-from>
<Hsp_hit-to>292</Hsp_hit-to>
<Hsp_query-frame>1</Hsp_query-frame>
<Hsp_hit-frame>1</Hsp_hit-frame>
<Hsp_identity>292</Hsp_identity>
<Hsp_positive>292</Hsp_positive>
<Hsp_gaps>0</Hsp_gaps>
<Hsp_align-len>292</Hsp_align-len>
<Hsp_qseq>
AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGAGTTTG
CGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCTGGTGGTGCATG
CCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCC
ACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA
</Hsp_qseq>
<Hsp_hseq>
AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGA
GTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCT
GGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGG
CTGTGGTGAGCCACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAAC
AAATAA
</Hsp_hseq>
<Hsp_midline>
|||||||||||||||||||||||||||||||||||||||||||||||||||||
|||||||||||||||||||||||||||||||||||||||||||||||||||||
|||||||||||||||||||||||||||||||||||||||||||||||||||||
|||||||||||||||||||||||||||||||||||||||||||||||||||||
|||||||||||||||||||||||||||||||||||||||||||||||||||||
|||||||||||||||||||||||||||
</Hsp_midline>
</Hsp>
</Hit_hsps>
</Hit>
.........................
.........................
.........................
</Iteration_hits>
<Iteration_stat>
<Statistics>
<Statistics_db-num>327</Statistics_db-num>
<Statistics_db-len>80506</Statistics_db-len>
<Statistics_hsp-lenv16</Statistics_hsp-len>
<Statistics_eff-space>21528364</Statistics_eff-space>
<Statistics_kappa>0.46</Statistics_kappa>
<Statistics_lambda>1.28</Statistics_lambda>
<Statistics_entropy>0.85</Statistics_entropy>
</Statistics>
</Iteration_stat>
</Iteration>
</BlastOutput_iterations>
</BlastOutput>
上述命令可以在 python 內部使用以下程式碼執行 -
>>> from Bio.Blast.Applications import NcbiblastnCommandline >>> blastn_cline = NcbiblastnCommandline(query = "search.fasta", db = "alun", outfmt = 5, out = "results.xml") >>> stdout, stderr = blastn_cline()
這裡,第一個是 blast 輸出的控制代碼,第二個是 blast 命令生成的可能的錯誤輸出。
由於我們已將輸出檔案作為命令列引數提供 (out = “results.xml”) 並將輸出格式設定為 XML (outfmt = 5),因此輸出檔案將儲存在當前工作目錄中。
解析 BLAST 結果
通常,BLAST 輸出使用 NCBIXML 模組解析為 XML 格式。為此,我們需要匯入以下模組 -
>>> from Bio.Blast import NCBIXML
現在,**使用 python open 方法直接開啟檔案** 並 **使用 NCBIXML parse 方法**,如下所示 -
>>> E_VALUE_THRESH = 1e-20
>>> for record in NCBIXML.parse(open("results.xml")):
>>> if record.alignments:
>>> print("\n")
>>> print("query: %s" % record.query[:100])
>>> for align in record.alignments:
>>> for hsp in align.hsps:
>>> if hsp.expect < E_VALUE_THRESH:
>>> print("match: %s " % align.title[:100])
這將產生如下輸出 -
query: gnl|alu|Z15030_HSAL001056 (Alu-J) match: gnl|alu|Z15030_HSAL001056 (Alu-J) match: gnl|alu|L12964_HSAL003860 (Alu-J) match: gnl|alu|L13042_HSAL003863 (Alu-FLA?) match: gnl|alu|M86249_HSAL001462 (Alu-FLA?) match: gnl|alu|M29484_HSAL002265 (Alu-J) query: gnl|alu|D00596_HSAL003180 (Alu-Sx) match: gnl|alu|D00596_HSAL003180 (Alu-Sx) match: gnl|alu|J03071_HSAL001860 (Alu-J) match: gnl|alu|X72409_HSAL005025 (Alu-Sx) query: gnl|alu|X55502_HSAL000745 (Alu-J) match: gnl|alu|X55502_HSAL000745 (Alu-J)