一句话加速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.233sgrep是處理匹配關系,獲得的是包含關鍵詞但不一定全等于關鍵詞,加一個-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.787shttps://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倍的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 获得诺贝尔奖之后影响力会下降?绘制精英科
- 下一篇: 诺奖奖金为何119年还没发完?