BLAST短序列比对没有结果

Bug

BLAST短序列比对hg19，没有比对结果。需要设置特定参数。

运行环境、问题说明

  • 版本：ncbi-blast-2.14.0+
  • 模块：blastn
  • 参考：blast db是用hg19建的
  • 问题：比对引物序列（21~23 nt）没有结果，没有任何报错信息
  • 测试：比对稍微长点的序列（130 nt）有结果
  • 测试：同样的引物短序列在UCSC的BLAT有比对结果

解决方法

简单说明

blastn比对短序列，要设置这3个参数：

  • -dust必须设置成no，否则可能没有结果
  • -word_size设置成query序列长度的1/2
  • -evalue设置大一点，例如1000，否则可能会漏掉部分结果

详细说明

参考1

BLAST(n): No hits found 

1
2
3
4
5
6
1. Masking
Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it.
2. Score/e-value thresholds
Another source of complication is that even if a match is found, that match will have very bad scores. Both the actual score of the alignment and the e-value will be very bad. Since this is such a simple sequence, it will always score badly.
3. Word-size
The way blast works is (simplifying a bit) by finding a match for N residues (the word size) and then extending that match if extending increases the score. If your query sequence is shorter than the word size, no match will ever be found.

以上3个可能原因：

  1. 不是，query序列不是低复杂度序列
  2. 不是，测试了-evalue 1、-evalue 10，没结果
  3. 不是，测试了-word_size 4，没结果

参考2（建议仔细看这个）

Searching Short Sequences
参考2的操作说明针对的是在线版本的blast；软件对应的参数就是这3个，具体说明看blastn -help。

最终设置如下3个参数，可以比对出大部分query：

  • -dust no

    • 是否用DUST算法过滤query序列，调用的应该是dustmasker。默认是开启DUST的，默认参数是level 20 window 64 linker 1。
    • 可能是在这个window参数下，我的query太短了，所以都被滤掉了，其实根本没有进行比对，所以只设置-word_size是无效的？

  • -word_size 11

    • word_size我的理解可能是相当于比对起始的种子序列长度？
    • blastn的help里只写了需要≥4，没写默认是多少。
    • 参考2里说blastn默认的word_size是11，但是测试了不设置word_size时没有结果，设置为11时有结果。
    • 参考2里说可以调小word_size，但建议最小是query长度的1/2。
    • 测试了4和11都有结果，4会慢很多。6线程，4跑了9m57.677s，11跑了0m1.174s。

  • -evalue 1000

    • 看到有些结果的Evalue是0.003，如果设置-evalue 1e-5会漏掉这部分结果

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
sample="Test"
query_fa="Input.fa"
threads="6"

database="/xxxx/balst_db"

identity="98"
cov_hsp="95"

# Long Query Sequence
# blastn -query $query_fa -db $database -out ${sample}.noheader -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -num_threads $threads -perc_identity $identity -qcov_hsp_perc $cov_hsp -evalue 1e-5

# Short Query Sequence
blastn -query $query_fa -db $database -out ${sample}.noheader -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -num_threads $threads -perc_identity $identity -qcov_hsp_perc $cov_hsp -evalue 1000 -dust no -word_size 11

echo -e "Query Seq-id\tSubject Seq-id\tIdentity\tAlignment length\tMismatch\tGap Open\tQuery Start\tQuery End\tSubject Start\tSubject End\tEvalue\tBitscore" > ${sample}.tsv

cat ${sample}.noheader >> ${sample}.tsv

rm -f ${sample}.noheader