「JCVI」如何筛选得到最优blast比对结果?
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之前,就已经进行了输入文件的排序。
排序后文件的第一行,对应了最终的过滤文件 —— 也就是说后续等参数的过滤都是基于排序后的文件。
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
题外话
办公时间段摸鱼。
JCVI这个名字,到底怎么取出来的?是唐老师自己取的吗?
打开bing一搜,最上面一条的搜索结果引起了我的注意,
这名字,我是不是在哪里看到过???
打开bilibili —— 打开鬼谷老师的频道 —— 点开【科学八卦史】他以一人之力碾压全球科学家,让资本沦为自己的提款机,却这期节目。
哇,Amazing哇,Amazing哇,Amazing
原来是发明鸟枪测序法的爷 —— John Craig Venter
「JCVI」如何筛选得到最优blast比对结果?相关推荐
- 干掉「卧槽」,就用她了!
将开源小分队设为星标 精品文章第一时间读 大家好,我是爱撸码的开源大叔! 平时碰到厉害.愤怒.惊叹的人和事情,大家是不是只会一句卧槽行天下? 别问大叔怎么知道,俺也一样! 今天,大叔就给大家推荐一个非 ...
- 「JCVI教程」如何基于物种的CDS的blast结果绘制点图(dotplot)
这是唐海宝老师GitHub上的JCVI工具的非官方说明书. 该工具集的功能非常多,但是教程资料目前看起来并不多,因此为了能让更多人用上那么好用的工具,我就一边探索,一边写教程 这一篇文章教大家如何利用 ...
- 性能调优之三十六计 —— 「取而代之」Echo/Json 篇
文章目录 性能调优之三十六计 -- 「取而代之」Echo/Json 篇 Echo 高性能.极简框架 C.JSON Json-iterator Q&A 附录 性能调优之三十六计 -- 「取而代之 ...
- 深度学习在人脸识别中的应用 —— 优图祖母模型的「进化」
原作者: 腾讯优图 | 来自: 机器之心 序言--「弱弱」的人工智能 说到人工智能(Artificial Intelligence, AI)人们总是很容易和全知.全能这样的词联系起来.大量关于人工智能 ...
- 「走过」微软、优步,老工程师告诉你哪些数据结构和算法最重要
数据结构和基础算法作为计算机科学的必学课程,近几年却关注度越来越少.但程序员真的不再需要这两门基础知识了吗?一位在 Uber 等科技公司工作过的开发者分享了他的一手经验,告诉你实际工作中会用到哪些数据 ...
- 解密「天池」:如何做好一场万人AI竞赛的「大后方」?
来源:机器之心本文约9800字,建议阅读10+分钟面对数据集保护.算力公平性.结果可复现性等诸多挑战,天池是如何克服的呢? 一场一万五千人的竞赛,如何确保比赛顺利进行?如何保证公平公正?在这场活动中, ...
- 「机器学习」机器学习算法优缺点对比(汇总篇)
作者 | 杜博亚 来源 | 阿泽的学习笔记 「本文的目的,是务实.简洁地盘点一番当前机器学习算法」.文中内容结合了个人在查阅资料过程中收集到的前人总结,同时添加了部分自身总结,在这里,依据实际使用中的 ...
- 让你的系统“坚挺不倒”的最后一个大招——「降级」
来源:跨界架构师 前面两篇我们已经聊过了「熔断」(如何在到处是"雷"的系统中「明哲保身」?这是第一招)和「限流」(想通关「限流」?只要这一篇),这次我们聊的就是「高可用三剑客」中剩 ...
- 怎样对流媒体进行压力测试_对node工程进行压力测试与性能分析「干货」
作者:小黎 转发链接:https://mp.weixin.qq.com/s/WBe7ZLoqFD9UqNusnv_IDA 前言 在系统上线前,为了看下系统能承受多大的并发和并发下的负载情况,常常会先进 ...
最新文章
- 移动互联网改变商业环境:商品的颠覆
- MySQL模拟Oralce闪回操作
- Java EE 6 Web配置文件。 在云上。 简单。
- 程序员面试金典 - 面试题 17.07. 婴儿名字(并查集)
- java坐标移动题目case_用java怎样编写一个二维坐标平移程序
- 网页mp3提取器_用Python写一个酷狗音乐下载器!
- android chrome 不支持 audio/video的autoplay 属性
- Copy as Markdown - 将页面链接按照 Markdown 格式copy
- python自动备份手机_python实现自动备份windows应用数据
- csp2020 j2民间数据下载_摊开母婴市场数据集看一看
- 初中生学计算机应用有什么好方面,计算机有哪些专业 初中毕业学习相关专业有发展吗...
- python 坦克大战
- dcs world f15c教学_【温故知新】DCS如何操作?看这篇就全懂了!
- [MATLAB]关于SOR迭代计算其次线性方程组的数值解
- [译文]三重缓冲:为什么我们爱它
- [C语言]常用库函数
- H3C静态路由与BFD联动(单跳检测)配置案例
- 小圆圈o表示的数学符号是复合映射或Hadamard积(矩阵元素一一对应相乘)
- 用bim建模和用传统的图纸有什么差别?什么bim软件能提高建模效率?
- 云计算对电子商务行业的影响
热门文章
- OpenWrt路由开启DDNS+端口转发进行外网访问
- excel计算二元线性回归_用人话讲明白梯度下降Gradient Descent(以求解多元线性回归参数为例)...
- 怎么在添加为知笔记编辑器/为知笔记怎么用其他编辑器编辑/为知笔记怎么才能用Word/notepad++编辑
- 烤仔创作者联盟 | NFT是市场的下一个答案?或迎来新一轮“造福潮”
- Google TensorFlow课程 编程笔记(10)———使用神经网络对手写数字进行分类
- AI Conference:2018, 不容错过的世界人工智能大会 | 抢票
- 【IntelliJ IDEA插件】Alibaba Cloud AI Coding Assistant
- (SEED-Lab)Buffer Overflow Vulnerability Lab缓冲区溢出实验
- 最新2019版个税计算器(5000起征点 + 个税专项扣除项)
- 宝山区助行业强主体稳增长若干政策措施的实施细则(20条)(征求意见稿)