都亟待读入文件然后赋值给变量,最后用RAV4语言data.table包举行拍卖

   
由于基因组数据过大,想进一步用奥迪Q3语言管理缅怀系统内部存款和储蓄器非常不足,因而想着将文件按染色体拆分,开掘python,awk,奥迪Q3语言都能够特别简单快速的兑现,那么速度是否有差别啊,因而在跑多少个50G的大文件在此之前,先用了244MB的多少对各种脚本进行测验,况兼将其速度实行自己检查自纠。

文本管理是生物音信学中的三个生死攸关步骤,但本人一般是不太关爱数据结构,也不关怀使用什么措施什么语言,只要能将文件中的新闻快狠准的提抽出来就行了,终究精力有限,数据的生物学意义远比数据提取方法主要。在文本处理中本人个人的习贯是怎么实惠用哪些,比如不规则文本提取用grep/awk正则相配提取,法规且量比较小的数字数据表用翼虎,大批量队列数据用python,啥福利啥来。原则有多个,代码越短小越好,提取速度越快越好,简单的讲就是追求短小快。个人认为在海洋生物音讯学里,未有最佳的言语,独有最合适的语言,种种脚本语言都得会那么一些。

先是是awk管理,awk实行的是逐行管理,具备温馨的语法,具备十分大的灵活性,一行代码消除,用时24S,

1. awk大法好

本人个人最欣赏awk的一些正是不需显式读入文件,而perl/python/PAJERO/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
上面是个管理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

说明:ca88亚洲城网站,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

 

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

 

ca88亚洲城网站 6

 

 

末尾用福特Explorer语言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 }

ca88亚洲城网站 7

   
用时10.6秒,开采刚刚读完数据,立时就管理和写出扫尾,管理和写出时间特别短,由此总体用时异常的短。

 

总结

    纵然都是逐行管理,但由上述结果猜度awk内部运行并未python快,但awk书写一行代码化解,书写速度快,至于python比data.table慢,猜度原因是PRADOdata.table用C语言写,并且利用八线程写出,hash读取,传地址各类措施优化速度的结果。当然,上述结果仅供参照他事他说加以考察。

 

相关文章