一个基于位数组的序列覆盖率分析器 (Was Re: [PerlChina] Re: 在我遇到的这种情况下散列和数组哪个快？)
agentzh
>
> 今天晚些时候，我会提供一个 C++ 实现 ;)
>
刚刚写了一个 ANSI C++ 实现，呵呵，也不过 180 行代码，在 Linux 和 Win32 上进行了简单的测试。代码我作为 appears 项目放在了 GitHub 上面：
http://github.com/agentzh/appears
如果没有 git 或者不熟悉 git，可以直接从下面的页面下载 v0.01 版的压缩包：
http://github.com/agentzh/appears/downloads
解压后的编译和用法见 README 和 Makefile.
简单说来，在 Win32 上使用
nmake -f NMakefile
来编译。然后用命令
appears.exe input.txt
来分析输入文件 input.txt
在 Linux 上则直接
make
程序用法是
./appears input.txt
为了测试方便，上面的项目构造过程还会生成 appears1 和 appears2 这两个程序，分别对应序列长度为 1 和 2 的情形（默认的 appears 对应 15）。
仅进行了简单的测试，欢迎同学们提供更彻底的基于 Perl 的测试集，哈哈！
在 Linux 利用 pmap 检查 appears 程序的内存占用情况：
total 134036K
果然只有 100 多 MB 的 RAM 占用，哈哈！
稍后，我得空了，再利用 Inline 模块封装为 Perl 可调用的形式，呵呵。
C++ 写得不好，欢迎大家指正 :)
Cheers,
-agentzh
Beckheng Lam
空格
On 9月23日, 下午6时04分, agentzh <agen...@gmail.com> wrote:
> 2009/9/23 agentzh <agen...@gmail.com>:
>
>
>
> > 今天晚些时候，我会提供一个 C++ 实现 ;)
>
> 刚刚写了一个 ANSI C++ 实现，呵呵，也不过 *180 行代码*，在 Linux 和 Win32 上进行了简单的测试。代码我作为 appears
> 项目放在了 GitHub 上面：
>
> http://github.com/agentzh/appears
>
> 如果没有 git 或者不熟悉 git，可以直接从下面的页面下载 v0.01 版的压缩包：
>
> http://github.com/agentzh/appears/downloads
>
> 解压后的编译和用法见 README 和 Makefile.
>
> 简单说来，在 Win32 上使用
>
> nmake -f NMakefile
>
> 来编译。然后用命令
>
> appears.exe input.txt
>
> 来分析输入文件 input.txt
>
> 在 Linux 上则直接
>
> make
>
> 程序用法是
>
> ./appears input.txt
>
> 为了测试方便，上面的项目构造过程还会生成 appears1 和 appears2 这两个程序，分别对应序列长度为 1 和 2 的情形（默认的
> appears 对应 15）。
>
> 仅进行了简单的测试，欢迎同学们提供更彻底的基于 Perl 的测试集，哈哈！
>
> 在 Linux 利用 pmap 检查 appears 程序的内存占用情况：
>
> total 134036K
>
> *果然只有 100 多 MB 的 RAM 占用*，哈哈！
空格
op - 22:26:25 up 13:09, 1 user, load average: 1.20, 1.03, 0.86
Tasks: 144 total, 2 running, 142 sleeping, 0 stopped, 0 zombie
Cpu(s): 47.2% us, 3.2% sy, 0.0% ni, 49.7% id, 0.0% wa, 0.0% hi,
0.0% si
Mem: 2019000k total, 705476k used, 1313524k free, 13848k
buffers
Swap: 20972816k total, 429208k used, 20543608k free, 270020k
cached
PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
10886 ribozyme 20 0 57716 7092 4464 R 93 0.4 7:06.10 konsole
23717 ribozyme 20 0 130m 128m 676 S 7 6.5 0:40.05 appears
3509 ribozyme 20 0 410m 86m 18m S 1 4.4 59:36.64 firefox
2275 root 20 0 399m 31m 5948 S 1 1.6 17:12.06 X
可以看到，appears占用了2G的内存的6.5%，正好是大约130M。和几位大大估算的非常一致。因为程序是把结果打印到屏幕的，所以
konsle有占用很多的cpu资源。
On 9月23日, 下午6时04分, agentzh <agen...@gmail.com> wrote:
> 2009/9/23 agentzh <agen...@gmail.com>:
>
>
>
> > 今天晚些时候，我会提供一个 C++ 实现 ;)
>
> 刚刚写了一个 ANSI C++ 实现，呵呵，也不过 *180 行代码*，在 Linux 和 Win32 上进行了简单的测试。代码我作为 appears
> 项目放在了 GitHub 上面：
>
> http://github.com/agentzh/appears
>
> 如果没有 git 或者不熟悉 git，可以直接从下面的页面下载 v0.01 版的压缩包：
>
> http://github.com/agentzh/appears/downloads
>
> 解压后的编译和用法见 README 和 Makefile.
>
> 简单说来，在 Win32 上使用
>
> nmake -f NMakefile
>
> 来编译。然后用命令
>
> appears.exe input.txt
>
> 来分析输入文件 input.txt
>
> 在 Linux 上则直接
>
> make
>
> 程序用法是
>
> ./appears input.txt
>
> 为了测试方便，上面的项目构造过程还会生成 appears1 和 appears2 这两个程序，分别对应序列长度为 1 和 2 的情形（默认的
> appears 对应 15）。
>
> 仅进行了简单的测试，欢迎同学们提供更彻底的基于 Perl 的测试集，哈哈！
>
> 在 Linux 利用 pmap 检查 appears 程序的内存占用情况：
>
> total 134036K
>
> *果然只有 100 多 MB 的 RAM 占用*，哈哈！
agentzh
从一个基因组中取了50M大小的一段序列测试了一下，比较成功。top显示的结果如下：
Yay! 我刚刚又对 appears 进行了一些改进，现发布 v0.02：
http://github.com/agentzh/appears/downloads
主要改动如下：
1. 代码命名调整：使用 token 来表示核苷，而非有歧议的 char 这个术语
2. 在代码中进一步去除“神秘的数”
3. 按 gcc 要求启用 large file 支持（原先的 v0.01 在 Linux 上读取 4 GB 的文件会报错： "Value too large for defined data type"）
4. 添加了 t/gensample.pl 和 t/gensample.c 这两个程序，用来生成指定大小的随机序列文件（默认为 4.8 G）。
gensample.pl 和 gensample.c 的功能相同，但性能有明显差别：Perl 实现大约 1 秒生成 1 MB 的序列，而后者，即 C 实现， 1 秒大约生成 20 MB 的序列。
我使用 gensample.c 程序生成了一个 4 GB 大小的随机序列文件，其输出如下：
$ time t/gensample > samples/big.txt
real 2m25.392s
user 2m11.432s
sys 0m6.148s
$ ls -lh samples/big.txt
-rw-r--r-- 1 agentz agentz 4.0G 2009-09-24 10:20 samples/big.txt
$ time ./appears samples/big.txt > /tmp/out.txt
INFO: Reading file samples/big.txt...
INFO: Using a bit array of length 1073741824 (sequence length: 15).
INFO: Searching missing combinations of sequences...
real 7m36.456s
user 6m59.882s
sys 0m9.977s
在我的 Core2Duo T9600, SATA 7600rpm disk 的本本上，一共只需要 7 分半钟的时间，呵呵。
可以看到，appears占用了2G的内存的6.5%，正好是大约130M。和几位大大估算的非常一致。因为程序是把结果打印到屏幕的，所以
konsle有占用很多的cpu资源。
结果最好重定向到磁盘文件，就像上面所演示的那样，呵呵。
-agentzh
空格
ftp://ftp.ensembl.org/pub/release-55/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.55.dna.toplevel.fa.gz
下载速度还是蛮快的。^_^
> 在我的 Core2Duo T9600, SATA 7600rpm disk 的本本上，*一共只需要 7 分半钟*的时间，呵呵。
空格
panic: realloc at ChrY-0.2-bit.pl line 23, <FA> line 42954286.
好象是什么"堆溢出的恐慌"。。。怎么会这样呢？
//bow
On 9月24日, 上午11时02分, agentzh <agen...@gmail.com> wrote:
> 2009/9/23 空格 <ribozyme2...@gmail.com>
> 在我的 Core2Duo T9600, SATA 7600rpm disk 的本本上，*一共只需要 7 分半钟*的时间，呵呵。
空格
的标题。下面才是序列的内容。另外一个建议就是真实的测序文件中会有N出现。这表示测序结果不确定的情况。在处理序列时可以把N分隔的序列逐段保存成一
个字符数组再逐个计数。
On 9月24日, 上午11时02分, agentzh <agen...@gmail.com> wrote:
> 2009/9/23 空格 <ribozyme2...@gmail.com>
>
> 在我的 Core2Duo T9600, SATA 7600rpm disk 的本本上，*一共只需要 7 分半钟*的时间，呵呵。
agentzh
另外，提几个小建议。一个是自动删除序列中以大于号开头的行。因为基因组序列保存一般都是fasta格式。就是第一行用一个大于号开头，后面是这个序列
的标题。下面才是序列的内容。另外一个建议就是真实的测序文件中会有N出现。这表示测序结果不确定的情况。在处理序列时可以把N分隔的序列逐段保存成一
个字符数组再逐个计数。
另一个值得一提的地方是，如果系统中使用的是带 pthread 支持的 glibc，则在 appears.cpp 中使用 getc_unlocked 代替 getc 可以观察到 70% 以上的总性能提升。我试了一下，处于 4 GB 文件的总用时从 7 分半下降到 2 分。对 t/gensample.c 施用类似的 s/putc/putc_unlocked/ 亦会看到近 70% 的提速。回头得空了，我再加一下 autoconf/automake 支持，以便能自动检测系统中是否有 *_unlocked 函数，呵呵。。。
Cheers,
-agentzh
agentzh
一个问题，请问您是怎样打开一个大文件的。我打开一个2.4G大小的基因组文件会报错：
panic: realloc at ChrY-0.2-bit.pl line 23, <FA> line 42954286.
好象是什么"堆溢出的恐慌"。。。怎么会这样呢？
//bow
你使用的 perl 5 解释器在编译时启用了 large file 支持了么？呵呵。
Cheers,
-agentzh
空格
另外，还有一个小建议，就是可不可以加一个参数。把这个寡核苷酸串的长度设定成按照用户需求设定的一个值。
On 9月24日, 下午4时51分, agentzh <agen...@gmail.com> wrote:
> 2009/9/24 空格 <ribozyme2...@gmail.com>
agentzh
另一个值得一提的地方是，如果系统中使用的是带 pthread 支持的 glibc，则在 appears.cpp 中使用 getc_unlocked 代替 getc 可以观察到 70% 以上的总性能提升。我试了一下，处于 4 GB 文件的总用时从 7 分半下降到 2 分。
Cheers,
-agentzh
agentzh
好象确实是这个大文件支持没有打开。正在研究perl解释器中，呵呵。
另外，还有一个小建议，就是可不可以加一个参数。把这个寡核苷酸串的长度设定成按照用户需求设定的一个值。
if (tokens_seen < APPEARS_SEQSIZE) {
tokens_seen++;
}
须替换为
tokens_seen++;
即无条件递增。
由于总共不到 200 行 C++ 代码，而且我在可读性上面下了些功夫，应该不难修改和调试的，呵呵。
Cheers,
-agentzh