[PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？

书上说， hash大了，用each 来遍历，  即

while ( ($key,$val) = each(%hash) ) {

 

是不是真的效率高很多？ 要多大的hash 用each 才有效果？

 

 

 

2009/9/25 空格 <ribozyme2004@gmail.com>

楼上的几位牛人说献丑是谦虚，我这里拿得出来的就真是“丑”了。
请各位高人指正。
按照目前的考虑，出现很少的次数（比如10次）和没有出现对后面的其他计算是一样的。所以这个程序仍然是计数的。另外，它可以按用户要求的长度计算寡核
苷酸串的出现频率。我正在用它算12个核苷酸的串的情况。
程序里用的还是一个散列，如前面各位说的，它对内存的消耗非常大。目前看，计算12核苷酸（12bp）长度的串使用了2G左右的内存。计算10bp串使
用大约800M内存。我还没有算13bp的情况。truncatei 说的调整为整型的那个我还没有加上去。

agentzh兄弟的那个appear.cpp程序里，有一点必须考虑的是寡核苷酸串的互补串的情况。这个说起来比较麻烦。
简单说，就是找到一个15bp串的时候，要同时把这个序列掉个，再按照AT互补GC互补的情况处理。生成的序列同样是应该考虑的一种情况。我最早说4G
的数据，其实一般的基因组都是2G多一点，就是因为考虑互补链才会加倍。



#!/usr/bin/perl -w
use warnings;
use strict;

my ($dir_a,$i,$file,@order,$m,%hash,@array,$temp,$line,@seq,$seq,
$str,@str,$index,$char,$n,$key,$val,$num);
print "Plz input the length of Oligonucleotide you want to count:";
my $len = <>;
chomp $len;

print "Input the file of genome sequence[G], or use default files
[D]...";
$i=<>;
chmop $i;
if ( $i eq 'D' ) {
#在缺省情况下，对某个目录中所有形如_temp_12.fa的文件进行计算。这是我把基因组文件切成50M左右的临时文件。
   $dir_a = "/home/ribozyme/data/genomes/panda/";
   opendir DIR, $dir_a or die;
   $i = 0;
   while ( $file = readdir(DIR) ) {
       if ( $file =~ m{_temp_\d+} ) {  $order[$i] = $dir_a.$file; $i++;  }
   }
   close DIR;
} elsif ( $i eq 'G' ) {
   $order[0]=$i;
} else {
   die "Plz input right choice.\n";

}

my %code_of = (
   A => 0,
   G => 1,
   C => 2,
   T => 3
   );

$m = 0;
$i = 0;

foreach $file (@order) {
   open FA,qq($file); # 打开基因组序列文件
   $temp = '';
   while ( $line = <FA> ) {

#如果是以大于号开头，则表示是fasta格式的序列标题行，用N代替之。一个文件中可能有多个fasta格式的序列
       if ( $line =~ m{^>} ) {
           $temp .= 'N';

# 如果不是以大于号开头，则表示是真正的序列，其中包含大写的ATGCN五种字母。
       } elsif ( $line !~ m{^>} ) {

# N表示测序不清楚，只能确定那里有一个核苷酸的位置。在基因组文件中很多的N ^_^
           chomp $line;
           $temp .= $line;
       }
       else {}
   }
   close FA;
   print length($temp)." raw data obtained for $file.\n";

   $i=0;
# 把所有连续的非N的长度在所要测的寡核苷酸串以上的序列都存入序列数组。
   while ($temp=~m{([^N]{$len,})}g) {  $seq[$i]=$1; $i++;  }
   print "$#seq sequences for the file...\n";

   foreach $seq (@seq) {
       $i = 0;
       do {
           $str[0] = uc( substr($seq,$i,$len) ); #从分段的基因组序列中截取一个短串
           $str[1] = lc( reverse($str[0]) );
           $str[1] =~ s{a}{T}g;
           $str[1] =~ s{t}{A}g;
           $str[1] =~ s{g}{C}g;
           $str[1] =~ s{c}{G}g;    #得到该短串的互补串。
           foreach $str (@str) {
# 计算一个15bp核苷酸串对应的下标，每个字符占用两个二进制位

               my $index = 0;
               for $char (split q[], $str) { # 按字符分开

                   $index = $index << 2; # 左移两位
                   $index += $code_of{$char}; # 将当前位加入末两位
               }
               if ( exists( $hash{$index} ) ) {
                   $hash{$index}++;  # 如hash中有键值则计数
               } else {
                   $hash{$index} = 1;  # 如hash中没有则添加键值
                   $m++;
               }
           }
           $i++;
       } until ( length($seq) < $i + $len );
   }
   print qq($i 序列已处理，散列中有\t $m 个记录。\n);
}

#遍历散列并转换数字下标为核苷酸串
$n=0;
while ( ($key,$val) = each(%hash) ) {
   $str='';
   do {
       $num = $key & 2;
       if ( $code_of{$char} == $char ) {
           $str .= $char;
       } else {
           print "impossible!!\n";
       }
       $key = $key >> 2;
   } until ( $key == 0 );
   print qq($str\t$val\n);
   $n++;
}
print "There are totally $m kinds of oligonucs.\n";
print qq(Normal quit!\n);

--
           Yours Sincerely
                   Zeng Hong

--~--~---------~--~----~------------~-------~--~----~
您收到此信息是由于您订阅了 Google 论坛"PerlChina Mongers 讨论组"论坛。
 要在此论坛发帖，请发电子邮件到 perlchina@googlegroups.com
 要退订此论坛，请发邮件至 perlchina+unsubscribe@googlegroups.com
 更多选项，请通过 http://groups.google.com/group/perlchina?hl=zh-CN 访问该论坛

-~----------~----~----~----~------~----~------~--~---