一句话加速grep近30倍-CSDN博客

本文链接：https://blog.csdn.net/qazplm12_3/article/details/107829342

生物信息学习的正确姿势

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 >1

real    0m14.569s
user    0m12.943s
sys    0m0.626s

实际上大小写可能不一致而需要转换，耗时17 s。

time awk '{if(tolower($1)=="tmsb4xp6") print $2;}' 334d41a7-e34a-4bab-841c-eb07bd84513f.txt >2

real    0m17.638s
user    0m17.031s
sys    0m0.595s

采用grep命令提取 (-i忽略大小写)，用时5 s。

time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | grep -i 'Tmsb4Xp6' >4

real    0m5.454s
user    0m5.134s
sys    0m1.272s

上面的grep是全句匹配，想着加上^匹配行首是否会减少匹配量，速度能快一些，效果不明显，用时4 s。

time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | grep -iP '^Tmsb4Xp6' >5

real    0m4.262s
user    0m3.984s
sys    0m1.233s

grep是处理匹配关系，获得的是包含关键词但不一定全等于关键词，加一个-w参数，匹配更精确些，耗时6.7 s。

time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | grep -iPw '^Tmsb4Xp6' >6

real    0m6.723s
user    0m6.390s
sys    0m1.348s

从上面来看，采用正则限定并不能提速，还是采用固定字符串方式提取，速度也差不多，耗时5 s。(fgrep等同于grep -F)

time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | fgrep -i 'Tmsb4Xp6' >7

real    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 >8

real    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.txt

real    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 >11

real    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' >10

real    0m11.281s
user    0m9.713s
sys    0m5.326s

time cat 334d41a7-e34a-4bab-841c-eb07bd84513f.txt | rg -iF -f a >12

real    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年，请忽略封面信息）