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

luffy.deng

原帖由 blackjimmy 于 2009-2-15 12:54 发表
        start_end_seq = sequence[start-1:end]
        return start_end_seq
        f.close()

先把这儿改一下吧。都return了还close什么？

luffy.deng

原帖由 blackjimmy 于 2009-2-15 12:54 发表
        for line in lists:
                temp = line.replace('\n', '').split('\t')
                chr = temp[0]
                start = temp[1]
                end = temp[2]
                name = temp[3]
                seq = get_oneseq(chr, int(start), int(end))
                print >>f2, '>'+name+'\t'+chr+':'+str(start)+'-'+str(end)+'\n'+seq.

chr start  end name 这些赋值操作和str(start)  str(end)根本没有必要。直接写成get_oneseq(temp[0], int(temp[1]), int(temp[2]))不是一样？如果你想处理一下int操作的异常倒是可以保留start和 end
try:
        start = int(temp[1])
        end = int(temp[2])
except:
        ....
                              
另外能不能把测试文本传上来，光看描述不直观。

[ 本帖最后由 luffy.deng 于 2009-2-15 15:47 编辑 ]

luffy.deng

fragment如果行数很多。chr*.fa文件被反复打开读取并没有修改，如果chr*.fa不是很大可以考虑全部打开读入，不要反复打开关闭。

luffy.deng

原帖由 blackjimmy 于 2009-2-15 19:12 发表
另外：人家的程序是c写的，并且预先将chr*.fa转成了二进制。

二进制 io是要快一点，python中照样可以使用二进制io。加上b标志就可以了。不知道“预先将chr*.fa转成了二进制“是什么意思？难道chr*.fa
行是定长的？那么可以直接seek()应该比 start_end_seq = sequence[start-1:end]快。

[ 本帖最后由 luffy.deng 于 2009-2-16 00:53 编辑 ]

luffy.deng

或者考虑io部分直接用c来写。

luffy.deng

即使是二进制io,下面这种操作执行几万遍也不会很快，
f = open(chr+'.fa', 'r')
sequence = ''.join(f.readlines()[1:]).replace('\n', '')
start_end_seq = sequence[start-1:end]
f.close()
尤其是chr文件有上百兆。因为只是从文件的start读取到end,感觉直接seek()应该可行，比如
    f=open('chr1.fa','rb')
    start=3
    end=6
    f.seek((start-1)*52+7)
    start_end_seq=f.read((end-start)*52)
    f.close()
我在window下用一个200多兆的chr试了一下还可以，不知道大牛们什么意见？

[ 本帖最后由 luffy.deng 于 2009-2-15 21:27 编辑 ]

luffy.deng

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

luffy.deng

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

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

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