將基因組數據分類並寫出文件，python，awk，R data.table速度PK

由於基因組數據過大，想進一步用R語言處理擔心繫統內存不夠，因此想著將文件按染色體拆分，發現python，awk，R 語言都能夠非常簡單快捷的實現，那麼速度是否有差距呢，因此在跑幾個50G的大文件之前，先用了244MB的數據對各個腳本進行測試，並且將其速度進行對比。

首先是awk處理，awk進行的是逐行處理，具有自己的語法，具有很大的靈活性，一行代碼解決，用時24S，

 1 #!/usr/bin/sh
 2 function main
 3 {
 4 start_tm=date
 5 start_h=`$start_tm +%H`
 6 start_m=`$start_tm +%M`
 7 start_s=`$start_tm +%S`
 8 awk -F $sep '{print $1","$2","$3 >> "'"$inputfile"'""_"$1}' $inputfile
 9 end_tm=date
10 end_h=`$end_tm +%H`
11 end_m=`$end_tm +%M`
12 end_s=`$end_tm +%S`
13 use_tm=`echo $end_h $start_h $end_m $start_m $end_s $start_s | awk '{ print ($1 - $2),"h",($3-$4),"m",($5-$6),"s"}'`
14 echo "Finished in "$use_tm
15 }
16 
17 
18 if [ $# == 2 ]; then
19 sep=$1
20 inputfile=$2
21 main
22 else
23 echo "usage: SplitChr.sh sep inputfile"
24 echo "eg: SplitChr.sh , test.csv"
25 fi

接下來是用python,python語言簡單，書寫方便。因此很快就實現了程序，同樣逐行處理，比awk添加了一點細節，只挑出需要的染色體。用時19.9秒。

 1 #!/usr/bin/python
 2 import sys
 3 import time
 4 def main:
 5     if len(sys.argv)!=3:
 6         print "usage : SplitChr sep inputfile eg: SplitChr ',' test.txt"
 7         exit
 8     sep=sys.argv[1]
 9     filename=sys.argv[2]
10     f=open(filename,'r')
11     header=f.readline
12     if len(header.split(sep))<2:
13         print "The sep can't be recongnized !"
14         exit
15     chrLst=range(1,23)
16     chrLst.extend(["X","Y"])
17     chrLst=["chr"+str(i) for i in chrLst]
18     outputdic={}
19     for chrI in chrLst:
20         output=filename+"_"+chrI
21         outputdic[chrI]=open(output,'w')
22         outputdic[chrI].write(header)
23     for eachline in f:
24         tmpLst=eachline.strip.split(sep)
25         tmpChr=tmpLst[0]
26         if tmpChr in chrLst:
27  outputdic[tmpChr].write(eachline)
28     end=time.clock
29     print "read: %f s" % (end - start)
30 
31 
32 
33 if __name__=='__main__':
34     start=time.clock
35     main

最後用R語言data.table包進行處理，data.table是data.frame的高級版，在速度上作了很大的改進，但是和awk和python相比，具有優勢嗎?

 1 #!/usr/bin/Rscript
 2 library(data.table)
 3 main <- function(filename,sep){
 4 started.at <- proc.time
 5 arg <- commandArgs(T)
 6 sep <- arg[1]
 7 inputfile <- arg[2]
 8 dt <- fread(filename,sep=sep,header=T)
 9 chrLst <- lapply(c(1:22,"X","Y"),function(x)paste("chr",x,sep=""))
10 for (chrI in chrLst){
11     outputfile <- paste(filename,"_",chrI,sep="")
12     fwrite(dt[.(chrI),,on=.(chr)],file=outputfile,sep=sep)
13 }
14 cat ("Finished in",timetaken(started.at),"\n")
15 }
16 
17 arg <- commandArgs(T)
18 if (length(arg)==2){
19 sep <- arg[1]
20 filename <- arg[2]
21 main(filename,sep)
22 }else{
23 cat("usage: SplitChr.R sep inputfile eg: SplitChr.R '\\t' test.csv","\n")
24 }

用時10.6秒，發現剛剛讀完數據，立刻就處理和寫出完畢，處理和寫出時間非常短，因此總體用時較短。

總結

雖然都是逐行處理，但由上述結果猜測awk內部運行並沒有python快，但awk書寫一行代碼搞定，書寫速度快，至於python比data.table慢，猜測原因是R data.table用C語言寫，並且運用多線程寫出，hash讀取，傳地址各種方式優化速度的結果。當然，上述結果僅供參考。