将基因组数据分类并写出文件

新葡萄娱乐,   
由于基因组数据过大,想进一步用Odyssey语言管理担忧系统内部存款和储蓄器缺乏,因而想着将文件按染色体拆分,发掘python,awk,大切诺基语言都能够非常简单飞速的贯彻,那么速度是还是不是有反差啊,由此在跑多少个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

 

 

末尾用路虎极光语言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慢,估量原因是Highlanderdata.table用C语言写,而且选拔多线程写出,hash读取,传地址种种方法优化速度的结果。当然,上述结果仅供参谋。