伟大的尝试，使用了27G内存的Perl程序
2009-12-24 17:35:34| 分类：
Bioinformatics
| 标签：
|举报
|字号大中小 订阅
题目描述：（有点小复杂）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
［叹］：真正体会到了用集群做生物信息的数据分析的必要性和好处！
附上关键脚本，供纪念：
1 #!/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
评论这张
转发至微博
转发至微博
关闭
玩LOFTER，免费冲印20张照片，人人有奖！
我要抢>
<#--最新日志，群博日志-->
<#--推荐日志-->
<#--引用记录-->
<#--博主推荐-->
<#--随机阅读-->
<#--首页推荐-->
<#--历史上的今天-->
<#--被推荐日志-->
<#--上一篇，下一篇-->
<#-- 热度 -->
<#-- 网易新闻广告 -->
<#--右边模块结构-->
<#--评论模块结构-->
<#--引用模块结构-->
<#--博主发起的投票-->
评论