NCBI官網：https://www.ncbi.nlm.nih.gov/
BioPython官方文檔：https://biopython-cn.readthedocs.io/zh_CN/latest/
bioPython是生物資訊領域用于基因序列比對的一個工具，在搭建本地化blast的時候很實用。
以下筆記内容是在開發過程中參考官方文檔 第七章：BLAST 後的相關問題及解決記錄。
首先，使用NCBIApplications編寫blast指令行
blastx_cline = NcbiblastnCommandline(cmd='blastn',query='file.txt',db='nt',evalue=0.01,outfmt=5,out='output.xml')
blastx_cline.set_parameter("task",'megablast')
blastx_cline.set_parameter("num_threads",100)
blastx_cline.set_parameter("num_alignments",100)
其中有個小問題，
-task
參數我一直加不進去，這個參數用于不同的查詢方式，即
-task megablast/dc-megablast/blastn
，運作後根據報錯資訊，我進了源碼檢視，結果發現它與衆不同……相關的set方法沒找到正确的将其置空的位置，于是我心一橫，直接改源碼（。
改完後正常了，由于目前應用有限，暫時沒有bug，這個等待後續更新了。
blast跑一次大概20s，我掐着計時器算的，同一個序列檔案不論是
megablast
還是
dc-megablast
或者是
blastn
都是差不多時間，這塊大概也需要後續更多的測試。
在跑完blast資料後，下一步就是解析了。從生成的xml檔案提取資料，其實一開始沒啥頭緒，看着那HTML标簽我甚至以為要搞爬蟲。。不過并不用，事實上，bioPython已經提供了相應的解析工具，而這裡我用到了兩個：
SearchIO
和
NCBIXML
。
results_io = SearchIO.parse('output.xml',"blast-xml")
for result in results_io:
seq_map_id = "".join([result.id ,";"])
for hsp in result.hits:
hsp_id = hsp.id
hsp_description = hsp.description
之是以使用兩個解析器，是因為
NCBIXML
解析出來的
list
并沒有按照
query id
羅列。在使用
blast
比對序列的時候使用的
fastq
檔案内其實很經常會有多個序列，如果沒有
query id
幫忙區分，會造成很大的困擾。相當于一排沒有标簽的書架，我知道我要的書就在裡面，但我沒有指引，隻能一排排書架去找，這很浪費時間。是以我先用
SearchIO
取出
query id
和
description
。其實可以的話，我應該也能夠在
SearchIO
裡找到我需要的屬性元素，然鵝大腦當機沒看源碼。。也沒找到詳細的UML圖給我取用他的元素。這個後面再研究下。
然後再用NCBIXML解析，前面說到的UML圖，其實BioPython官方文檔裡的也有點小坑，不過不妨礙事，源碼一看即得x（是以說上面的
SearchIO
也完全可以看源碼取屬性的。。大意了。。
result_handle = open('output.xml')
blast_records = NCBIXML.parse(result_handle)
for blast_record in blast_records:
for alignment in blast_record.alignments:
title = alignment.title.split(" ")[0]
length = alignment.length
accession = alignment.accession
hit_def = alignment.hit_def
hit_id = alignment.hit_id
for hsp in alignment.hsps:
bits = str('%.1f' % hsp.bits)
score = str('%.1f' % hsp.score)
identities = hsp.identities
positives = str(hsp.positives)
gaps = str(hsp.gaps)
strand = str(hsp.strand)
frame = str(hsp.frame)
query = str(hsp.query)
query_start = str(hsp.query_start)
query_end = str(hsp.query_end)
match = hsp.match
sbjct = str(hsp.sbjct)
sbjct_start = str(hsp.sbjct_start)
sbjct_end = str(hsp.sbjct_end)
num_alignment = str(hsp.num_alignments)
align_len = hsp.align_length
# 這個percentage對應NCBI上的identities百分比
percentage = '%.2f' % (float(identities)/align_len*100)