一句话加速grep近30倍
生物信息学习的正确姿势
NGS系列文章包括NGS基础、转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这)、ChIP-seq分析 (ChIP-seq基本分析流程)、单细胞测序分析 (重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程 (原理、代码和评述))、DNA甲基化分析、重测序分析、GEO数据挖掘(典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集)等内容。
最近做一个项目,需要从表达矩阵中提取单个特定基因的表达值。最开始时文件比较小,使用awk
单个读取处理也很快,但后来数据多了,从一个1.2 G
的文件中提取单个基因的表达需要30 s
,用grep
来写需要25 S
,这在平时写程序是可以接受的,但在网站上是接受不了的。所以就想着如何优化一下。
探索下来优化也很简单,把grep
换为LC_ALL=C grep
再加其它参数速度就快了近30倍,把时间控制在1 s
左右。
下面是整个探索过程 (写这篇总结文章是在早晨,服务器不繁忙,所以下面的示例中只能看出来快了5
倍左右。这也表明不加LC_All=C
时grep
受服务器负载影响较大,加了之后则几乎不受影响。)
获取单基因表达量
查看下文件大小
ls -sh 334d41a7-e34a-4bab-841c-eb07bd84513f.txt# 1.2G 334d41a7-e34a-4bab-841c-eb07bd84513f.txt
查看下文件内容
head 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | cut -f 1,2# Rnu2-1 -0.52
# Tmsb4Xp6 11.81
# S100A14 1.99
# Krt17 1.26
# Aldh1A1 6.92
# Fxyd3 0.56
# Rnu2-2P 0.35
# Rarres1 6.03
# Rnvu1-7 9.53
# Lcn2 3.44
假设基因名字大小写一致时使用awk
提取其表达信息,用时14 s
。
time awk '{if($1=="Tmsb4Xp6") print $2;}' 334d41a7-e34a-4bab-841c-eb07bd84513f.txt >1real 0m14.569s
user 0m12.943s
sys 0m0.626s
实际上大小写可能不一致而需要转换,耗时17 s
。
time awk '{if(tolower($1)=="tmsb4xp6") print $2;}' 334d41a7-e34a-4bab-841c-eb07bd84513f.txt >2real 0m17.638s
user 0m17.031s
sys 0m0.595s
采用grep
命令提取 (-i
忽略大小写),用时5 s
。
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | grep -i 'Tmsb4Xp6' >4real 0m5.454s
user 0m5.134s
sys 0m1.272s
上面的grep
是全句匹配,想着加上^
匹配行首是否会减少匹配量,速度能快一些,效果不明显,用时4 s
。
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | grep -iP '^Tmsb4Xp6' >5real 0m4.262s
user 0m3.984s
sys 0m1.233s
grep
是处理匹配关系,获得的是包含关键词但不一定全等于关键词,加一个-w
参数,匹配更精确些,耗时6.7 s
。
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | grep -iPw '^Tmsb4Xp6' >6real 0m6.723s
user 0m6.390s
sys 0m1.348s
从上面来看,采用正则限定并不能提速,还是采用固定字符串方式提取,速度也差不多,耗时5 s
。(fgrep
等同于grep -F
)
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | fgrep -i 'Tmsb4Xp6' >7real 0m5.496s
user 0m5.128s
sys 0m1.366s
主角出场,加上LC_ALL=C
后,速度明显提升了,只需要1 s
时间。
time LC_ALL=C fgrep -i '^Tmsb4Xp6' 334d41a7-e34a-4bab-841c-eb07bd84513f.txt >8real 0m1.027s
user 0m0.671s
sys 0m0.355s
多次测试下来,发现添加LC_ALL=C
后grep
命令快了很多,而且多次测试速度都很稳定 (不论服务器是繁忙还是空闲)。这里面的原理是涉及字符搜索空间的问题,我们操作的文件只包含字母、字符、数字,没有中文或其它复杂符号时都是适用的,具体原理和更多评估可查看文末的两篇参考链接,了解更多信息。
为了简化应用,我们可以alias grep='LC_ALL=C grep'
(把这句话放到~/.bashrc
或~/.bahs_profile
里面(具体用法见:PATH和path,傻傻分不清)),后续再使用grep
时就可以直接得到速度提升了。
time grep -F -i '^Tmsb4Xp6' 334d41a7-e34a-4bab-841c-eb07bd84513f.txtreal 0m1.013s
user 0m0.679s
sys 0m0.334s
那如果获取多个基因怎么操作呢?
一个方式是使用正则表达式,多个基因一起传递过去,分别匹配,耗时4.6 s
。
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | LC_ALL=C grep -iP 'Tmsb4Xp6|Sox1|Sox2|Sox3'real 0m4.654s
user 0m4.366s
sys 0m1.227s
或者还是使用固定字符串查找模式,把所有基因每行一个写入文件a
,然后再去匹配,耗时2.5 s
,且测试发现在基因数目少于10
时(这是通常的应用场景),基因多少影响不大 (这也说明能用固定字符串查找时最好显示指定)。
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | LC_ALL=C fgrep -i -f a >11real 0m2.539s
user 0m2.191s
sys 0m1.249s
这里还比较了另外2个号称比grep
快的命令ag
和rg
在这个应用场景没体现出性能优势。
time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | LC_ALL=C ag -i '^Tmsb4Xp6|Sox1|Sox2|Sox3' >10real 0m11.281s
user 0m9.713s
sys 0m5.326stime cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | rg -iF -f a >12real 0m4.337s
user 0m3.444s
sys 0m2.787s
https://www.inmotionhosting.com/support/website/speed-up-grep-searches-with-lc-all/
https://stackoverflow.com/questions/42239179/fastest-way-to-find-lines-of-a-file-from-another-larger-file-in-bash#
封面来源于:https://image.shutterstock.com/image-photo/conceptual-image-methaphor-busy-fast-260nw-641184331.jpg
最后来个促销吧,有缘者得(课程有效期到2025年,请忽略封面信息)
往期精品(点击图片直达文字对应教程)
后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集
一句话加速grep近30倍相关推荐
- 30倍加速,3毫秒极速识别,人、车、OCR等9大识别任务一网打尽
本文已在飞桨公众号发布,查看请戳链接: 30倍加速,3毫秒极速识别,人.车.OCR等9大识别任务一网打尽 人脸.车辆.人体属性.卡证.交通标识等经典图像识别能力,在我们当前数字化工作及生活中发挥着极其 ...
- java 下载加速_Java实现大文件下载,提速30倍!想学?我教你啊!
前言 兄弟们看到这个标题可能会觉得是个标题党,为了解决疑虑,我们先来看下最终的测试结果: 测试云盘下载的文件 46M,自己本地最大下载速度 2M 1. 单线程下载,总耗时: 603s img 2. 多 ...
- 互联网日报 | 蔚来汽车股价年内涨幅近十倍;小米成立可穿戴部;恒大物业拟12月2日港交所上市...
今日看点 ✦ 恒大物业拟12月2日港交所上市,至多募资158亿港元 ✦ 蔚来汽车股价年内涨幅近十倍,市值已超宝马.通用等传统汽车巨头 ✦ 小米最新组织调整:成立可穿戴部,高原为总经理 ✦ 中国移动发布 ...
- 调研发现,近30年来中国ICT产业发展关键是产业良好的营商环境
罗兰贝格与法国里昂商学院4月16日联合发布<中国ICT产业营商环境白皮书>,全面总结了中国ICT产业巨大的规模及高速增长态势,指出这是归功于持续改进的良好的营商环境,并进行了具体分析,且向 ...
- 上海外资银行总资产达1.5万亿元 较中国入世初期增近7倍
中新社北京1月24日电 (记者 王恩博)上海银保监局党委委员周文杰24日在北京透露,截至2018年末,该机构辖内外资银行总资产达1.5万亿元(人民币,下同),较2001年末中国加入世贸组织初期增长近7 ...
- IB网络用户数量超过私有网络近4倍
世界领先的数据中心端到端互连方案提供商Mellanox(纳斯达克交易所代码: MLNX)今天宣布,在最新发布的全球超级计算机Top500榜单中,InfiniBand再次延续了其在互连方案上的绝对领先地 ...
- 14.7倍推理加速、18.9倍存储节省!北航、商汤、UCSD提出首个点云二值网络 | ICLR 2021...
允中 编辑整理 量子位 报道 | 公众号 QbitAI 编者按: 无论是在自动驾驶场景中,还是在手持移动设备上,基于点云的深度学习模型应用越来越广泛. 但这些离线边缘场景自身的限制,给模型的推理.存储 ...
- Ian Thiel:靠这 3 点,实现 30 倍增长,从不盈利到营收 5.5 亿
以全球视野萃取本土智慧.继精益数据分析鼻祖 Alistair Croll 做客神策 2017 数据驱动大会后,被誉为"所到之处必增长"的硅谷增长专家 lan Thiel 的国内首次 ...
- 预览速度提升30倍,这是什么黑科技?(天猫618之3D渲染篇)
简介: 天猫618宣布的 3D 购物时代,相信有很多小伙伴好奇,这背后有哪些"黑科技"?橙子从以下三点为你揭秘--3D实景复刻.3D渲染.3D算法,上周讲了<天猫618宣布开 ...
最新文章
- IDE---Gvim支持php的函数自动补全
- jsp项目在idea需要导入什么依赖_Java开发工具IntelliJ IDEA配置项目系列教程(五):模块依赖...
- 鸡蛋中营养和脂质含量与降低LDL的食物
- 使用trackBy启动流程
- 铁拳nat映射_铁拳如何重塑我的数据可视化设计流程
- 协方差意味着什么_“零”到底意味着什么?
- Java并发编程实战~Guarded Suspension模式
- 设计师对孟菲斯设计风还不了解?
- javascript中五种常见的DOM方法
- 第四十二节,configparser特定格式的ini配置文件模块
- 理解 Python 中的线程
- iOS开发-iOS学习完整路线
- 清华大学计算机杜瑜皓,立足改革、多元择优为清华选拔有理想有才能的优秀学子-清华大学新闻网...
- java traingdx函数实现_提取伪彩色图像的信息
- TypeError: Cannot set property ‘styles‘ of undefined
- Word中插入参考文献 自动管理
- Codeforces Round #507 (Div. 2, based on Olympiad of Metropolises) B. Shashlik Cooking
- AD(十九)class、设计参数、规则的创建
- android开发中磁场传感器,Android开发获取传感器数据的方法示例【加速度传感器,磁场传感器,光线传感器,方向传感器】...
- 二叉搜索树(BST)——基本概念及基本实现代码