在我遇到的这种情况下散列和数组哪个快？
空格
性。我想统计一下，这个大字符串中是否包含了所有的长度为15的可能的字串。如果没有包含全部，那么有哪些字串的出现次数为零。
为此，我想需要建立一个很大的表，然后从那个超大的字符串中逐个取出长度为15的字串，然后在表中统计其出现次数。这样可以得到结果。
我的问题是，这样大的表格，用散列写好还是用二维数组写比较好？或者有什么别的方式实现更可行一些。
谢谢各位～
Michael Zeng
--
Yours Sincerely
Zeng Hong
Qiang (James) Li
> 有一个长度为4.8G的字符串，其中只有四种字母ATGC。按照排列组合数，这四个字母组成的长度为15字符串总共有1｀073｀741｀824种可能
> 性。我想统计一下，这个大字符串中是否包含了所有的长度为15的可能的字串。如果没有包含全部，那么有哪些字串的出现次数为零。
> 为此，我想需要建立一个很大的表，然后从那个超大的字符串中逐个取出长度为15的字串，然后在表中统计其出现次数。这样可以得到结果。
不知道 bioperl 里是否有现成的工具，偶不是搞生物的。
不过你匹配的时候把满足要求的用 $hash{...}++ 直接统计就可以了，不用把不符合要求的也放到 hash 里。
> 我的问题是，这样大的表格，用散列写好还是用二维数组写比较好？或者有什么别的方式实现更可行一些。
hash 即可。
Qiang
Jester
这种hash直接处理当然是最简单的，不过我觉得可能需要考虑一下内存的问题。
我遇到过处理比较大的序列时，一不小心就out of memory了。:(
不知道哪位有比较节约内存的方法？
Jester，jes...@perlchina.org
2009-09-22
----- Original Message -----
From: Qiang (James) Li
To: perlchina
Sent: 2009-09-22, 00:03:19
Subject: [PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？
>2009/9/21 空格 :
黄浩
agentzh
如果一定要用 table 去统计各种 15 字节的字符串个数，可以考虑使用 TokyoCabinet (TC)，并预留大约 10 GB 的磁盘空间，并将 bucket number 置为 15! 的 2 到 3 倍。扫描字符串需要做成流式的，比如几 KB 几 KB 或者几 MB 几 MB 地从文件读取。照你所说，扫描过程中，把 15 字节字符串在 table （即 tokyocabinet）中 set value，最后看 TC 数据库中一共有多少个 key-value 对，如果不足 15!，则回答了你的第一个问题“这个大字符串中是否包含了所有的长度为15的可能的字串”。而至于找出没有出现过的 15 长字串的所有实例，似乎也比较花时间，需要遍历 15! 种 15 长度的串，分别测试它是否已出现在了 TC 数据库里。
TC 占用的空间大约可以估计为 15! * (15 + 4) * 2，大约是 40 GB 吧，需要预留好磁盘空间，这个倒是可以接受的。
但如果假设 TC 平均读写速度为 0.1ms，则遍历一遍 15! 个组合，貌似需要 3404 年！呃。。。好像还没有生物有那么长寿吧？呵呵。当然了，如果找几百台机器实现几千个 TC 操作并行起来，可以缩减到 1 年，还是有些长。。。那再多找些盘吧。。。
上面是蛮力法。更聪明的做法可能需要一些组合数学方面的规律了。。。我数学不太好，还需要多想想。。。
当然了，如果有许多台机器，总 RAM 达到几十 GB 的量级的话，就不必用 TC，也不必走磁盘了。以 RAM 10G/s 的读写速度算，在运算时间上还是可以接受的。。。但最好不要用 Perl 哈，直接上 C/C++ 吧！呵呵！Lua
Cheers,
-agentzh
msmouse
不知道我的分析对不对
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
agentzh
没有15！那么大的量级吧 应该是4的15次幂，也就是2的30次幂，大约1G种组合 ，用一个int32数组计数即可，散列的话考虑到很多种组合是不存在的 有2G内存也是可能一次跑成的
嗯嗯嗯，忘了是有重复元素的列表了。。。哈哈，多谢指正。总排列数确实是 4^15 :) 每一个位子都只有 4 种可能性，便是 4*4*4*...*4 这 15 个 4 相乘。为避免存储各个 key，将 key 作 hash 到一个 int32 整数，便只有 4 个字节。如果不保存计数值，用位数组来代替哈希表的话，应该更省空间，1 GB RAM 作数据本身的存储，再留个几 MB RAM 作其他计算用 :) 保存次数并用数组的话，用 int32 作为数组元素则需要 4 GB 了吧？不知道这一次我有没有算错，呵呵 :P
Cheers,
-agentzh
agentzh
嗯嗯嗯，忘了是有重复元素的列表了。。。哈哈，多谢指正。总排列数确实是 4^15 :) 每一个位子都只有 4 种可能性，便是 4*4*4*...*4 这 15 个 4 相乘。为避免存储各个 key，将 key 作 hash 到一个 int32 整数，便只有 4 个字节。
- 事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。
- 然后把当前的15长度的串顺序地两个比特两个比特地编码，最终得到的值便是我们要的 hash 后的数值。
- 接着用它去下标访问我们的 1G 大小的位数组，这个寻址当是极快的。置比特也当是极快的。
- 最后，我们遍历这个 1G 位数组，找出比特为 0 的下标，再还原为 ATGC 形式的碱基序列，便是从未出现过的 15 长度的串了，呵呵。
不知这个是否靠谱？哈哈？
Cheers,
-agentzh
agentzh
- 接着用它去下标访问我们的 1G 大小的位数组，这个寻址当是极快的。置比特也当是极快的。
msmouse
如果要perl并且需要计数的话，考虑用一个带有cache的嵌入式数据库即可。。前面也提到了吧？
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
agentzh
agentzh
嗯，１
该死的 Gmail 快捷键。。。唉。。。1G 个比特便是 125 M 个字节，确实 125MB 就够了 :)
Cheers,
-agentzh
空格
今天试着写了一下，确实对内存的消耗很大。我正在想办法。。。
On 9月22日, 上午10时20分, 黄浩 <huanghao19...@gmail.com> wrote:
> 可以每次处理一部分， 分而治之。
> 为了提高运行效率，Perl的实现采取是典型的用空间换时间的做法。所以对内存消耗特别大，而如果用hash的话，消耗就会更大。所以只能每次读一部分到内存，一段段的处理。
>
> 2009/9/22 Jester <jes...@perlchina.org>
>
>
>
> > bioperl好像没有这样的用法......
> Email：huang...@nrchpc.ac.cn
空格
15长度的字串是在序列上依次取的。因而一个长度为N bp 的DNA 序列会有(100－15)*2 种字串，当然可能有重复的。。。
On 9月21日, 下午10时48分, Michael Zeng <galaxy2...@gmail.com> wrote:
> 嗬嗬， 原来是搞生物的，
>
> 用hash好了
>
> $hash{ $key} ++ ;
>
> 不知道你那个15长度的字符串怎么取的
>
> 2009/9/21 空格 <ribozyme2...@gmail.com>
空格
有两个问题我不太明白：
“为避免存储各个 key，将 key 作 hash 到一个 int32 整数，便只有 4个字节。”
这个是不是说，我把四个字母换成数字就可以了？
我尝试着用正则匹配把字符串换成01的数字串，就像这样：$tempa=~s{A}{00}g;
但是似乎并没有很节省内存啊。
您说的位数组我不了解。好像在perl里不能实现，只能在C／C++里实现是么？
On 9月22日, 下午12时52分, agentzh <agen...@gmail.com> wrote:
> 2009/9/22 msmouse <msmo...@gmail.com>
空格
On 9月22日, 下午1时08分, msmouse <msmo...@gmail.com> wrote:
> 同意 C++的用一个vector<bool>或者bitset就搞定了，不需要自己做位计算，只需要125M内存 哈哈
>
> 如果要perl并且需要计数的话，考虑用一个带有cache的嵌入式数据库即可。。前面也提到了吧？
>
> ----------------------------------
> msmo...@ir.hit.edu.cn
> msmo...@gmail.com
>
> 2009/9/22 agentzh <agen...@gmail.com>
>
> > 2009/9/22 agentzh <agen...@gmail.com>
>
> >> 嗯嗯嗯，忘了是有重复元素的列表了。。。哈哈，多谢指正。总排列数确实是 4^15 :) 每一个位子都只有 4 种可能性，便是 4*4*4*...*4
> >> 这 15 个 4 相乘。为避免存储各个 key，将 key 作 hash 到一个 int32 整数，便只有 4 个字节。
>
> > 当然，我在写上一封邮件的时候其实并不知道这样的 hash 函数是否好写。。。呵呵，应该是好写的：
>
> > 1. 事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。
> > 2. 然后把当前的15长度的串顺序地两个比特两个比特地编码，最终得到的值便是我们要的 hash 后的数值。
> > 3. 接着用它去下标访问我们的 1G 大小的位数组，这个寻址当是极快的。置比特也当是极快的。
> > 4. 最后，我们遍历这个 1G 位数组，找出比特为 0 的下标，再还原为 ATGC 形式的碱基序列，便是从未出现过的 15 长度的串了，呵呵。
msmouse
这个想法就是把这个内存放不下的数组放在硬盘上。由于这1G种不同组合很可能不是均匀分布的，也就是说其中某些更经常出现，那么如果存储引擎有cache的话，就可以加快访问。但是实际的效果怎么样没有试就不好说，因为在这个问题里数据访问频度虽然很可能不均匀，但是也是没有什么局部性的，如果是按块的cache（比如OS本身的disk cache）对这种应用的效果可能就不大理想，因为相邻访问之间涉及的地址不是连续的，cache必须不停换页。
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
Michael Zeng
空格
我的电脑是4G内存的，下午测试发现最多只能装下大约3*10^8种15bp的字串。这个比1*10^9种情况可是少很多了。是不是hash表不可能更
节省了。我根据这个数目正在做9bp长度的字串分析。
如果要再长，是不是就只能从硬盘上分段去做了？
On 9月22日, 下午10时09分, Michael Zeng <galaxy2...@gmail.com> wrote:
> 要实在提高效率 只能用 C/C++ 语言比较好吧
> 他们有很多字符串函数
>
> perl 只是写起来方便， 其本身是 c语言写出来的，效率怎么能高过c呢
>
> - 显示引用的文字 -
空格
基因组字串的分布应该还是比较均匀的。因为它很接近随机序列。cache的效果如何很难说。。。
我原来想的就是用简单的文本文件处理。
因为hash表会按键值排序，所以相似的字串在文本上比较也会比较容易些吧。
On 9月22日, 下午6时44分, msmouse <msmo...@gmail.com> wrote:
> 比如sqlite。。
>
> 这个想法就是把这个内存放不下的数组放在硬盘上。由于这1G种不同组合很可能不是均匀分布的，也就是说其中某些更经常出现，那么如果存储引擎有cache的话，-就可以加快访问。但是实际的效果怎么样没有试就不好说，因为在这个问题里数据访问频度虽然很可能不均匀，但是也是没有什么局部性的，如果是按块的cache（比-如OS本身的disk
> cache）对这种应用的效果可能就不大理想，因为相邻访问之间涉及的地址不是连续的，cache必须不停换页。
> ----------------------------------
> msmo...@ir.hit.edu.cn
> msmo...@gmail.com
>
> 2009/9/22 空格 <ribozyme2...@gmail.com>
>
>
>
> > 抱歉你说的这个"带cache的嵌入式数据库"我不懂。具体是什么，在cpan上有包么？
>
> > On 9月22日, 下午1时08分, msmouse <msmo...@gmail.com> wrote:
> > > 同意 C++的用一个vector<bool>或者bitset就搞定了，不需要自己做位计算，只需要125M内存 哈哈
>
> > > 如果要perl并且需要计数的话，考虑用一个带有cache的嵌入式数据库即可。。前面也提到了吧？
>
> > > ----------------------------------
> > > msmo...@ir.hit.edu.cn
> > > msmo...@gmail.com
>
> > > 2009/9/22 agentzh <agen...@gmail.com>
>
> > > > 2009/9/22 agentzh <agen...@gmail.com>
>
> > > >> 嗯嗯嗯，忘了是有重复元素的列表了。。。哈哈，多谢指正。总排列数确实是 4^15 :) 每一个位子都只有 4 种可能性，便是
> > 4*4*4*...*4
> > > >> 这 15 个 4 相乘。为避免存储各个 key，将 key 作 hash 到一个 int32 整数，便只有 4 个字节。
>
> > > > 当然，我在写上一封邮件的时候其实并不知道这样的 hash 函数是否好写。。。呵呵，应该是好写的：
>
> > > > 1. 事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。
> > > > 2. 然后把当前的15长度的串顺序地两个比特两个比特地编码，最终得到的值便是我们要的 hash 后的数值。
> > > > 3. 接着用它去下标访问我们的 1G 大小的位数组，这个寻址当是极快的。置比特也当是极快的。
> > > > 4. 最后，我们遍历这个 1G 位数组，找出比特为 0 的下标，再还原为 ATGC 形式的碱基序列，便是从未出现过的 15
> > 长度的串了，呵呵。
>
> > > > 这里我们就不关心出现的 15 长度序列的具体次数了。只记录出现过或者未出现过。这样 2 GB RAM 的机器是很够的了 ;)
>
> > > > 不知这个是否靠谱？哈哈？
>
> > > > Cheers,
> > > > -agentzh- 隐藏被引用文字 -
>
> - 显示引用的文字 -
Jester
我感觉前面有人讨论的建议的使用位运算应该是会比较节约的做法，但具体的实现我也不会。
我有个不太成熟的想法，不太知道数组的存储是否会比hash节约？如果是，到底能节约多少？（请高手指教）
我的想法是把ATCG用0123代替，然后看作是一个4进制的数，这样可以很容易转换成一个数值，也就是数组的下标，如此对应。
当然，如果数组的开销还是很大，就不行了......
我不知道你在哪个实验室，国内应该有不少做类似分析的人，郝先生今年刚发在NAR webserverissue上的CVTree本质上就是分析这个。
说到底用C还是最高效的，不过用perl尝试一下也挺好。我也是丢了C很多年了，执着的perler ^_^
请高手再指点吧，我也想提高一下coding技巧。
Jester，jes...@perlchina.org
2009-09-23
From: 空格
To: PerlChina Mongers 讨论组
Sent: 2009-09-22, 23:33:47
Subject: [PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？
>原来是这个意思。。。。
>基因组字串的分布应该还是比较均匀的。因为它很接近随机序列。cache的效果如何很难说。。。
>我原来想的就是用简单的文本文件处理。
>因为hash表会按键值排序，所以相似的字串在文本上比较也会比较容易些吧。
>
>> 比如sqlite。。
>>
>> 这个想法就是把这个内存放不下的数组放在硬盘上。由于这1G种不同组合很可能不是均匀分布的，也就是说其中某些更经常出现，那么如果存储引擎有cache的话，-就可以加快访问。但是实际的效果怎么样没有试就不好说，因为在这个问题里数据访问频度虽然很可能不均匀，但是也是没有什么局部性的，如果是按块的cache（比-如OS本身的disk
>> cache）对这种应用的效果可能就不大理想，因为相邻访问之间涉及的地址不是连续的，cache必须不停换页。
>> ----------------------------------
>> msmo...@ir.hit.edu.cn
>> msmo...@gmail.com
>>
>> 2009/9/22 空格
>>
>>
>>
>>
>> > > 同意 C++的用一个vector或者bitset就搞定了，不需要自己做位计算，只需要125M内存 哈哈
>> > > 如果要perl并且需要计数的话，考虑用一个带有cache的嵌入式数据库即可。。前面也提到了吧？
>>
>> > > ----------------------------------
>> > > msmo...@ir.hit.edu.cn
>> > > msmo...@gmail.com
>>
>> > > 2009/9/22 agentzh
>>
msmouse
”事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。“
即用两个位表示一个字符，15个字符即是30位。30个二进制位即是15个4进制位
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
agentzh
> 呵呵，有四五年没有碰C了,现在看到没有$的变量名都不习惯。还是在perl里倒腾吧。
今天晚些时候，我会提供一个 C++ 实现 ;)
Stay tuned!
Cheers,
-agentzh
PIG
在 09-9-23，agentzh<age...@gmail.com> 写道：
刘鑫
"基因组字串的分布应该还是比较均匀的。因为它很接近随机序列"----这样的说法很不严谨！我相信这不是你真实的想法，但是这样的表述会误导很多人。任何一个物种的基因组序列都不会是ATCG的随机排布！
我感觉前面有人讨论的建议的使用位运算应该是会比较节约的做法，但具体的实现我也不会。
我有个不太成熟的想法，不太知道数组的存储是否会比hash节约？如果是，到底能节约多少？（请高手指教）
我的想法是把ATCG用0123代替，然后看作是一个4进制的数，这样可以很容易转换成一个数值，也就是数组的下标，如此对应。
当然，如果数组的开销还是很大，就不行了......
我不知道你在哪个实验室，国内应该有不少做类似分析的人，郝先生今年刚发在NAR webserverissue上的CVTree本质上就是分析这个。
说到底用C还是最高效的，不过用perl尝试一下也挺好。我也是丢了C很多年了，执着的perler ^_^
请高手再指点吧，我也想提高一下coding技巧。
位运算我刚基于Postgres做了一个数据库实现，基本上放几个亿还是没有问题的：）
--
光见贼吃肉，没见贼挨打。
……
劉鑫
March.Liu
Jester
From: 刘鑫
To: perlchina
Sent: 2009-09-23, 14:17:22
Subject: [PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？
>
>位运算我刚基于Postgres做了一个数据库实现，基本上放几个亿还是没有问题的：）
>
空格
作为整体尺度上承受选择压力的作用而偏离随机分布的量是很小的。
我们用互信息建立了一个衡量标准。当考虑双核甘酸串时，对于无限长的随机序列是零值，而基因组序列的这个值基本在0.02-0.05左右。对于更长的核
甘酸串，这个值更低。这个方法最早是罗辽复等提出的。按照罗的计算结果，对于编码序列，当寡核甘酸串长度在3以上时，其互信息值基本都在随机波动阈值以
下。我们正在做的就是考虑全基因组的这个互信息值的情况。
On 9月23日, 上午11时02分, "Jester"<jes...@perlchina.org> wrote:
> "基因组字串的分布应该还是比较均匀的。因为它很接近随机序列"----这样的说法很不严谨！我相信这不是你真实的想法，但是这样的表述会误导很多人。任何一个物种的基因组序列都不会是ATCG的随机排布！
>
空格
就是说，
当取到一个字串赋给$key之后，
$key=~s{A}{0}g;
$key=~s{G}{1}g;
$key=~s{T}{2}g;
$key=~s{C}{3}g;
再用 exists(%hash{$key}); 检查是否有特定键值。是这样的吧？
On 9月23日, 上午11时17分, msmouse <msmo...@gmail.com> wrote:
> 就是这个意思。前面agentzh说：
> ”事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。“
> 即用两个位表示一个字符，15个字符即是30位。30个二进制位即是15个4进制位
> ----------------------------------
空格
算了有1.6G的样子。内存消耗维持在60%占用的水平上。我的内存是4G的，估计这个hash大小有2.4G左右的样子吧。12bp字串的种类有
16'777'216种。不过我用的是十位数字的0,1,10,11来表示的ATGC。昨晚犯的错，呵呵。
郝先生的文章我拜读过，其实我的思路也有些是受他的启发。所不同的是，CVTree是做系统发育树的。基本不会用到9bp以上的寡核甘酸，而且目前
CVTree也主要是构建原核生物树，只有到了很少的几种真核生物基因组。
On 9月23日, 上午11时02分, "Jester"<jes...@perlchina.org> wrote:
msmouse
my @count_of;
# 字符和数字的转换表
my %code_of = (
A=>0,
G=>1,
C=>2,
T=>3
);
# 计算一个15字符的字串对应的下标，每个字符占用两个二进制位
my $index = 0;
for $char (split q[], $str) { # 按字符分开
$index *= 4; # 左移两位
$index += $key{$char}; # 将当前位加入末两位
}
# 计数：
$count_of[$index] += 1;
打印出来需要通过一个循环取末两位、右移两位的方法取得下标对应的实际字符串
不知道这样说清楚没有 ：）
但是如昨天的讨论所说 这个数组本身就至少占4G内存（假定每项是一个4字节的int），因此在perl中完全在内存中进行是不大可行的。
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
truncatei
2位表示一个字符
00 A
01 G
10 T
11 C
这样下来一个15个字符的序列一共30位，一个int足矣，必用字符串省很多内存
取各个位
如：
$char = $key & 0b000011;
$char = $key & 0b001100;
$char = $key & 0b110000
...
取出每个位以后再做判断
$char == 0;
$char == 1;
$char == 2;
$char == 4;
空格
算了有1.6G的样子。内存消耗维持在60%占用的水平上。我的内存是4G的，估计这个hash大小有2.4G左右的样子吧。12bp字串的种类有
16'777'216种。不过我用的是十位数字的0,1,10,11来表示的ATGC。昨晚犯的错，呵呵。
郝先生的文章我拜读过，其实我的思路也有些是受他的启发。所不同的是，CVTree是做系统发育树的。基本不会用到9bp以上的寡核甘酸，而且目前
CVTree也主要是构建原核生物树，只有到了很少的几种真核生物基因组。
On 9月23日, 上午11时02分, "Jester"<jes...@perlchina.org> wrote:
空格
On 9月23日, 下午5时40分, msmouse <msmo...@gmail.com> wrote:
> # 用一个1G项的数据计数
> my @count_of;
> # 字符和数字的转换表
> my %code_of = (
> A=>0,
> G=>1,
> C=>2,
> T=>3
> );
>
> # 计算一个15字符的字串对应的下标，每个字符占用两个二进制位
> my $index = 0;
> for $char (split q[], $str) { # 按字符分开
> $index *= 4; # 左移两位
> $index += $key{$char}; # 将当前位加入末两位
>
> }
>
> # 计数：
> $count_of[$index] += 1;
>
> 打印出来需要通过一个循环取末两位、右移两位的方法取得下标对应的实际字符串
>
> 不知道这样说清楚没有 ：）
> 但是如昨天的讨论所说 这个数组本身就至少占4G内存（假定每项是一个4字节的int），因此在perl中完全在内存中进行是不大可行的。
>
> ----------------------------------
> msmo...@ir.hit.edu.cn
> msmo...@gmail.com
>
> 2009/9/23 空格 <ribozyme2...@gmail.com>
空格
在perl里如何实现你说的这个二进制的bit运算？抱歉我不知道用什么语句能实现这个“2位表示一个字符”。
另外这个：
$char = $key & 0b000011;
On 9月23日, 下午5时49分, truncatei <trunca...@gmail.com> wrote:
> bit运算应该是这样：
> 2位表示一个字符
> 00 A
> 01 G
> 10 T
> 11 C
>
> 这样下来一个15个字符的序列一共30位，一个int足矣，必用字符串省很多内存
>
> 取各个位
> 如：
> $char = $key & 0b000011;
> $char = $key & 0b001100;
> $char = $key & 0b110000
> ...
> 取出每个位以后再做判断
> $char == 0;
> $char == 1;
> $char == 2;
> $char == 4;
>
> 2009/9/23 空格 <ribozyme2...@gmail.com>
Beckheng Lam
以上纯属猜测，如有雷同，实属撞。。。
--
Perl乐事 -- http://www.perlersh.org
我的博客 -- http://www.perlersh.org/blog.html
msmouse
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
msmouse
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
msmouse
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
msmouse
他演示的是从一个计算出来的index里面取两个二进制位的方法
看一下agentzh的c++程序吧，看懂了你就明白了。。
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
Michael Zeng
msmouse
跟c对比觉得不方便的时候就是它不能按字符索引操作而已，或者说跟c的区别在于perl字符串不是字符数组
就现在的问题来看，操作一个128M的字符串，当做一个1G下标的位数组来用，其实也没有太不方便了。。当然要慢很多，毕竟不是数组
总之 这个事情perl不好做是因为资源不够用。。而不是处理字符串不强。。否则一个hash不就搞定了
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
truncatei
一个ASCII字符一般是一byte是占用8个bit (二进制)
agent大侠的算法应该是用位实现的(我今天没时间down下来研究)
这里有关于perl位运算的文档
http://perldoc.perl.org/perlop.html#Bitwise-And
空格
我之前一直想的是数字，而不是数位。
尤其这句：
$index *= 4; # 左移两位
希望我这次是真明白点了。。。 //blush...
On 9月23日, 下午9时50分, msmouse <msmo...@gmail.com> wrote:
> 不是32位是30位，，写顺手了。。
> ----------------------------------
> msmo...@ir.hit.edu.cn
> msmo...@gmail.com
>
> 2009/9/23 msmouse <msmo...@gmail.com>
>
> > 不是15位十进制数，是32位二进制 只是我没把1234写成 0b00 0b01 0b10 0b11，没把$index*4写成$index<<2而已。。
> > ----------------------------------
truncatei
agentzh
真的猜对了
一个ASCII字符一般是一byte是占用8个bit (二进制)
agent大侠的算法应该是用位实现的(我今天没时间down下来研究)
http://github.com/agentzh/appears/blob/master/appears.cpp
真的要感谢 DNA/RNA 世界是“四进制”的，于是可以很容易地和计算机的二进制系统对齐（毕竟是 2 的幂），这给核苷 GATC 到二进制比特的正反向编码映射，带来了极大的便利。如果遗传物质是“三进制”的话，appears.cpp 真的会复杂许多，呵呵。
Cheers,
-agentzh
刘鑫
还是perl 社区的牛人多，那个C++实现太赞了。
我前两周的时候，也是正好遇到一个做生物的朋友提这个问题，所以写了一个基于Postgres的数据库bimap映射实现，还很简陋，正在整理资料。
数据结构很简单：
create table bitmap(bigint idx, segment bit(32));
其它计算方面的东西相信就都很好写了，这里我献丑放一个记录覆盖的：
create or replace function segment_true(s bigint, e bigint) returns void as $$
declare
sidx bigint;
eidx bigint;
d_data bit(32) := 4294967295::bit(32); -- pow(2, 32) -1
begin
if currval('bitmap_idx_seq') < (e/32)+1 then
for i in currval('bitmap_idx_seq')..e/32+1 loop
insert into bitmap(segment) select 0::bit(32);
end loop;
end if;
sidx := s/32;
eidx := e/32;
if sidx = eidx then
update bitmap set segment = (d_data << cast((32-e%32) as int))&(d_data >> cast((s%32) as int)) where idx = sidx;
else
update bitmap set segment = segment | (d_data >> cast((s%32) as int)) where idx = sidx;
update bitmap set segment = d_data where sidx < idx and idx < eidx;
update bitmap set segment = segment | (d_data << cast((e%32) AS INT)) where idx = eidx;
end if;
end;
$$ language plpgsql;
说明一下，PG的整数没有无符号类型，所以操作起来比较讨厌，我干脆直接用它的bit类型，这个类型最多可以支持2^31位，所有基础的位运算它都支持。通过idx可以将分片的bit记录连成一个大的bitmap。具体bit应该多少位，idx应该用int还是bigint，我想应该视数据量和业务来调整，不必太死板，能达到好用，快速就好。
数据库实现的好处是可以在数据量和性能以及并发性上有个比较折中的实现，而且因为都是现成的功能，组合起来比较快，质量也可靠。
空格兄的问题很有意思，我想写写试试。
空格
请各位高人指正。
按照目前的考虑，出现很少的次数（比如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);
Michael Zeng
空格
# 寡核苷酸串长度
len of oligonuc string = 12
#来源文件目录
directory of genome sequences = /home/ribozyme/data/genomes/galga/
#是否遍历该目录
traversal the directory = N
遍历目录时所用的正则信息
key words of the files = _temp_\d+\.fsa
#不遍历目录时，给定的文件名
name of the files = Gallus_gallus.WASHUC2.55.dna_rm.chromosome.28.fa
#结果文件名称
result file name = /home/ribozyme/bin/ChrY/hash_12bp.dat
# Finish the ctrl file.
Perl程序文件中也做了修改，如下：
#!/usr/bin/perl -w
use warnings;
use strict;
my (@ctl,$len,$dir_a,$i,$file,@order,$m,%hash,@array,$temp,$line,@seq,
$seq,$str,@str,$index,$char,$n,$key,$val,$num);
open FA,q(/home/ribozyme/bin/ChrY/OligoNuc.ctl) or die "Can't find the
control file\n";
$i = 0;
while ( $line = <FA> ) {
next if $line=~m{^#};
$line =~ m{\s\=\s(.+)};
$ctl[$i]=$1;
$i++;
}
close FA;
$len = $ctl[0];
$dir_a = $ctl[1];
opendir DIR, $dir_a or die "Directory can't be opened!\n";
if ( $ctl[2] =~ m{^Y$}i ) {
$i = 0;
while ( $file = readdir(DIR) ) {
if ( $file =~ m{$ctl[3]} ) { $order[$i] = $dir_a.$file; $i++; }
}
} elsif ( $ctl[2] !~ m{^Y$}i ) {
$order[0] = $dir_a.$ctl[4];
}
close DIR;
my %code_of = (
A => 0,
G => 1,
C => 2,
T => 3
);
%hash = ();
$m = 0;
$i = 0;
foreach $file (@order) {
open FA,qq($file); # 打开基因组序列文件
$temp = '';
while ( $line = <FA> ) {
if ( $line =~ m{^>} ) { #如果是以大于号开头，则表示是fasta格式的序列标题行，用N代替之。一个文件中可能有多个
fasta格式的序列
$temp .= 'N';
} elsif ( $line !~ m{^>} ) { # 如果不是以大于号开头，则表示是真正的序列，其中包含大写的ATGCN五种字母。
chomp $line; # N表示测序不清楚，只能确定那里有一个核苷酸的位置。
$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 "$i 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);
}
#遍历散列并转换数字下标为核苷酸串
open FB,qq(>$ctl[5]);
$n=0;
while ( ($index,$val) = each(%hash) ) {
$str='';
$i=0;
do {
$num = $index & 2;
foreach $key ( keys %code_of ) {
if ( $code_of{$key} == $num ) {
$str .= $key;
}
}
$index = $index >> 2;
$i++;
} until ( $i == $len );
print FB qq($str\t$val\n);
$n++;
}
close FB;
print "There are totally $m kinds of oligonucs.\n";
print qq(Normal quit!\n);
On 9月25日, 下午1时56分, 空格 <ribozyme2...@gmail.com> wrote:
空格
On 9月25日, 下午11时04分, Michael Zeng <galaxy2...@gmail.com> wrote:
> 书上说， hash大了，用each 来遍历， 即
> while ( ($key,$val) = each(%hash) ) {
>
> 是不是真的效率高很多？ 要多大的hash 用each 才有效果？
>
>
Michael Zeng
agentzh
简单说，就是找到一个15bp串的时候，要同时把这个序列掉个，再按照AT互补GC互补的情况处理。生成的序列同样是应该考虑的一种情况。
这个可以在向 BitArray 注册序列时，把其“同构异形体也一并注册了 :) 在输出比特为 0 的对应序列的过程中，亦可将当前输出序列的同构异形体也在 BitArray 中注册为比特 1，这样后面就不会重复输出了。所以应当是加几行 C++ 代码的交情，哈哈！不知道我说的对不对？
-agentzh
agentzh
程序里用的还是一个散列，如前面各位说的，它对内存的消耗非常大。目前看，计算12核苷酸（12bp）长度的串使用了2G左右的内存。计算10bp串使
用大约800M内存。我还没有算13bp的情况。truncatei 说的调整为整型的那个我还没有加上去。
Cheers,
-agentzh
TBY
需要把4.8g的字串作两次转换。分别将奇数和偶数元素转换成字符串，然后约定各种8字符字符串去找吧。
On 9月23日, 上午11时17分, msmouse <msmo...@gmail.com> wrote:
> 就是这个意思。前面agentzh说：
> "事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。"
> 即用两个位表示一个字符，15个字符即是30位。30个二进制位即是15个4进制位
> ----------------------------------
> msmo...@ir.hit.edu.cn
> msmo...@gmail.com
>
> 2009/9/23 Jester <jes...@perlchina.org>
>
>
>
>
>
> > "基因组字串的分布应该还是比较均匀的。因为它很接近随机序列"----这样的说法很不严谨！我相信这不是你真实的想法，但是这样的表述会误导很多人。任何一个-物种的基因组序列都不会是ATCG的随机排布！
>
> > 我感觉前面有人讨论的建议的使用位运算应该是会比较节约的做法，但具体的实现我也不会。
> > 我有个不太成熟的想法，不太知道数组的存储是否会比hash节约？如果是，到底能节约多少？（请高手指教）
> > 我的想法是把ATCG用0123代替，然后看作是一个4进制的数，这样可以很容易转换成一个数值，也就是数组的下标，如此对应。
> > 当然，如果数组的开销还是很大，就不行了......
>
> > 我不知道你在哪个实验室，国内应该有不少做类似分析的人，郝先生今年刚发在NAR webserverissue上的CVTree本质上就是分析这个。
> > 说到底用C还是最高效的，不过用perl尝试一下也挺好。我也是丢了C很多年了，执着的perler ^_^
> > 请高手再指点吧，我也想提高一下coding技巧。
>
> > Jester，jes...@perlchina.org
> > 2009-09-23
> > ----- Original Message -----
> > From: 空格
> > To: PerlChina Mongers 讨论组
> > Sent: 2009-09-22, 23:33:47
> > Subject: [PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？
>
> > >原来是这个意思。。。。
> > >基因组字串的分布应该还是比较均匀的。因为它很接近随机序列。cache的效果如何很难说。。。
> > >我原来想的就是用简单的文本文件处理。
> > >因为hash表会按键值排序，所以相似的字串在文本上比较也会比较容易些吧。
>
> > >On 9月22日, 下午6时44分, msmouse wrote:
> > >> 比如sqlite。。
>
> > 这个想法就是把这个内存放不下的数组放在硬盘上。由于这1G种不同组合很可能不是均匀分布的，也就是说其中某些更经常出现，那么如果存储引擎有cache的话，--就可以加快访问。但是实际的效果怎么样没有试就不好说，因为在这个问题里数据访问频度虽然很可能不均匀，但是也是没有什么局部性的，如果是按块的cache（-比-如OS本身的disk
> > >> - 显示引用的文字 -- 隐藏被引用文字 -
>
> - 显示引用的文字 -
TBY
有4.8g的碱基组合。
On 9月27日, 下午3时56分, TBY <tangbo...@hotmail.com> wrote:
> 其实我觉得可以直接用16进制来计算，1个16进制数可以代表相邻2个碱基的所有组合可能，15个字符，可以转化为8个0-F的字符串。麻烦的地方在于
> 需要把4.8g的字串作两次转换。分别将奇数和偶数元素转换成字符串，然后约定各种8字符字符串去找吧。
>
> On 9月23日, 上午11时17分, msmouse <msmo...@gmail.com> wrote:
>
>
>
> > 就是这个意思。前面agentzh说：
> > "事先约定 A,T,G,C 分别对应 00, 01, 10, 11，即 2 个比特的数值。"
> > 即用两个位表示一个字符，15个字符即是30位。30个二进制位即是15个4进制位
> > ----------------------------------
> > msmo...@ir.hit.edu.cn
> > msmo...@gmail.com
>
> > 2009/9/23 Jester <jes...@perlchina.org>
>
> > > "基因组字串的分布应该还是比较均匀的。因为它很接近随机序列"----这样的说法很不严谨！我相信这不是你真实的想法，但是这样的表述会误导很多人。任何一个--物种的基因组序列都不会是ATCG的随机排布！
>
> > > 我感觉前面有人讨论的建议的使用位运算应该是会比较节约的做法，但具体的实现我也不会。
> > > 我有个不太成熟的想法，不太知道数组的存储是否会比hash节约？如果是，到底能节约多少？（请高手指教）
> > > 我的想法是把ATCG用0123代替，然后看作是一个4进制的数，这样可以很容易转换成一个数值，也就是数组的下标，如此对应。
> > > 当然，如果数组的开销还是很大，就不行了......
>
> > > 我不知道你在哪个实验室，国内应该有不少做类似分析的人，郝先生今年刚发在NAR webserverissue上的CVTree本质上就是分析这个。
> > > 说到底用C还是最高效的，不过用perl尝试一下也挺好。我也是丢了C很多年了，执着的perler ^_^
> > > 请高手再指点吧，我也想提高一下coding技巧。
>
> > > Jester，jes...@perlchina.org
> > > 2009-09-23
> > > ----- Original Message -----
> > > From: 空格
> > > To: PerlChina Mongers 讨论组
> > > Sent: 2009-09-22, 23:33:47
> > > Subject: [PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？
>
> > > >原来是这个意思。。。。
> > > >基因组字串的分布应该还是比较均匀的。因为它很接近随机序列。cache的效果如何很难说。。。
> > > >我原来想的就是用简单的文本文件处理。
> > > >因为hash表会按键值排序，所以相似的字串在文本上比较也会比较容易些吧。
>
> > > >On 9月22日, 下午6时44分, msmouse wrote:
> > > >> 比如sqlite。。
>
> > > 这个想法就是把这个内存放不下的数组放在硬盘上。由于这1G种不同组合很可能不是均匀分布的，也就是说其中某些更经常出现，那么如果存储引擎有cache的话，---就可以加快访问。但是实际的效果怎么样没有试就不好说，因为在这个问题里数据访问频度虽然很可能不均匀，但是也是没有什么局部性的，如果是按块的cache-（-比-如OS本身的disk
agentzh
其实我觉得可以直接用16进制来计算，1个16进制数可以代表相邻2个碱基的所有组合可能，15个字符，可以转化为8个0-F的字符串。麻烦的地方在于
需要把4.8g的字串作两次转换。分别将奇数和偶数元素转换成字符串，然后约定各种8字符字符串去找吧。
一个 16 进制数可以表示 4 个比特，即 2 个碱基。而 ascii 形式的 16 进制数则需要一个字节，即 8 个比特。从空间效率上讲，多用了一倍的空间，呵呵。而且 ascii 形式的序列在这里确实没有原始二进制的序列处理起来方便 ;)
Cheers,
-agentzh
Huangj
分而治之...
在2009-09-21，"Michael Zeng" <galax...@gmail.com> 写道： 嗬嗬， 原来是搞生物的，
用hash好了
$hash{ $key} ++ ;
不知道你那个15长度的字符串怎么取的
2009/9/21 空格 <ribozy...@gmail.com>
有一个长度为4.8G的字符串，其中只有四种字母ATGC。按照排列组合数，这四个字母组成的长度为15字符串总共有1｀073｀741｀824种可能性。我想统计一下，这个大字符串中是否包含了所有的长度为15的可能的字串。如果没有包含全部，那么有哪些字串的出现次数为零。
为此，我想需要建立一个很大的表，然后从那个超大的字符串中逐个取出长度为15的字串，然后在表中统计其出现次数。这样可以得到结果。我的问题是，这样大的表格，用散列写好还是用二维数组写比较好？或者有什么别的方式实现更可行一些。谢谢各位～-- Yours Sincerely Zeng Hong
"中国制造",讲述中国60年往事
空格
比如我们设定散列：
my %code_of = (
A => 0,
G => 1,
C => 2,
T => 3
);
这样确实可以满足A（00）和T（11）互补以及G（01）和C（10）互补。但是核甘酸序列的书写是有方向的,生物化学中的叫法是从5撇端向3撇端方
向。彼此为互补链的两条序列书写时的方向要反过来。
比如一个四核甘酸：ATAC。对应的互补序列是GTAT而不是TATG。因为他的互补时的方向是这样的：
|-ATAC->
||||
<-TATG-|
这里“|-” 表示DNA链的开始方向，箭头“->”表示链终止方向。竖线表示上下两个核甘酸可以按照碱基互补原则配对。
专业一点的写法是这样：
5`-ATAC-3`
||||
3`-TATG-5`
为了省事，大家约定书写核甘酸序列都按从5撇到3撇的方向写。所以3`-TATG-5`的简写就是GTAT，完整写出来是5`-GTAT-3`
所以，ATAC的位运算表示是00110010，其互补链GTAT是01110011
这个和你说的同构异形还要加一个字符串反向的步骤。
在我写的那个脚本里是按照“先反向再替换”来得到互补链的：
$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; #以上以此替换，得到互补链
On 9月27日, 下午12时35分, agentzh <agen...@gmail.com> wrote:
> 2009/9/25 空格 <ribozyme2...@gmail.com>
空格
On 9月27日, 下午12时46分, agentzh <agen...@gmail.com> wrote:
> 2009/9/25 空格 <ribozyme2...@gmail.com>
>
空格
上去右拐服务台对面第三个还是第四个书架上就是。大概位置是这样：
｜
| PERL！
| 过
--------- 书架
道
书架
-----
|电梯 | ____ 过
| | 服 |
| | 务 |
| | 台 | 道
On 9月26日, 下午9时47分, Michael Zeng <galaxy2...@gmail.com> wrote:
> 那本书有电子版么， 发给我看好么
>
> 2009/9/26 空格 <ribozyme2...@gmail.com>
msmouse
不过15位长的情况下确实还是能够在内存里面做到的，如果是只统计覆盖，不统计具体数目的话，内存还十分够用。。
----------------------------------
msm...@ir.hit.edu.cn
msm...@gmail.com
Michael Zeng
空格
写得还算详细全面。被我当成Perl的字典来用。
On 9月28日, 下午8时27分, Michael Zeng <galaxy2...@gmail.com> wrote:
> 这样的，这本书主要讲啥
>
> 什么技术内幕，高级编程么？