将基因组数据分类并写出文件,python,awk,途锐 data.table速度PK

   
由于基因组数据过大,想进一步用奥德赛语言处理担心系统内存不够,因而想着将文件按染色体拆分,发现python,awk,Sportage语言都可以卓殊容易快速的兑现,那么速度是或不是有反差啊,由此在跑多少个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

 

图片 1

 

接下去是用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()

 

图片 2

 

 

末尾用奥迪Q5语言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 }

图片 3

   
用时10.6秒,发现刚刚读完数据,立刻就处理和写出了事,处理和写出时间相当长,因而全体用时较短。

 

总结

    即使如此都以逐行处理,但由上述结果揣度awk内部运营并没有python快,但awk书写一行代码化解,书写速度快,至于python比data.table慢,估计原因是Muranodata.table用C语言写,并且使用二十四线程写出,hash读取,传地址种种法子优化速度的结果。当然,上述结果仅供参考。