最近我被老板的绝招打败了,心服口服。
有感之余,决定总结一把自己在用的,一级刚刚学到手的各类效率秘籍,算是个总结笔记吧。如果对其他人有帮助,那本鼠就更开心了。同时欢迎各位来本博客做做,磕磕瓜子,交流交流生物信息的各种心得,事实证明:生物信息跟拿基因枪一样,都是手艺活嘛.
我只臭屁的希望若有哥们觉得有点用,然后转载时,请注明出处:即科学网本鼠的博问链接,哦哈哈哈.
1. hash的搜索效率
先说说老板如何搞定我的吧,就是这个hash..
其实早在N年前(本科)啃小骆驼的时候,就读到过hash的搜索是经过优化的,本质类似于一个二维数组。潜意识里,这就给了我一个误区,还是个二维数组嘛,快也就是50%效率之类的吧。事实证明我想当然了。前几日需要整理一个project的结果,就写了类似如下的perl脚本(摘录)。
#re_structure the FKPM file
my @fkpm;
for(my $i=0;$i<scalar(@array);$i++){
my @line=split(/\s+/,$array[$i]);
$fkpm[$i][0]=$line[0];#id
$fkpm[$i][1]=$line[6];#FKPM
$fkpm[$i][2]=$line[9];#cover rate
}
my $id;
foreach $a(@name){ #第一个loop
if($a =~ /(Th.*?)\s+(.*?)$/){
$id=$1;
}
my $mark=0;
for(my $i=0;$i<scalar(@fkpm);$i++){ #第二个loop
if($id eq $fkpm[$i][0]){
print "$a\t$fkpm[$i][1]\t$fkpm[$i][2]\n";
$mark=1;
last;
}
}
if($mark == 0){
print "$a\t\-\t\-\n";
}
}
写完后估计了一下时间,于是跟老板说,估计得跑个把小时,吃饭去了。老板说,用hash最多就是一分钟的事儿,我还不服气的辩解了几句,找了几个借口(此处省略),在老板提高嗓门后觉得:呃,应该是自己错了。
吃饭大概半小时,回来后果然:这段脚本只处理了三成的数据。于是乎修改成了hash,又跑了一遍,用了15秒左右.....程序就完成了...15秒左右!!我囧了...以下是我修改的代码:
#re_structure the FKPM file
my %fkpm;
for(my $i=0;$i<scalar(@array);$i++){
my @line=split(/\s+/,$array[$i]);
$fkpm{$line[3]}="$line[6]\t$line[9]";
#print "$fkpm[$i][0] $fkpm[$i][1] $fkpm[$i][2]\n";die;
}
my $id;
foreach $a(@name){ #第一个loop
if($a =~ /(Th.*?)\s+(.*?)$/){
$id=$1;
}
if(!$fkpm{$id}){
print "$a\t\-\t\-\n";
}
else{
print "$a\t$fkpm{$id}\n"; #hash
}
}
还是知其然,再知其所以然吧。我再也不敢像原来那样囫囵吞枣了...
hash是如何做到快速搜索定位的呢?查查定义。
若结构中存在关键字和K 相等的记录,则必定在f (K )的存储位置上。由此,不需比较便可直接取得所查记录。称这个对应关系f 为散列函数 (Hash function),按这个思想建立的表为散列表。 好嘛,人家压根就不找存储地址,直接根据key把地址算出来,怪不得。
2。善用shell命令
我一度做东西有误区,就是一到文本操作就用perl,哗啦哗啦写一大堆,耽误时间,效率低下。
不久前才改过来。以下是个例子。
我需要把一个fasta文件的部分信息提出来,其文件格式如下:
>Th0000001 old_id=63_1_CFGU_CFGW_CFHG_CFHH_CFHI
CAAACAACATGCTAAGATCTGACCTTGTTCATCTCCACATCTTTTCAACGACAGATTAAA
ACTGAAAA.....
>Th0000001 old_id=.....
....
....
想要获得的信息为:
Th0000001 63_1_CFGU_CFGW_CFHG_CFHH_CFHI
Th0000002
半年前我会写个perl,估计十几分钟就没了,再出个小错,说不定半个小时就搭进去了. 其实只需要在shell里打了一行命令而已,耗时可以用秒计算吧:
grep ">" Th_unigene_assembled_ESTs.fasta |cut -c 2-|sed "s/\told_id\=/\t/" >name_list
# Th_unigene_assembled_ESTs.fasta 源文件
# name_list 输出文件
第二个例子
有时候一个程序要对不同的数据集都跑一遍,这些命令行有一个特点:非常相似!
难道就为了这,在perl里面加个loop,那也太...那..那怎么办?
唉,我曾经不胜其烦的按着上下键(调出历史命令)修改一点儿,运行。
一遍还好,有时候上下游的数据出了问题,需要程序重跑就无语了,以下YY情节:
某同学:耗子兄,你这个EST放错行了嘛,我们没法用。
本鼠:哦,好,就改一行代码嘛,容易。
一刻钟后....
某同学:那个,改好了么?我们老板要看呢
本鼠:....
旁白:一只耗子正在电脑前勤奋的按着上下键和回车....
其实嘛,磨刀不误砍柴工,写个shell脚本就好了。不过请注意,这个脚本是可以不用手写的 ,是可以再写段代码生成的..某些大侠一眼就能看出下文可以用for loop生成...
#!/usr/bin/bash
perl add_FRAG_hash.pl est.list 0_best_blast/1.blastn 5_count_frag_for_each_dataset/1.2.frag.count >1_added_frag
perl add_FRAG_hash.pl 1_added_frag 0_best_blast/2.blastn 5_count_frag_for_each_dataset/2.2.frag.count >2_added_frag
perl add_FRAG_hash.pl 2_added_frag 0_best_blast/3.blastn 5_count_frag_for_each_dataset/3.2.frag.count >3_added_frag
perl add_FRAG_hash.pl 3_added_frag 0_best_blast/4.blastn 5_count_frag_for_each_dataset/4.2.frag.count >4_added_frag
perl add_FRAG_hash.pl 4_added_frag 0_best_blast/5.blastn 5_count_frag_for_each_dataset/5.2.frag.count >frag_added.list
rm -f ?_added_frag#删除中间文件
什么什么?才几行还不如直接敲命令?重跑也认了?!
那要是,几百行呢?
此处罗列一个本鼠用来生成shell脚本的命令,不解释了..只是说一下,这个 run_header_remove.sh一共882行..
ls all_out/ | sed 's/^.*$/perl header_remove.pl \<all_out\/& >clean_header_psl\/&/' >run_header_remove.sh
3. $$的妙用
感谢现在发达的科技,包括我们实验室在内的大多数服务器升级到了N核(N >= 2)..不怀好意的想占用多个CPU的同学远不止本鼠。八仙过海,各显神通,有用perl的多线程包的,有喜欢shell下直接一个 “&” 省事儿的,爽啊,原来需要1个CPU跑7天的工作,现在就7个CPU跑1天好了..美哉,美哉..
然而这对于本鼠有1个bug,那就是:要是写的程序早在运行中存取临时文件怎么办...
7个CPU,7个线程,跑的是不同的数据,要是线程1拿了线程2的中间文件,线程2找不到,顺手牵走了线程5的,导致线程7找不到文件报错...这个...呃...
小A说,那就创建7个工作文件夹,在里面分别把程序跑一遍。
小C说,不就是中间文件的文件名嘛,我把程序修改下,弄个程序1,程序2...
这两种想法直观,却不符合生物信息的精神,小A和小C勤奋的跑出结果,牛了,也只是NA和NC...什么什么,我说大了,好吧,至少不符合perl的编程精神。那perl的编程精神是什么呢?七个字:随心所欲的懒人.....囧
懒人们只需要在文件名后面加上$$就好了,它有个响当当的大名:Perl解释器的进程ID(PID). $$是perl内置变量,只要运行perl就自动赋值,那么中间文件temp就会变成temp12345, temp44857...一般的服务器的PID至少要一天才能转一圈吧(99999个呢)基本保险了.
照例再贴出本鼠的一个小脚本做个例子:
my $input_file=$ARGV[0];
`grep ">" $input_file |sed "s/\>//" >temp$$`; #$$
open FILE, "<temp$$" or die"hell:$!"; #读入临时文件$$
my @array=<FILE>;
close FILE;
chomp @array;
`rm temp$$`; #$$哥们你可以功成身退了(删除)
my %frag_count;
foreach $a(@array){
my @line = split(/\./,$a);
$frag_count{$line[0]}=$line[1];
}
my @array2=keys(%frag_count);
for(my $i=0;$i<scalar(@array2);$i++){
print "$array2[$i]\t$frag_count{$array2[$i]}\n";
}
PS: 如果同时服务器上有个程序不断放出子线程,几个小时把PID用一圈的话..呃,相信同学们的智慧是无穷的。
友情提示:换个系统时间变量用用吧,看100年内它能转一圈不?
4.注意多loop下的代码行数
事实证明一个细节会被大量重复来放大,良好的编程习惯很重要.
这个就直接上例子了,结果是提升了约4倍的效率.
#get the line numbers for calculation my $num1=scalar(@native_one); my $num2=scalar(@native_two); my @residue=residue(\@native_one,\@native_two,$num1,$num2,$range); return @residue; }
#calculate the distance and find the residues sub residue{ my($first,$second,$num1,$num2,$range)=@_; my $distance; my @result_first; my @result_second; for(my $i=0;$i<$num1;$i++){ for(my $j=0;$j<$num2;$j++){ #双loop $distance=distance($$first[$i][2],$$first[$i][3],$$first[$i][4],$$second[$j][2],$$second[$j][3],$$second[$j][4]);
if($distance <= $range){ my $temp="$$first[$i][0]"; push @result_first,$temp; my $temp2="$$second[$j][0]"; push @result_second,$temp2; } } } #delete the reduntant residues my %hash=(); my @result1 = grep{$hash{$_}++ <1} @result_first; my %hash2=(); my @result2 = grep{$hash2{$_}++ <1} @result_second; my @result_final=(@result1,"divide",@result2); return @result_final; } #Function for distance sub distance{ #双loop下调用的子函数 1 my($x1,$y1,$z1,$x2,$y2,$z2)=@_; my $square=($x1-$x2)**2+($y1-$y2)**2+($z1-$z2)**2; my $result=sqrt($square); return $result; }
修改子函数为:
sub distance{ #双loop下调用的子函数 2 my $result=sqrt( ($_[0]- $_[3])**2 + $_[1]- $_[4])**2 + $_[3]- $_[5])**2 );
}
行了,就到这儿吧,天色已晚,本鼠该回窝了..各位看官,来磕瓜子啊..
转载本文请联系原作者获取授权,同时请注明本文来自陈昊科学网博客。 链接地址: https://blog.sciencenet.cn/blog-395566-420073.html
上一篇:
怎样才能避免不谨慎和想当然 下一篇:
perl的代码阅读-1