生物信息學習的正確姿勢

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倍的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。