做生信的人基本都繞不開序列比對。要麼用網頁版BLAST要麼自己裝個本地BLAST。我記得第一次跑本地BLAST是為了批量比對幾百條16S序列網頁版一次只能提交幾條提交多了服務器直接報錯等半天結果還不一定能下載全。後來咬咬牙在服務器上裝了BLAST從NCBI拉了數據庫幾分鐘就把幾百條序列全部比對完了。從那以後本地BLAST就成了我日常生信分析裏最常用的工具之一。這篇文章不是照搬說明文檔而是把我從零開始摸索本地BLAST的經驗整理出來包括為什麼要用本地版、怎麼裝、怎麼建庫、怎麼跑命令、怎麼調參、怎麼避坑。不管你是剛接觸生信的新手還是只想把序列比對這件事做得更順手的老手應該都能從裏面找到有用的東西。1. 為什麼有網頁BLAST還要折騰本地BLAST1.1 網頁版很好用但是它有三個短板NCBI的網頁BLAST確實是絕大多數人接觸序列比對的第一站。界面友好輸入序列就能出結果還有可視化圖譜對單條序列、偶爾跑一跑的場景來說完全夠用。但真正把BLAST當成生產力工具之後網頁版的限制就非常明顯了。首先是批量查詢的限制。網頁BLAST本質上是一次性任務你上傳幾條序列還可以如果手裏有幾百條候選序列要逐一比對網頁版就非常痛苦。你得分批提交每批還要注意格式結果還要手動下載、手動整理整個流程重複、機械、容易出錯。其次是數據庫限制。網頁BLAST比對的是NCBI官方維護的數據庫比如nr、nt、refseq這些。但很多時候我們需要比對自己課題組的私有數據比如某個物種的全基因組、某個轉錄組組裝結果、一批自己測序得到的序列。這些數據不可能上傳到公共網站也不應該上傳。在這種場景下本地BLAST是唯一合理選擇。第三是安全和隱私。有些序列數據可能涉及未發表的研究成果、臨床樣本、商業合作項目按規定不能上傳到外部服務器。把數據留在本地在自己的機器上跑比對既符合數據管理要求又能避免不必要的風險。1.2 本地BLAST真正解決的是什麼本地BLAST說白了就是兩件事把數據庫放在自己電腦或服務器上把BLAST程序安裝在自己環境裏然後用命令行自由控制比對過程。它解決的不是能不能比對的問題而是能不能按你想要的方式、在你能掌控的範圍內、可重複地進行比對的問題。舉個例子我曾經要做一批基因家族成員的鑑定需要把某個物種的全部蛋白序列和擬南芥的已知蛋白家族序列進行雙向比對。用網頁版做雙向比對先要準備查詢序列和目標數據庫一次提交幾百條然後等待、下載、整理來回折騰好幾天。用本地BLAST寫個循環腳本半小時就能跑完而且結果格式統一方便後續用Python或R批量處理。另外本地BLAST的可重現性也是網頁版比不了的。你可以在命令行裏精確記錄用了哪個版本的BLAST、哪個版本的數據庫、什麼參數組合這樣別人復現你的分析時結果完全一樣。學術發表時這點特別重要。1.3 哪些人需要用本地BLAST說實話不是所有人都需要本地BLAST。如果你只是偶爾查一條序列確認一下它是什麼基因那網頁版足夠了沒必要折騰。但如果你符合下面任何一條我建議你早日擁抱本地BLAST每天或每週都要做序列比對且查詢序列不止一兩條需要比對自建數據庫比如物種特有基因、轉錄本序列、變異位點附近的序列課題涉及大量序列篩選例如基因家族鑑定、系統發育分析前的序列提取對結果格式有要求比如需要表格形式批量處理參與的項目對數據隱私敏感不方便上傳到公共服務器。一句話網頁版適合偶爾查一下本地BLAST適合長期、批量、定製化。2. 從零搭一個本地BLAST環境版本選擇與安裝2.1 BLAST和舊版BLAST的關係先說個容易搞混的地方。NCBI在2009年前後發布了BLAST套件替代了早期的legacy BLAST。現在提到本地BLAST默認都是指BLAST這一套命令行工具通常是blastn、blastp、blastx、tblastn、tblastx、makeblastdb等。舊版的blastall和formatdb已經被淘汰了新手完全沒必要去接觸。BLAST的核心優勢是模塊化。不同的比對任務用不同的程序參數更清晰輸出格式更多樣包括JSON、XML、CSV等而且支持多線程。我們日常用的主要是以下幾個blastn核酸序列 vs 核酸數據庫blastp蛋白序列 vs 蛋白數據庫blastx核酸序列翻譯成蛋白 vs 蛋白數據庫tblastn蛋白序列 vs 核酸數據庫翻譯成的蛋白tblastx核酸序列翻譯 vs 核酸數據庫翻譯最慢一般少用makeblastdb把fasta格式的序列文件變成BLAST可查詢的數據庫安裝前先想清楚你的操作系統和包管理方式。絕大多數生信環境是Linux服務器但macOS和Windows也可以裝後面分開說。2.2 Linux環境安裝apt/conda在Linux上安裝BLAST有兩種主流方式系統包管理器和conda。我最推薦conda因為它不依賴系統權限可以安裝指定版本環境隔離乾淨後續如果要裝其他生信工具也方便。用conda安裝只需要一行命令conda create -n blast -c bioconda blast conda activate blast這會創建一個名為blast的虛擬環境並安裝最新版的BLAST。如果你想要特定版本可以這樣指定conda create -n blast -c bioconda blast2.14.0如果你用的是Ubuntu/CentOS這類系統也可以用apt或yum直接裝但版本通常比較舊不建議。比如Ubuntu 20.04默認倉庫裏的blast版本可能是2.9.0而NCBI早就更新到2.14甚至更高了。舊版本在參數、輸出格式、性能上都有差異能用新版盡量用新版。不過有些集群環境不允許聯網或者管理員不讓你用conda那就只能手動下載。NCBI的FTP站點上有預編譯的Linux二進制包下載解壓後把bin目錄加入PATH即可。wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast/LATEST/ncbi-blast-2.14.1-x64-linux.tar.gz tar zxvf ncbi-blast-2.14.1-x64-linux.tar.gz echo export PATH$PATH:/path/to/ncbi-blast-2.14.1/bin ~/.bashrc source ~/.bashrc2.3 macOS/Windows怎麼裝macOS上可以用conda也可以用Homebrew。個人推薦conda因為和Linux風格一致命令完全通用。brew install blast或者conda create -n blast -c bioconda blastWindows用戶稍微麻煩一點因為BLAST原本是在Unix環境下開發的。現在最常見的辦法是在Windows Subsystem for Linux裏面安裝比如Ubuntu子系統然後用和Linux一樣的方式裝。如果你不想用WSL也可以直接下載NCBI提供的Windows安裝包.exe那種安裝後在命令提示符裏運行但路徑處理、多線程支持方面會有一些小坑我建議能用WSL就別用原生Windows後續寫腳本也更順手。2.4 安裝後第一件事確認PATH裝好後別急著建庫先確認一下程序是否可用。打開終端輸入blastn -version正常的話會輸出類似下面的信息blastn: 2.14.1 Package: blast 2.14.1, build Sep 9 2023如果提示command not found說明bin目錄沒有加入PATH。要麼重新檢查安裝路徑要麼手動export PATH。如果是在conda環境裏激活環境後應該自動解決這個問題。另外我習慣檢查一下makeblastdb是否也存在因為它和blastn在同一個包裏但有時手殘安裝了不完整的二進制版本就會缺。檢查命令which makeblastdb確認這兩個命令都能跑環境就準備好了。3. 建立自己的比對數據庫makeblastdb細節3.1 數據庫是什麼為什麼不能直接用fasta很多新手第一次建庫時會有個疑問我明明有一個fasta文件裏面全是序列為什麼還要建庫直接把這個fasta當數據庫不行嗎不行。BLAST在比對時需要對數據庫建立索引這樣才能快速定位到可能的匹配區域而不是把每條查詢序列和數據庫裏每一條序列逐個做全局比對。makeblastdb做的事情就是把fasta文件裏的序列拆分成可搜索的索引結構生成一系列後綴名為.nhr、.nin、.nsq、.phr、.pin、.psq等的文件。這些文件就是BLAST的數據庫。所以建庫之後原來的fasta文件可以留作備份但BLAST程序直接讀取的是那些索引文件。移動數據庫目錄時要把整個目錄一起移動或者用絕對路徑否則容易出現找不到數據庫的錯誤。3.2 構建蛋白庫和核酸庫的參數差異makeblastdb的用法很簡單但要注意區分核酸庫和蛋白庫因為選項不同。構建核酸數據庫makeblastdb -in nucleotide.fasta -dbtype nucl -out my_nucl_db -parse_seqids構建蛋白數據庫makeblastdb -in protein.fasta -dbtype prot -out my_prot_db -parse_seqids其中-in指定輸入fasta文件-dbtype必須明確寫nucl或prot這是最容易出錯的地方如果你告訴它是nucl但文件實際是蛋白序列程序會報錯-out指定輸出數據庫的前綴名這個前綴名後面會被blastn/blastp當作-db參數-parse_seqids這個參數建議加上它會解析fasta文件中後面的序列ID方便後續用blastdbcmd按ID提取序列如果你不加序列ID會被忽略後續檢索會有問題。還有一個常用選項是-title可以給數據庫加一個描述性的標題方便自己在多個數據庫之間區分。3.3 常用數據來源NCBI下載、自建物種庫建庫的第一步是先有序列數據。常見來源有這麼幾種NCBI的refseq/nr/nt等公共數據庫可以在NCBI的FTP站點上直接下載已經格式化好的BLAST數據庫或者下載fasta文件自己建庫。下載預建庫的好處是省去建庫時間但是佔空間大而且更新需要重新下載。Ensembl/UCSC等基因組數據庫通常提供物種全基因組、全轉錄組、蛋白組的fasta文件下載後自己建庫。自己的測序結果無論是拼接出來的contig還是註釋出來的蛋白序列想比對先建庫。以人類基因組為例如果你只想比對人類的cDNA或蛋白序列可以去Ensembl下載Homo_sapiens.GRCh38.cdna.all.fa和Homo_sapiens.GRCh38.pep.all.fa然後用makeblastdb建庫。下載時注意選擇合適的版本和物種。ncbi下載預建庫的方式是update_blastdb.pl --decompress nt但你首先得安裝update_blastdb.pl腳本它包含在BLAST安裝包中。不過這個腳本在國內網絡環境下可能速度不理想如果只需要特定物種個人建議直接下載fasta自己建庫更可控。3.4 實操用一個小fasta文件建庫我在這裏做個完整的演示方便你直接跟著操作。假設我有一個名為test_sequences.fasta的文件內容是幾條短序列seq1 ACGTACGTACGTACGTACGT seq2 ACGTACGTACGTACGTACGTACGTACGT seq3 TTTTACGTACGTACGTACGTACGT現在把它建成核酸庫makeblastdb -in test_sequences.fasta -dbtype nucl -out test_db -parse_seqids運行後如果看到類似下面輸出就表示建庫成功Building a new DB, current time: ... New DB name: test_db Type: Nucleotide Added 3 sequences in 0.01 seconds.然後我們可以用blastdbcmd查看數據庫裏的內容blastdbcmd -db test_db -info輸出會顯示數據庫名稱、序列條數、總長度等信息這是驗證數據庫是否正確的重要手段。4. 跑一次真實的本地BLAST常用命令與輸出解讀4.1 blastn和blastp基本語法數據庫建好之後就可以跑比對了。最基礎的blastn命令長這樣blastn -query my_query.fasta -db test_db -out result.txt其中-query是查詢序列文件fasta格式-db是數據庫名就是makeblastdb時-out指定的前綴-out是結果文件名。蛋白比對換成blastpblastp -query protein_query.fasta -db my_prot_db -out protein_result.txt注意查詢序列的類型要和程序匹配不能用蛋白序列跑blastn也不能用核酸序列跑blastp。如果你有一個基因序列想搜索蛋白庫可以用blastx它會先把核酸翻譯成六個閱讀框的蛋白序列再比對。4.2 輸出格式怎麼選pairwise/table/json等默認輸出格式是pairwise就是那種傳統的帶比對對齊的文本結果適合人眼閱讀但不適合程序處理。如果要做批量分析我推薦用表格格式或JSON格式。用-outfmt參數控制輸出格式。常用的格式代碼有0pairwise默認帶具體對齊信息6tabular製表符分隔的表格每一行是一條命中記錄7帶註釋的tabular比6多了幾行說明5XML格式13JSON格式。比如我想輸出成表格方便用Excel或Python讀取blastn -query query.fasta -db test_db -out result.tsv -outfmt 6如果還想讓表頭帶上列名可以用blastn -query query.fasta -db test_db -out result.tsv -outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore這個寫法等價於-outfmt 6但顯式列出了14列順序依次是查詢序列ID、數據庫序列ID、百分比一致性、比對長度、錯配數、缺口數、查詢起始、查詢結束、數據庫序列起始、數據庫序列結束、E值、分數。這樣列名清晰後續處理不容易搞混。4.3 解析關鍵欄位identity, e-value, score, coverage不管用哪種格式最終都要看幾個核心指標來判斷比對結果是否可信。Percentage identitypident查詢序列和數據庫序列在比對區域的鹼基/氨基酸一致程度。一般來說同源基因的pident越高親緣關係越近。但這個值不能孤立看因為它只統計比對區域如果比對長度很短哪怕identity很高也不能說明問題。E-value期望值表示在隨機情況下你觀察到這個得分或更高得分的比對次數期望。E-value越小越顯著一般低於1e-5可以算顯著。E-value受到數據庫大小影響數據庫越大同樣得分對應的E-value越大越不顯著所以跨數據庫比較E-value要小心。Bit score標準化的得分不受數據庫大小影響但受比對矩陣和gap罰分影響。通常分數越高越可靠。兩個不同比對之間的bit score可以直接比較這點比E-value好。Coverage覆蓋率即查詢序列被比對上的區域佔總長度的比例。這在檢查是否只有部分序列比對上時很有用。比如你有一條全長基因結果顯示只有100bp比對到數據庫即使那100bp identity是100%也不能宣稱這就是全長基因。這些指標在表格輸出裏都有對應的列看懂它們是解讀BLAST結果的基本功。4.4 實例鑑定一個未知序列我來演示一個最常見場景拿到一條未知的DNA序列想知道它是什麼基因。假設序列存在unknown.fasta裏本地的nr庫或自建庫叫做refseq_nucl。blastn -query unknown.fasta -db refseq_nucl -out unknown_vs_refseq.tsv -outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore stitle -evalue 1e-5 -num_threads 4我加了stitle列這樣結果裏會直接顯示數據庫序列的描述信息方便一眼看出它是什麼。運行後用sort按evalue排序就能看到最顯著的命中sort -k12,12g unknown_vs_refseq.tsv | head -20注意上面的列我加了stitle之後它變成第15列所以排序的k參數要改成15。實際按你輸出的列數來調整。如果排在最前面的結果描述是某個已知基因且evalue小於1e-50identity超過90%那基本可以鎖定序列的歸屬。5. 查詢序列多、數據庫大時怎麼調優5.1 多執行緒參數num_threads的正確用法BLAST支持多線程加-num_threads參數即可。比如blastn -query many_queries.fasta -db big_db -num_threads 8 -out result.tsv -outfmt 6這裡要注意-num_threads的作用是加速單次查詢任務而不是簡單地把多個查詢分配到不同線程。BLAST在處理大批量查詢時如果-num_threads設置得過大反而可能因為線程間通信開銷導致性能下降。我個人的經驗是對於比較大的數據庫線程數設置為CPU物理核數的50%到75%比較合理。例如8核機器用4到6個線程16核機器用8到12個。另外如果你是並行跑多個BLAST任務比如用xargs或snakemake那就不要每個任務都開太多線程避免CPU超訂閱。通常每個任務2到4個線程同時跑多個任務整體吞吐量更高。5.2 拆分查詢與task選擇當查詢序列數量特別多比如上萬條直接一次性跑可能耗時很長。一個實用的優化方法是把查詢文件拆分成小塊然後並行運行。split -l 1000 many_queries.fasta query_part_這樣會生成query_part_aa、query_part_ab等文件每個文件1000條序列。然後用一個簡單的for循環批量跑for f in query_part_*; do blastn -query $f -db big_db -out $f.tsv -outfmt 6 -num_threads 4 done wait最後把結果cat起來cat query_part_*.tsv all_result.tsv這樣做的好處是即使某個任務意外失敗也不會影響其他任務方便重跑。另外blastn還有一個-task參數可以選擇不同比對模式。比如用blastn默認模式或者用megablast模式適合高度相似的序列速度快或者用dc-megablast適合跨物種的保守區域搜索。如果你的目的是尋找近緣序列-task megablast往往比默認blastn快很多。blastn -query query.fasta -db db -task megablast -out result.tsv -outfmt 65.3 -evalue和-word_size的影響E值這類參數不僅影響結果過濾還影響運行速度。提高E值閾值比如從默認10降到1e-5可以過濾掉大量隨機匹配減少輸出量但並不一定能顯著加速搜索因為BLAST在種子搜索階段並不完全依賴E值過濾。影響速度更明顯的是-word_size。blastn的默認word_size是11鹼基如果你確定查詢序列和數據庫序列非常相近可以增大word_size比如設為16或20這樣可以更快速跳過不相關的區域大幅縮短搜索時間。但代價是可能漏掉一些較遠的匹配。所以這個參數要根據你的實際場景權衡。蛋白比對的word_size默認是3或6不同任務不同調大也會加速但降低敏感性。5.4 用dust/seg遮蔽重複序列重複區域比如基因組裏的轉座子、低複雜度區域在比對時會產生大量非特異性的高分數命中浪費計算資源還會誤導結果解讀。BLAST默認會用低複雜度過濾blastn對核酸序列使用DUST程序blastp對蛋白序列使用SEG程序。這些過濾默認是開啟的但在某些情況下你會想關閉它比如研究序列本身就包含重複單元。關閉方式是在命令行加-dust no或-seg noblastn -query repetitive_query.fasta -db db -dust no -out result.tsv -outfmt 6如果保留默認遮蔽輸出結果裏的那些匹配區域可能都被截斷了你看到的identity會比真實值低因為遮蔽區域不參與比對。在做變異檢測之類的敏感分析時需要特別注意過濾對結果的影響。6. 我踩過的坑與解決方案6.1 數據庫名稱找不到db path寫錯這個錯誤是我見過最多的。明明makeblastdb成功了跑blastn卻提示BLAST Database error: No alias or index file found for nucleotide database test_db in current directory.原因通常是-out指定了帶路徑的名稱但blastn運行時卻用了相對路徑。比如你在/data/blast目錄下執行makeblastdb -in sequences.fasta -dbtype nucl -out /data/blast/my_db然後你回到家目錄跑blastn -query q.fasta -db my_db -out out.tsv系統就在當前目錄找my_db自然找不到。正確做法是在-db參數裏給出完整路徑blastn -query q.fasta -db /data/blast/my_db -out out.tsv另外如果數據庫文件存在但缺少某些後綴文件比如.nin或.nsq也會報錯。這種情況通常是因為建庫過程被打斷或者手動刪了某個文件。重新跑一次makeblastdb即可。6.2 記憶體不足報錯BLAST本身不算特別吃內存但當數據庫特別大比如整個nt庫或者查詢序列很長時仍然可能出現內存不足。常見報錯包括Out of memory或者std::bad_alloc解決思路有幾個。第一如果同時跑了很多個BLAST任務先減小並行任務數避免內存疊加。第二用-max_target_seqs參數限制每個查詢輸出的比對數量比如設為1或5可以大幅減少內存中保存的候選匹配。注意-max_target_seqs在BLAST 2.14之後有一些行為變化它不再嚴格限制最終報告的比對數但對性能優化還是有效的。blastn -query q.fasta -db big_db -max_target_seqs 5 -out out.tsv -outfmt 6第三如果查詢序列本身特別長比如幾百kb可以考慮分段比對或者切分查詢序列。BLAST對長序列的支持雖然沒有問題但計算量會隨長度增加。6.3 查不到任何結果時先檢查這三個參數跑出來的結果文件是空的或者只有一行表頭這種情況很多人遇到過。別慌按順序檢查三個地方。第一檢查數據庫類型是否正確。如果你用blastp查詢蛋白序列但數據庫是核酸庫那肯定沒結果。反之一樣。用blastdbcmd -db db_name -info可以查看數據庫類型。第二檢查evalue閾值是否太嚴格。默認evalue是10如果你手動設成1e-100那幾乎所有結果都會被過濾掉。調大一點再試。第三檢查word_size是否太大。如果你設了word_size 20但查詢序列和數據庫之間相似性不高種子就找不到自然就沒有後續的延伸比對。調回默認值試試。還有一個容易忽略的點查詢序列裏如果含有特殊字符比如*或-可能導致序列無法解析。檢查fasta文件格式確保只有標準的ACGT或氨基酸字母。6.4 奇奇怪怪的格式輸出問題有人用-outfmt 6時結果文件裏某些論文描述是從NCBI下載的序列自帶的比如lcl|前綴這沒問題。但有時你想把序列description也輸出到表格裏就需要手動加stitle列。問題來了如果description裏有製表符或換行符會破壞表格結構。我遇到過某條序列的description裏居然有回車符導致輸出的表格歪掉。解決辦法是在建庫時不要保留過長的description或者下載數據後先清理一下。可以用awk或sed把description裏的特殊字符去掉。另外有些人運行blastn時想要XML格式-outfmt 5但解析XML發現裏面的Query長度、Hsp長度這些標籤不熟悉。我建議想進一步處理結果的話直接用JSON格式-outfmt 13配合Python的json模塊非常方便。示例import json with open(result.json) as f: data json.load(f) for result in data[BlastOutput2][report][results][search]: for hit in result[hits]: for hsp in hit[hsps]: print(result[query_id], hit[description][0][title], hsp[identity], hsp[evalue])這樣提取信息非常靈活不用再去跟製表符搏鬥。最後再分享一個小技巧本地BLAST用得越久越能體會到工欲善其事必先利其器這句話。如果只是偶爾跑一次前面提到的基本命令完全夠用但如果要頻繁使用我強烈建議你把自己的常用數據庫整理成一個目錄結構比如~/blastdb/下按物種或按類型分好子目錄再寫幾個簡單的shell腳本或者用Snakemake/Nextflow把建庫、比對、提取結果的流程固化下來。這樣不管是三個月前還是三年後的自己只要看到你的腳本都能立即重現整個分析過程。我自己的習慣是每建完一個數據庫就順手寫一行README說明數據來源、版本、下載日期、建庫命令放在同一目錄下。這個習慣在論文回顧或者課題組交接的時候幫我節省了不少時間。BLAST本身是一個非常成熟的老工具但用好它仍然需要你在自己的工作場景裏反覆錘煉。希望這篇文章能幫你少走彎路。