书上说， hash大了，用each 来遍历， 即
while ( ($key,$val) = each(%hash) ) {
是不是真的效率高很多？ 要多大的hash 用each 才有效果？
2009/9/25 空格 <[email protected]>
> 楼上的几位牛人说献丑是谦虚，我这里拿得出来的就真是“丑”了。
> 请各位高人指正。
> 按照目前的考虑，出现很少的次数（比如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 讨论组”论坛。
要在此论坛发帖，请发电子邮件到 [email protected]
要退订此论坛，请发邮件至 [email protected]
更多选项，请通过 http://groups.google.com/group/perlchina?hl=zh-CN 访问该论坛
-~----------~----~----~----~------~----~------~--~---
