「JCVI」如何筛选得到最优blast比对结果？

2022-07-16 本文已影响0人陈有朴

JCVI，包含了太多的功能，但是我感觉好像又没有一个特别好的说明文档（小声bb，感谢唐老师的开发的好用工具）

blast比对

未过滤的blast比对结果，所使用参数是：-outfmt 6 -max_hsps 1 -max_target_seqs 500 -num_threads 16 -evalue 1e-5 -perc_identity 70 -task megablast

Note：-max_target_seqs 500为默认参数，代表query序列最多可以保留500条不同ref序列的比对结果。

查看哪些query序列有多个比对结果

这个标题是句废话，因为我已经设置了-max_target_seqs 500。

但是由于我构建的数据库只有1800左右的序列，输入序列有1700条，因此每一条query序列也不太可能出现500条这么多的结果。

# 输出有多个比对结果query序列ID
awk '{print $1}' raw.blast.txt | sort | uniq -d

现在想要得到最优的blast比对结果，我应该怎么操作？

自行编写脚本，根据identity、e-value、score、pairing length等参数（×）
使用JCVI（√）

过滤blast比对结果

使用如下命令，就可以得到最优blast比对结果，

python -m jcvi.formats.blast best -n 1 raw.blast.txt

JCVI过滤的标准是什么？

（1）解读源码

Produce a new blast file and filter based on:
- score: >= cutoff
- pctid: >= cutoff
- hitlen: >= cutoff
- evalue: <= cutoff
- ids: valid ids

也就是说，是基于1）比对得分、2）identity、3）联配长度、4）e-value以及自行定义的5）序列ID来得到过滤后的blast文件的。

（2）用几条序列测试一下

根据源码，已经可以得到JCVI是根据1）

那当e-value和比对得分以及比对上的长度都是相等的时候？JCVI是怎么提取对应结果的？

运气很好的是，我试了上一个部分输出的几个有多个比对结果的query ID，就找到了可以用于解析JCVI如何过滤blast比对结果的query序列。

将对应query序列的结果grep出来，保存下列文件，命名为query1.blast.txt

query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 219 767 0.0 640
query1 Os01t0624000-04-E6 87.796 549 64 1 353 898 276 824 0.0 640
query1 Os01t0624000-02-E10 87.796 549 64 1 353 898 292 840 0.0 640

运行JCVI，python -m jcvi.formats.blast best -n 1 query1.blast.txt（结果文件为query1.blast.txt.best）

query1  Os01t0624000-01-E10     87.80   549     64      1       353     898     219     767     0       640

然后我的疑问出现了：为什么3条序列e-value和Score都一样？为什么JCVI选择了第一条？

Note：上述3种比对结果的联配长度是一样的，均为548

与比对到参考序列上的起始位置有关，即选择起始位置最前的比对结果
与参考序列的ID有关

测试1：起始位置

文件修改为如下，保存为test1.blast.txt

query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 219 767 0.0 640
query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 192 740 0.0 640

运行JCVI，python -m jcvi.formats.blast best -n 1 test1.blast.txt，查看该命令的结果文件：

query1  Os01t0624000-01-E10     87.80   549     64      1       353     898     192     740     0       640

需要注意的是，终端输出的命令为：sort -k1,1 -k12,12gr test2.blast.txt -o test1.blast.txt，

再使用JCVI之前，就已经进行了输入文件的排序。

排序后文件的第一行，对应了最终的过滤文件 —— 也就是说<font color='yellow'>后续等参数的过滤都是基于排序后的文件</font>。
query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 192 740 0.0 640
query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 219 767 0.0 640

测试2：序列ID

文件修改为如下格，保存为test2.blast.txt，

query1 Os01t0624000-04-E6 87.796 549 64 1 353 898 219 767 0.0 640
query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 219 767 0.0 640

使用JCVI过滤得到的重排后的文件以及结果如下，

# sorted file
query1 Os01t0624000-01-E10 87.796 549 64 1 353 898 219 767 0.0 640
query1 Os01t0624000-04-E6 87.796 549 64 1 353 898 219 767 0.0 640
# filter result
query1  Os01t0624000-01-E10     87.80   549     64      1       353     898     219     767     0       640