楼上的几位牛人说献丑是谦虚，我这里拿得出来的就真是“丑”了。
请各位高人指正。
按照目前的考虑，出现很少的次数（比如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);
--~--~---------~--~----~------------~-------~--~----~
您收到此信息是由于您订阅了 Google 论坛“PerlChina Mongers 讨论组”论坛。
要在此论坛发帖，请发电子邮件到 [email protected]
要退订此论坛，请发邮件至 [email protected]
更多选项，请通过 http://groups.google.com/group/perlchina?hl=zh-CN 访问该论坛
-~----------~----~----~----~------~----~------~--~---
