用python提取人类基因组序列片段（请帮修改的高效些）

efhilt

C和Py的效率相比差10倍。

此外，如果你的startnum，endnum表示第几个字符的话，应该用用file object的seek，以及os.read来读文件，这样效率应该高很多。

最后，这是一个很简单的文件操作问题。不要像有些人那样，听说过一个bio perl的名字，就以为可以包打天下了，连问题都没搞清楚就推荐给别人。

luffy.deng

原则上我同意efhilt 的观点，楼主程序最主要的问题应该不在于用字符io还是用二进制io,而在于把大量的时间浪费在把一个上百兆文件内容生成到一个列表里然后切片，而且这种操作还要进行上万次。可以参照我在9楼给你的示例，直接 seek()试一下。

liukaiyi

我在 lz 上加点哈！

1. 使用 一个 散列表 ，把 “fragment片段文件” 按顺序排列（此hashmap 当成 监听类表）。
      如果 fragment 文件也过大，那就 linux  sort 下 ，分批放到 内存的 散列中

2.建立多个 文件游标  类表（可认为简单的线程池，每个就文件打开一次） ，按 散列移动到下一个 key startnum 上开始记录。
     如果有多fragment 重叠位置，这也是优化的地方，目前没想到什么好办法。
   
这样就应该块点了

方法当然也很多，比如
   1.文件 每100个 char 都用 替换为    $(chr1.fa,25600){TTAA......}
       然后合并 这 22 个文件 ，然后逻辑也水到渠成了，当然没此有改动就要合并，做成增量合并，一游标  块状（100 优化值可调） 读取从头到脚来此 ok ，速度上去了(复杂度也上了)
   2. fragment 按文件 类型  linux sort 开个文件 来个游标 全部此文件问题解决关闭，下个文件接着
   

关于什么语言，我就 不掺和了
1.python 感觉对面向对象，图片处理，web服务，文件同步等还是不错的。
2.perl 对文本处理不错的（正则在这语言上登峰造极了，如果你的正则不好也大可不使用了）

[ 本帖最后由 liukaiyi 于 2009-2-16 10:23 编辑 ]

luffy.deng

原帖由 liukaiyi 于 2009-2-16 09:56 发表
我在 lz 上加点哈！
1. 使用 一个 散列表 ，把 “fragment片段文件” 按顺序排列（此hashmap 当成 监听类表）。
      如果 fragment 文件也过大，那就 linux  sort 下 ，分批放到 内存的 散列中

2.建立多个 文件游标  类表（可认为简单的线程池，每个就文件打开一次） ，按 散列移动到下一个 key startnum 上开始记录。
     如果有多fragment 重叠位置，这也是优化的地方，目前没想到什么好办法。
     

第1条你认真看看程序再说，fragment片段文件一次性全部读取逐行处理的，散列有啥实际意义？
第2条倒是可行，可以省去反复打开关闭ch*.fa的时间。

caesarok

start = int(temp[1])
IndexError: list index out of range

reiase

生物学里用perl很多，Gentoo的portage tree里很多生物学软件都是perl写的，我们就暂且相信用perl不用C是个比较明智的选择
python的优势是容易借助C提速，有很多软件，比如Cython能够把部分python代码翻译成C语言，从而提高速度，LZ可以从这方面考虑

lenin

darkduke

没想到居然还有同行了。
一般fragment这种文件我都先用shell 里面的sort先按照chromosome排序，这样每次切换的只是一个fasta文件。
具体的读fasta的我用的函数是这样的：

def getChromSeq(chrom_name):
    seq = '';
    for line in fileinput.input([chrom_name]):
        if '>' not in line:
            line = string.strip(line);
            line = string.upper(line);
            seq = seq+line;
    return seq;

因为我的repeat一般都是另外用annotation的，所以这里都换成大写了就是。
然后在main里，如果当前fasta文件不是前一个文件，则另外读入新的文件。
其实我一般都是拆成24个任务，然后每个任务提交到集群节点上去运行的，呵呵。

我python是自学的，也请大家指教了。

zhouyubio

Two alternative ways:

1. use UCSC genome browser tool, twoBitToFa.
2. use pygr library. You just need to build one seqdb such as BlastDB('hg18.fa'), then you can get the sequence for any genomic interval quickly.

That's also what I am using. Good luck.