BLASTN参数研究记之max_hsps
最近频繁使用到BLASTN，BLASTN有很多参数。有一些显而易见的参数，例如：
-query-out-outfmt-db
也有一些参数不那么常用，今天主要研究研究这些不常用的参数。
1. 准备数据
主要用的数据有：
假定这些文件都放置在当前文件夹，并且当前文件夹有一个ret的子文件夹，例如，运行命令：
$ blastn -query otu.fa -outfmt '6 qseqid sseqid qstart qend sstart send evalue' -db indexes/ref -out ret/normal.txt-out-fmt '6 qseqid sseqid qstart qend sstart send evalue的意思是：
- 结果文件以
tab键分割 - 每行共有
7列 - 每一列分别是
query序列ID(qseqid)、subject序列ID(sseqid)、query序列匹配的开始位置(qstart)、query序列匹配的终止位置(qend)、subject序列匹配的开始位置(sstart)、subject序列匹配的终止位置(send)、这次匹配的e-value(evalue)
把这个命令的结果normal.txt作为标准，用来检测添加其他参数的影响。
normal.txt的内容如下：
ST-E00126:289:H3727ALXX:5:1106:8521:58761 23 122 276 661 815 2.13e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 24 122 275 661 814 7.66e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 15 122 275 661 814 7.66e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 11 122 278 651 808 7.66e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 22 122 268 661 807 2.76e-49这里只粘贴了前5行，原始文件共有254行。
那么当前文件夹的目录结构应该是：
test-blastn
├── indexes
│ ├── ref.nhr
│ ├── ref.nin
│ ├── ref.nog
│ ├── ref.nsd
│ ├── ref.nsi
│ └── ref.nsq
├── otu.fa
└── ret
└── normal.txt2. 开始研究
2.1 NCBI给出的解释
Set maximum number of HSPs per subject sequence to save for each query意思是对于每个subject，保留多少HSP。
2.2 测试
运行命令：
$ blastn -query otu.fa -outfmt '6 qseqid sseqid qstart qend sstart send evalue' -db indexes/ref -out ret/max-hsps-1.txt -max_hsps 1max-hsps-1.txt共有244行，前5行如下：
ST-E00126:289:H3727ALXX:5:1106:8521:58761 23 122 276 661 815 2.13e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 24 122 275 661 814 7.66e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 15 122 275 661 814 7.66e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 11 122 278 651 808 7.66e-50
ST-E00126:289:H3727ALXX:5:1106:8521:58761 22 122 268 661 807 2.76e-492.3 分析
直观结论是，因为加入-max_hsps 1这个参数，结果被过滤掉了10行。那么过滤掉的是哪些行呢？
normal.txt第88、89、90行：
ST-E00126:289:H3727ALXX:5:1213:10835:15654 28 13 59 495 541 1.62e-15
ST-E00126:289:H3727ALXX:5:1213:10835:15654 28 176 203 661 688 5.87e-10
ST-E00126:289:H3727ALXX:5:1213:10835:15654 23 13 59 495 541 1.62e-15max-hsps-1.txt的88、89行：
ST-E00126:289:H3727ALXX:5:1213:10835:15654 28 13 59 495 541 1.62e-15
ST-E00126:289:H3727ALXX:5:1213:10835:15654 23 176 203 661 688 5.87e-102.4 结论
可见，normal.txt的第89行，因为query ID和subject ID都与第88行相同（只是匹配上的位置有点不同），被过滤掉了，那么我们就可以得出结论了：
由于同一个subject和同一个query可能在不同的位置上都有匹配，-max-hsps这个参数就是为了设置同一个query和同一个subject的匹配结果最多可以在结果文件中保留多少个。2.5 疑问
保留下来的是匹配度最好的吗？如果是的话，那么保留下来的这个的evalue一定得是最小的。2.6 验证
如果刚刚结论和疑问都成立的话，那么如下结果也是成立的：
normal.txt中，如果去掉query ID和subject ID都相同的行并且同时保留这些evalue最小的行，那么剩下的行应该和max-hsps-1.txt中是一样的。class BlastNRecord(object):
def __init__(self, qseqid, sseqid, qstart, qend, sstart, send, evalue):
self.qseqid = qseqid
self.sseqid = sseqid
self.qstart = qstart if isinstance(qstart, int) else int(qstart)
self.qend = qend if isinstance(qend, int) else int(qend)
self.sstart = sstart if isinstance(sstart, int) else int(sstart)
self.send = send if isinstance(send, int) else int(send)
self.evalue = evalue if isinstance(evalue, float) else float(evalue)
@classmethod
def parse_string(cls, string):
return cls(*string.rstrip().split('\t'))
def __eq__(self, other):
if isinstance(other, BlastNRecord):
if self.qseqid == other.qseqid and self.sseqid == other.sseqid:
return True
return False
def main(normal_file, max_hsps_1_file):
normal_dict = {}
max_hsps_1_dict = {}
with open(normal_file) as normal_f, open(max_hsps_1_file) as max_hsps_1_f:
# 把normal file读进字典中
for line in normal_f:
bn = BlastNRecord.parse_string(line)
if bn.qseqid in normal_dict:
if bn.sseqid in normal_dict[bn.qseqid]:
# 如果normal_dict存在和bn相同qseqid、sseqid则比较两者的evalue，保留较小的那个
normal_dict[bn.qseqid][bn.sseqid] = min([normal_dict[bn.qseqid][bn.sseqid], bn],
key=lambda _: _.evalue)
else:
# 如果normal_dict中存在该bn.qseqid，但是normal_dict[bn.qseqid]中不存在bn.sseqid
# 也添加进去
normal_dict[bn.qseqid][bn.sseqid] = bn
else:
# 如果normal_dict中不存在该bn，则添加到字典中去
normal_dict[bn.qseqid] = {bn.sseqid: bn}
# 把max hsps 1读进字典中
for line in max_hsps_1_f:
bn = BlastNRecord.parse_string(line)
if bn.qseqid in max_hsps_1_dict:
max_hsps_1_dict[bn.qseqid][bn.sseqid] = bn
else:
max_hsps_1_dict[bn.qseqid] = {bn.sseqid: bn}
# 第一步：检查normal.txt去重完成后是否只剩244行了
normal_cnt = 0
for qseqid in normal_dict:
normal_cnt += len(normal_dict[qseqid])
print('去除query ID和subject ID完全相同的行后，%s还剩下%s行。' % (normal_file, normal_cnt))
# 第二步：检查max-hsps-1.txt读入字典是不是也是有244行
# 如果是244行的话说明，max-hsps-1.txt中没有qseqid、sseqid完全相同的两行
max_hsps_1_cnt = 0
for qseqid in max_hsps_1_dict:
max_hsps_1_cnt += len(max_hsps_1_dict[qseqid])
print('读入%s后，还剩%s行。' % (max_hsps_1_file, max_hsps_1_cnt))
# 第三步，检查normal_dict是否和max_hsps_1_dict相同
for qseqid in normal_dict:
if qseqid in max_hsps_1_dict:
for sseqid in normal_dict[qseqid]:
if sseqid in max_hsps_1_dict[qseqid]:
if not normal_dict[qseqid][sseqid].evalue == max_hsps_1_dict[qseqid][sseqid].evalue:
print('发生错误，两个文件中的qseqid%s、sseqid%s的evalue不同' % (qseqid, sseqid))
else:
print('发生错误，%s中的qseqid%s、sseqid%s在%s中不存在' % (normal_file, qseqid, sseqid, max_hsps_1_file))
else:
print('发生错误，%s中的qseqid%s在%s中不存在' % (normal_file, qseqid, max_hsps_1_file))
sys.exit(2)
print('验证通过，%s去重后和%s完全相同' % (normal_file, max_hsps_1_file))
if __name__ == '__main__':
import sys
try:
main(sys.argv[1], sys.argv[2])
except (ValueError, IndexError):
print('使用方法: python %s normal-file max-hsps-1-file' % __file__)
sys.exit(1)运行：
$ python test_blastn.py ret/normal.txt ret/max-hsps-1.txt输出结果：
去除query ID和subject ID完全相同的行后，ret/normal.txt还剩下244行。
读入ret/max-hsps-1.txt后，还剩244行。
验证通过，ret/normal.txt去重后和ret/max-hsps-1.txt完全相同验证结果，结论是正确的，而且疑问也是正确的。