注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Q超越兔子的蜗牛O--逸云沙鸥Linux

飘飘何所似,天地一沙鸥;落霞与孤鹜齐飞,秋水共长天一色~~

 
 
 

日志

 
 

伟大的尝试,使用了27G内存的Perl程序  

2009-12-24 17:35:34|  分类: Bioinformatics |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

题目描述:(有点小复杂)lxc想用matlab来处理一批数据,chr1.depth 中存有位点坐标(有跳跃) 和相应的匹配覆盖度,chr1.methy_rate中含有位点坐标和该点的甲基化概率,先要统计chr1.methy_rate中概率大于0.7的为正样本,为零的为负样本;再将正负样本中的位点坐标到chr1.depth中匹配,查看其相应位点处 上下50bp位置都连续的点的覆盖度,并且按照:

1)正样本匹配结果 +1 1:2 2:2 3:2 4:2 5:2 6:2 7:2 8:2 9:2 10:2 ...

2)负样本匹配结果 -1 1:2 2:2 3:2 4:2 5:2 6:2 7:2 8:2 9:2 10:2 ...

的模式,将每个样本点采集100个数据为一行,输出到结果文件中。麻烦的是 以上两个文件数据量很大,都在700Mb 左右。为了用大内存加速计算,我直接采用了hash,一次性将相关文件读入到内存。写好脚本文件,成功投递到集群机,花了约半个小时,完成计算,正负样本采集数据量如下:

(老l用C++标准库写的程序,计算了16个小时产出数据目前都才十几Mb)

124M Dec 24 15:46 sample.po

1.3G Dec 24 15:46 sample.ng

[叹]:真正体会到了用集群做生物信息的数据分析的必要性和好处!

附上关键脚本,供纪念: 

#!/usr/bin/perl -w
 2 #compare_qunero.pl
 3 
 4 if (4 != @ARGV){
 5     die " usage:compare_qunero.pl depth_file rate_file sample_p sample_n!";
 6 }
 7 
 8 my @files = @ARGV;
 9 my $depth_file =  shift @files;
10 my $rate_file =  shift @files;
11 my $pos_file = shift @files;
12 my $neg_file = shift @files;
13 print "$depth_file $rate_file $pos_file $neg_file \n";
14 
15 open DEPTH,"<","$depth_file";
16 open RATE,"<","$rate_file";
17 open POS,">","$pos_file";
18 open NEG,">","$neg_file";
19 
20 my %index;
21 my %depth;
22 $loc=0;
23 while(<DEPTH>){
24     ($dep,$val)=split;
25     $index{$dep}=$loc++;
26     $depth{$dep}=$val;
27 }
28 
29 my %rate_p;
30 my %rate_n;
31 my $value;
32 my $key;
33 my $i;
34 my $start;
35 while(<RATE>){
36     ($key,$value)=split;
37     if($value==0){
38         $rate_n{$key}=0;
39     }elsif($value>0.7){
40         $rate_p{$key}=$value;
41     }#else do nothing
42 }
43 
44 foreach $key(sort keys %rate_p){
45     $start=$key-51;
46     if(100 == $index{$key+50}-$index{$key-50}){
47         print POS "+1 ";        
48         for($i=1 ; $i<=100 ; $i++){
49             print POS "$i:$depth{$start+$i";
50         }
51         print POS "\n";
52     }
53 }
54 
55 foreach $key(sort keys %rate_n){
56     $start=$key-51;
57     if(100 == $index{$key+50}-$index{$key-50}){
58         print NEG "-1 ";
59         for($i=1 ; $i<=100 ; $i++){
60             print NEG "$i:$depth{$start+$i";
61         }
62         print NEG "\n";
63     }
64 }
65 
  评论这张
 
阅读(747)| 评论(5)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017