大部分時候，我們都是看著別人的教程吹埠，然后嘗試處理自己的數(shù)據(jù)第步，結(jié)果跑完了，如果和預期相符合就不會懷疑這個工具有啥問題缘琅。如果你要學習生物信息學粘都，那么有一個信條一定要記住，不要盲目相信軟件的輸出結(jié)果刷袍，多多檢查翩隧。

BLAST作為最常用的一個序列比對工具，對于大多數(shù)用戶而言呻纹，都是用于找找一個基因在其他物種的同源基因堆生，或者看看一條引物在基因組上是否有多處匹配专缠，大多不會限制輸出的匹配數(shù)目。

但是對于一些分析流程淑仆，你可能希望只要輸出一個最佳匹配就好涝婉，這個時候你可能會用到兩個參數(shù)，一個是-max_target_seqs蔗怠，另一個是-num_alignments, 它們的參數(shù)說明如下

 -max_target_seqs <Integer, >=1>                                                                                                                
   Maximum number of aligned sequences to keep                                                                                                  
   Not applicable for outfmt <= 4
  * Incompatible with:  num_descriptions, num_alignments
 -num_alignments <Integer, >=0>                                                                                                                 
   Number of database sequences to show alignments for
   * Incompatible with:  max_target_seqs

如果你認為" aligned sequences" 或者"database sequences" 指的每個序列只會輸出一個匹配的話墩弯，那么你掉進坑了。

首先給定一個檢索序列, 命名為"query.fa"

>query1
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC

選擇擬南芥的前1000行序列作為 測試集1

head -n 1000 Athaliana.fa > test1.fa
makeblastdb -in  test1.fa -out test/test1  -dbtype nucl

直接搜索寞射，找不到任何的結(jié)果

blastn -query query.fa -db test/test1 -outfmt 6

將搜索的序列加入 測試集1渔工，得到 測試集2

cat query.fa test1.fa > test2.fa
makeblastdb -in test2.fa -out test/test2  -dbtype nucl

此時搜索，可以找到一個結(jié)果

blastn -query query.fa -db test/test2 -outfmt 6

如果將同一段序列怠惶，重復多次涨缚，中間加入一些無意義的序列（如下），之后再加入到 測試集1 中策治，得 測試集3

>query2
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGCCCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT
GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGTCTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC
CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT
GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGT
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC

此時搜索脓魏，會得到三個結(jié)果

cat query2.fa test1.fa > test3.fa
makeblastdb -in test3.fa -out test/test3 -dbtype nucl
blastn -query query.fa -db test/test3 -outfmt 6
# 結(jié)果如下
query1  query2  100.000 100 0   0   1   100 1   100 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 812 911 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 1623    1722    8.11e-50    185

加入我希望只要一個結(jié)果，根據(jù)我之前的經(jīng)驗通惫，我加入一個參數(shù)-num_alignments 1, 幫助里寫的是 "Number of database sequences to show alignments for" 我一直認為是茂翔，只返回第一個找到的序列，但是運行如下代碼履腋，結(jié)果出乎我意料

blastn -query query.fa -db test/test3 -outfmt 6 -num_alignments 1
# 結(jié)果如下
query1  query2  100.000 100 0   0   1   100 1   100 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 812 911 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 1623    1722    8.11e-50    185

他居然會得到三個結(jié)果. 為什么會出現(xiàn)這種情況珊燎？

我的猜想是：這里的 Number of database sequences 的 sequences 不是匹配上的序列，而是匹配上的序列對應所在整條序列遵湖。那么根據(jù)我的猜想悔政， 當數(shù)據(jù)庫中中有A,B,C三條序列，其中A有三個潛在匹配延旧，B有兩個潛在匹配谋国，C沒有匹配，A, B, C在建立搜索數(shù)據(jù)庫的順序會影響最終的輸出結(jié)果迁沫，也就是當-num_alignments 1時芦瘾，blastn搜索到A后，就會輸出A里面的三個檢索結(jié)果就不再繼續(xù)集畅，如果先搜索到B近弟，就會輸出一個匹配，不再繼續(xù)挺智。當-num_alignments 1會輸出四個匹配祷愉。

為了驗證假設，我將query.fa 和 query2.fa 以不同排列順序構(gòu)建測試數(shù)據(jù)集，發(fā)現(xiàn)輸出結(jié)果的確和順序有關(guān)谣辞。

#順序1
cat query.fa query2.fa test1.fa > test4.fa
makeblastdb -in test4.fa -out test/test4 -dbtype nucl
blastn -query query.fa -db test/test4 -outfmt 6 -num_alignments 1
# 結(jié)果如下
query1  query2  100.000 100 0   0   1   100 1   100 8.12e-50    185
query1  query2  100.000 100 0   0   1   100 812 911 8.12e-50    185
query1  query2  100.000 100 0   0   1   100 1623    1722    8.12e-50    185
# 順序2 
cat   query2.fa  query.fa   test1.fa > test5.fa
makeblastdb -in test5.fa -out test/test5 -dbtype nucl
blastn -query query.fa -db test/test5 -outfmt 6 -num_alignments 1
# 結(jié)果如下
query1  query1  100.000 100 0   0   1   100 1   100 8.12e-50    185

這個情況是我用一條短序列搜索全基因組數(shù)據(jù)時出現(xiàn)(如果是CDS序列庫則不會出現(xiàn)一條序列出現(xiàn)多個匹配的情況), 在你設計引物的時候迫摔，如果在一條染色體上能夠找到多個相似序列，那么這種引物序列就是不合格的泥从。

但是如果你一定要保證每條序列在每個染色體都只要第一個高質(zhì)量匹配，那么如何解決呢沪摄？

我們用到一個不怎么起眼的參數(shù)躯嫉，-max_hsps 1， 根據(jù)說明"Set maximum number of HSPs per subject sequence to save for each query", 也就是保留query中最佳的前幾個匹配

繼續(xù)測試杨拐，隨機刪除query2的一些序列得到query3祈餐，保證最后一個是最佳匹配，用來驗證返回的是一個序列中最好連配哄陶，而不是位置相關(guān)結(jié)果

>query3
CTCAAGAAAAGGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGCCCCTAAACCCTAAACCCTAAACCCT    AAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT]GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGTCTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTC    TCTGTGAATTTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC
CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT
GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGT
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC

構(gòu)建索引數(shù)據(jù)庫帆阳，進行檢索，會找到每個序列中的最好匹配屋吨，如果你設置了-num_alignments 1蜒谤，那么最后結(jié)果就真的只有一個。

#順序1
cat query.fa query3.fa test1.fa > test5.fa
makeblastdb -in test5.fa -out test/test5 -dbtype nucl
blastn -query query.fa -db test/test5 -outfmt 6 -max_hsps 1
# 結(jié)果如下
query1  query3  100.000 100 0   0   1   100 1606    1705    8.11e-50    185
query1  query1  100.000 100 0   0   1   100 1   100 8.11e-50    185

最后的結(jié)論如下：