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

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

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

#!/usr/bin/sh
function main()
{
start_tm=date
start_h=`$start_tm +%H`
start_m=`$start_tm +%M`
start_s=`$start_tm +%S`
awk -F $sep '{print $1","$2","$3 >> "'"$inputfile"'""_"$1}' $inputfile
end_tm=date
end_h=`$end_tm +%H`
end_m=`$end_tm +%M`
end_s=`$end_tm +%S`
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"}'`
echo "Finished in "$use_tm
}


if [ $# == 2 ]; then
sep=$1
inputfile=$2
main
else
echo "usage: SplitChr.sh sep inputfile"
echo "eg: SplitChr.sh , test.csv"
fi

 

 

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

#!/usr/bin/python
import sys
import time
def main():
    if len(sys.argv)!=3:
        print "usage : SplitChr sep inputfile eg: SplitChr ',' test.txt"
        exit()
    sep=sys.argv[1]
    filename=sys.argv[2]
    f=open(filename,'r')
    header=f.readline()
    if len(header.split(sep))<2:
        print "The sep can't be recongnized !"
        exit()
    chrLst=range(1,23)
    chrLst.extend(["X","Y"])
    chrLst=["chr"+str(i) for i in chrLst]
    outputdic={}
    for chrI in chrLst:
        output=filename+"_"+chrI
        outputdic[chrI]=open(output,'w')
        outputdic[chrI].write(header)
    for eachline in f:
        tmpLst=eachline.strip().split(sep)
        tmpChr=tmpLst[0]
        if tmpChr in chrLst:
            outputdic[tmpChr].write(eachline)
    end=time.clock()
    print "read: %f s" % (end - start)



if __name__=='__main__':
    start=time.clock()
    main()

 

 

 

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

#!/usr/bin/Rscript
library(data.table)
main <- function(filename,sep){
started.at <- proc.time()
arg <- commandArgs(T)
sep <- arg[1]
inputfile <- arg[2]
dt <- fread(filename,sep=sep,header=T)
chrLst <- lapply(c(1:22,"X","Y"),function(x)paste("chr",x,sep=""))
for (chrI in chrLst){
    outputfile <- paste(filename,"_",chrI,sep="")
    fwrite(dt[.(chrI),,on=.(chr)],file=outputfile,sep=sep)
}
cat ("Finished in",timetaken(started.at),"\n")
}

arg <- commandArgs(T)
if (length(arg)==2){
sep <- arg[1]
filename <- arg[2]
main(filename,sep)
}else{
cat("usage: SplitChr.R sep inputfile eg: SplitChr.R '\\t' test.csv","\n")
}

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

 

總結

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