博文
小鼠的生物信息实验操作总结--调侃效率
||||
最近我被老板的绝招打败了，心服口服。
#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]);
https://m.sciencenet.cn/blog-395566-420073.html
上一篇：怎样才能避免不谨慎和想当然
下一篇：perl的代码阅读-1
有感之余，决定总结一把自己在用的，一级刚刚学到手的各类效率秘籍，算是个总结笔记吧。如果对其他人有帮助，那本鼠就更开心了。同时欢迎各位来本博客做做，磕磕瓜子，交流交流生物信息的各种心得，事实证明：生物信息跟拿基因枪一样，都是手艺活嘛.
我只臭屁的希望若有哥们觉得有点用，然后转载时，请注明出处：即科学网本鼠的博问链接，哦哈哈哈.
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;
}
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;
}
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);
my $result=sqrt( ($_[0]-$_[3])**2 + $_[1]-$_[4])**2 +$_[3]-$_[5])**2);
｝
行了，就到这儿吧，天色已晚，本鼠该回窝了..各位看官，来磕瓜子啊..
https://m.sciencenet.cn/blog-395566-420073.html
上一篇：怎样才能避免不谨慎和想当然
下一篇：perl的代码阅读-1
扫一扫，分享此博文
全部作者的其他最新博文
全部精选博文导读
相关博文
- • Minerals线下恳谈会：履践致远、与时偕行——对话中国科学院广州地球化学研究所期刊合作学者
- • 聚英才 建高地 | 北京理工大学“特立青年学者”全球招聘开启
- • 700年后日本或濒临灭绝？日本学者推算预测：届时或仅剩1名15岁以下孩子
- • [转载]【同位素视角】非英语母语学者如何区分’e.g.’, ‘i.e.’, ‘namely’与‘such as’等混淆难题
- • 美国佐治亚大学等机构学者：刈割策略对Bulldog 805紫花苜蓿+Tifton 85狗牙根混播草地产量及品质的影响
- • 美国堪萨斯州立大学、密苏里大学等机构学者研究成果：土壤水分管理策略和品种多样性对紫花苜蓿产量、营养品质和农场盈利能力的影