说到底用R语言data.table包举办拍卖,awk举办的是逐行处理

   
由于基因组数据过大,想进一步用R语言处理担心系统内存不够,由此想着将文件按染色体拆分,发现python,awk,R
语言都可以相当简单疾速的兑现,那么速度是不是有差距啊,因而在跑几个50G的大文件从前,先用了244MB的多少对一一脚本举行测试,并且将其速度举行相比较。

   
由于基因组数据过大,想进一步用R语言处理担心系统内存不够,由此想着将文件按染色体拆分,发现python,awk,R
语言都可以分外简单快捷的兑现,那么速度是不是有差异啊,由此在跑多少个50G的大文件往日,先用了244MB的数目对一一脚本举行测试,并且将其速度举行相比。

文本处理是生物音讯学中的一个最紧要步骤,但我一般是不太关心数据结构,也不关注使用什么办法什么语言,只要能将文件中的音信快狠准的领到出来就行了,毕竟精力有限,数据的生物学意义远比数据提取方法紧要。在文本处理中我个人的习惯是什么便宜用咋样,比如不规则文本提取用grep/awk正则匹配提取,规则且量不大的数字数据表用R,大量系列数据用python,啥福利啥来。原则有六个,代码越短小越好,提取速度越快越好,可想而知就是追求短小快。个人觉得在海洋生物音讯学里,没有最好的言语,只有最合适的语言,各样脚本语言都得会那么一些。

第一是awk处理,awk举办的是逐行处理,具有温馨的语法,具有很大的油滑,一行代码解决,用时24S,

第一是awk处理,awk举行的是逐行处理,具有自己的语法,具有很大的八面玲珑,一行代码解决,用时24S,

1. awk大法好

我个人最欣赏awk的某些就是不需显式读入文件,而perl/python/R/java,都需要读入文件然后赋值给变量,还得声明变量类型,处理好了还得写个出口,麻烦的很。但awk
不需要啊,数据可以来自专业输入(stdin)、一个或两个公文,或任何命令的出口,完事后直接”>”输出到文件,省事。再搭配awk内置的变量/函数数组,以及与shell其他命令搭配,简直没有解决不了的文本处理问题。

 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 #!/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
下面是个处理bed格式的实例:
# 首先安装bedtools为的是获得随机bed文件:
sudo apt-get install bedtools
# 使用bedtools的randomBed获得1000条随机bed 位置:
randomBed -seed 2 -n 1000 -l 400 -g ./hg19.chrom_24.sizes >./a.bed

说明:hg19.chrom_24.sizes是基因组size文件,是采纳UCSC
的fetchChromSizes(从此间下载:http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86\_64/
程序下载或者自己编写也得以,格式如下:

ca88亚洲城网站 1

#查看随机序列:
head a.bed

说明:前三列是岗位信息,第四列是命名,那里是遵照行号来定名的,第五列是长度,最终一列是正负链。

ca88亚洲城网站 2

下边举办操作,首先是排序:

# 使用sort首先按照正负链,其次是染色体,最后是位置信息排序:
sort -k 6,6 -k 1,1 -k 2n,2 ./a.bed -o ./a.bed

说明:sort
函数有无数参数,-k是按域排序,先按照第六列链排序,然后是染色体第一列,然后是岗位信息第二列,第二列是以数字排序,所以要加n。结果如下:

ca88亚洲城网站 3

利用awk 提取指定染色体消息:

# 只提取chr1的信息:
awk 'BEGIN{FS="\t";OFS="\t"}{if($1=="chr1")print }' a.bed|tail

说明:BEGIN是文件处理前的操作,FS表示定义文件分割,这里是制表符,OFS是概念输出的分割符;$1代表第一列数据。

# 只提取chr1正链的信息:
awk 'BEGIN{FS="\t";OFS="\t"}{if($1=="chr1" && $6=="+")print }' a.bed|tail

说明:多加个规格就行了。


上边简单的领到功用awk
就是单排实现,而其它的如perl/python脚本需要或多或少行才能落实。后边介绍对gtf文件的拍卖以及多文文件比对处理,perl/python处理要复杂很多,而awk只需一行。
更多原创美观内容敬请关注生信随想

ca88亚洲城网站 4

 

 

ca88亚洲城网站 5

ca88亚洲城网站 6

 

 

接下去是用python,python语言简练,书写方便。由此连忙就落实了先后,同样逐行处理,比awk添加了一些细节,只挑出需要的染色体。用时19.9秒。

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

 

 

ca88亚洲城网站 7

ca88亚洲城网站 8

 

 

 

 

终极用R语言data.table包举办拍卖,data.table是data.frame的高档版,在进度上作了很大的改正,不过和awk和python相相比较,具有优势呢?

最后用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 }
 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 }

ca88亚洲城网站 9

ca88亚洲城网站 10

   
用时10.6秒,发现刚刚读完数据,登时就处理和写出了事,处理和写出时间十分短,由此全体用时较短。

   
用时10.6秒,发现刚刚读完数据,立时就处理和写出得了,处理和写出时间非凡短,由此全体用时较短。

 

ca88亚洲城网站, 

总结

总结

    虽然都是逐行处理,但由上述结果算计awk内部运转并不曾python快,但awk书写一行代码搞定,书写速度快,至于python比data.table慢,估算原因是R
data.table用C语言写,并且应用多线程写出,hash读取,传地址各类措施优化速度的结果。当然,上述结果仅供参考。

    即便都是逐行处理,但由上述结果估算awk内部运转并没有python快,但awk书写一行代码搞定,书写速度快,至于python比data.table慢,推测原因是R
data.table用C语言写,并且使用多线程写出,hash读取,传地址各样法子优化速度的结果。当然,上述结果仅供参考。

 

 

相关文章