BLASTN参数研究记之max_target_seqs和num_alignments

BLASTN参数研究记之max_target_seqs和num_alignments

1. 准备数据

见BLASTN参数研究记之max_hsps。

2. 开始研究

2.1 NCBI给出的解释

1. max_target_seqs

   Number of aligned sequences to keep. Use with report
   formats that do not have separate definition line and
   alignment sections such as tabular (all outfmt > 4). Not
   compatible with num_descriptions or num_alignments.

2. num_alignments

   Show alignments for this number of database sequences.

好像都是设置对于每个query保留多少匹配数目的，这两者有什么区别呢？

2.2 网友的答案

见于biostars

It depends on what value you set for the option outfmt (options range from 0 to 12) and you can only use either one of the max_target_seqs  or num_alignments options.
If you want only one hit:
    If you set outfmt to any option <=4
        set both -num_descriptions and -num_alignments to 1
    If you set outfmt to any option >4
         set -max_target_seqs to 1
I hope you are aware that if you want only the top hit and if your query sequence is in the database you are searching against ('nr' for example), your query sequence might be the top hit.

大概意思是-max_target_seqs和-num_alignments是功能是一样的，只是针对不同的输出格式。当-outfmt小于等于4时使用-num_alignments，当-outfmt大于4时，使用-max_target_seqs。

感觉不太对，如果是这样的话NCBI根本就没有原因设置两个参数来完成这个目的。

2.3 测试

在当前目录下新建两个文件夹 - max_target_seqs、num_alignments，为了比较这两个参数的不同，我们分别把-max_target_seqs和-num_alignments从1到500，分别运行500次，比较这500个结果之间的异同。

2.3.1 -max_target_seqs

$ for i in `seq 1 500`; do blastn -query otu.fa -outfmt '6 qseqid sseqid qstart qend sstart send evalue' -db indexes/ref -out max_target_seqs/max-target-seqs-$i.txt -max_target_seqs $i;
done

2.3.2 -num_alignments

$ for i in `seq 1 500`; do blastn -query otu.fa -outfmt '6 qseqid sseqid qstart qend sstart send evalue' -db indexes/ref -out num_alignments/num-alignments-$i.txt -num_alignments $i; done

2.4 分析

2.4.1 这两个参数是不是设置对于每个query保留多少条比对结果?

再仔细观察一下normal.txt文件，首先确定下，有多少个query有比对结果，运行：

$ awk '{print $1}' ret/normal.txt | sort -u

结果为：

ST-E00126:289:H3727ALXX:5:1106:8521:58761
ST-E00126:289:H3727ALXX:5:1111:1509:3717
ST-E00126:289:H3727ALXX:5:1112:28747:61064
ST-E00126:289:H3727ALXX:5:1117:19735:63173
ST-E00126:289:H3727ALXX:5:1119:9140:36381
ST-E00126:289:H3727ALXX:5:1201:11556:44187
ST-E00126:289:H3727ALXX:5:1213:10835:15654
ST-E00126:289:H3727ALXX:5:1213:8450:4280
ST-E00126:289:H3727ALXX:5:1215:23602:59780
ST-E00126:289:H3727ALXX:5:1219:8339:56475
ST-E00126:289:H3727ALXX:5:1220:26758:44011
ST-E00126:289:H3727ALXX:5:2102:21531:22932
ST-E00126:289:H3727ALXX:5:2110:14600:54700
ST-E00126:289:H3727ALXX:5:2120:29823:26782
ST-E00126:289:H3727ALXX:5:2121:4249:41304
ST-E00126:289:H3727ALXX:5:2123:10693:52432
ST-E00126:289:H3727ALXX:5:2219:23876:59938
ST-E00126:289:H3727ALXX:5:2224:11201:38438

发现有18条query有比对结果，normal.txt有254条结果，这其中都是一个query有多条比对结果。

那么如果-max_target_seqs和-num_alignments都是设置每个query可以保留多少条匹配结果的话，如果-max_target_seqs或者-num_alignments设置为1的话，结果肯定都是18条。

$ wc -l max_target_seqs/max-target-seqs-1.txt num_alignments/num-alignments-1.txt

结果为：

18 max_target_seqs/max-target-seqs-1.txt
18 num_alignments/num-alignments-1.txt
36 total

使用这个思路分析其他文件，发现也是成立的，那么这两个参数应该就是：

对于每个query设置可以在结果中保留多少条匹配。

2.4.2 两个参数的结果是否一样？

计算每个结果文件的md5值发现，当-max_target_seqs和-num_alignments设置为同一个值时，两个结果文件的md5值时完全一样的。

是不是可以得出结论，两个参数完全一样了呢？我觉得还不一定。可能我的测试文件特殊或是什么原因导致结果一样，如果你知道结果可以给我留言。

2.5 结论

-max_target_seqs和-num_alginments都是设置每个query保留多少条匹配结果，目前来看二者的结果是一样的。

2.6 疑问

比如我把-max_target_seqs设置为1，保留下来的是最佳匹配吗？

也就是说，BLAST去除的时候是按照匹配程度去除的吗？答案是肯定的，有兴趣的可以验证下。