對測序返還結果進行處理，並進行批量blast

引言

背景介紹

由於實驗室用的seqHunter年代太老了，並且最近須要進行不少擴增，測序，因此考慮搭建一個簡單的批量化處理流程。

軟件

• python
• Windows下的cmd命令
• ncbi的blast+

---

步驟

blast+的本地化搭建

1. 在ncbi網站上下載blast+（好像blast+的運行須要JAVA環境）
2. 安裝
   • 新建db文件夾做爲數據庫文件夾
   • input和output做爲輸入輸出文件夾
     
3. 添加環境變量
   • 變量名：BLASTDB ，變量值：\path\ncbi\blast-2.7.1\db
   • 變量名：Path， 變量值：\path\ncbi\blast-2.7.1\bin
4. 構建數據庫
   • 數據準備
     • ncbi上能夠下載別人上傳的
     • 也能夠選擇把本身的幾條須要變成一個fasta文件，從而做爲建庫文件
   • 格式化數據庫(這裏我僞裝有把我全部的481同源序列變成了一個fasta文件)

   cd path\blast-2.7.1\db # 把當前目錄切到db下，這些命令都是在cmd下敲的。    
   makeblastdb -in 481.fasta -dbtype nucl -parse_seqids -hash_index -out 481_nucl

   參數說明

   -in 須要建庫的fasta序列  
   -dbtype 若是是蛋白庫則用prot，核酸庫用nucl  
   -out 所建數據庫的名稱  
   -parse_seqids 和 -hash_index 通常都加上，解釋以下:  
   -parse_seqids => Parse Seq-ids in FASTA input         # 我也不懂這是幹啥  
   -hash_index => Create index of sequence hash values

對公司返還的測序數據進行整理

1. 去除ab1結尾的文件

from path import Path
import os
a = input('input your path')
d = Path(a)
for f in d.walk():   # 遍歷你選擇的目錄下全部的文件名
    if f.endswith('ab1'): # 若是你的文件名後綴是ab1
        os.remove(f) # 移除文件

1. 將後綴名爲seq的改成fasta（由於好像要比對的文件必須是fasta格式）

for f in d.walk():
    if f.endswith('seq'):
        portion = os.path.splitext(f) # 這一步會把文件名拆成後綴加文件名的列表
        newname = portion[0] + '.fasta' # 選擇其中的文件名，而後加上fasta
        os.rename(f, newname)

利用cmd的for in循環批量進行blast

for %i in (*fasta) do blastn -query %i -db ..\\db\481_nucl -evalue 1e-3 -outfmt 7 -out ..\output\%i.txt  
# 通過實踐，最好是cd 到input文件夾下，即%i產生的文件名不要帶有路徑，否則最後的output沒法輸出

參數說明

-query 輸入文件名，也就是須要比對的序列文件  
-db 格式化後的數據庫名稱  
-evalue 設定輸出結果中的e-value閾值  
-out 輸出文件名  
-num_alignments 輸出比對上的序列的最大值條目數  
-num_threads 線程數  
此外還有：  
-num_descriptions 對比對上序列的描述信息，通常跟tabular格式連用  
-outfmt      
 0 = pairwise,  
 1 = query-anchored showing identities,  
 2 = query-anchored no identities,  
 3 = flat query-anchored, show identities,  
 4 = flat query-anchored, no identities,  
 5 = XML Blast output,  
 6 = tabular,  
 7 = tabular with comment lines,  
 8 = Text ASN.1,  
 9 = Binary ASN.1  
 10 = Comma-separated values

利用python自動對沒有比對上的輸出文件進行剔除

利用文件中包含有'0 hits found'

c = Path('F:/科研經常使用軟件/ncbi/blast-2.7.1+/output')
for f in c.walk():
    with open(f) as fh:
        a = fh.read()
        if '0 hits found' in a:
            fh.close()
            # with open 應該是等整個執行完後自動close
            # 但咱們這裏因爲程序在未結束的狀況下，就要remove文件
            # 因此會出現程序被佔用，沒法結束的狀況
            # f.close強制關閉文件，而後執行remove
            os.remove(f)

---

小結

• 聽說有個叫ncbiwww的包能夠進行批量化，下次能夠試試
• blast產生的結果仍是要好好研究下