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

blackjimmy

问题描述：
我这里有人类基因组的染色体序列共24个文件，分别是chr1.fa, chr2.fa, ..., chr22.fa, chrX.fa, chrY.fa， 他们在当前文件夹下，格式是这样的：

1. >chrM
   
2. GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT
   
3. TTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTG
   
4. GAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATT
   
5. CTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACCTACTA
   
6. AAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAAT
   
7. GTCTGCACAGCCGCTTTCCACACAGACATCATAACAAAAAATTTCCACCA
   
8. AACCCCCCCCTCCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGC
   
9. CAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAAT
   
10. TTTATCTTTAGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACA
    
11. ......

复制代码

我需要作的就是按照三个参数来提取序列，即chr*, startnum, endnum。 就是取子序列拉。

另外，我还有一个按照chr, startnum, endnum, name组成的一个fragment片段文件。例如fragment_file.txt, 他的内容如下：

1. Chr        Mstart        Mends        NAME
   
2. chrX        1234        12345        a
   
3. chr1        127104        127121        b
   
4. chr2        12961        13025        c
   
5. chr3        138015        138125        d
   
6. chr4        128035        128071        e
   
7. chr5        128063        128120        f
   
8. chr6        128026        128169        g
   
9. chr7        128036        128218        h

复制代码

当然这个fragment_file.txt的文件有成千上万行。

我需要按照这个文件，并根据chr*.fa文件将序列提取出来，并输入到另外一个文件中，格式如下：

1. >a        chrX:12340-12345
   
2. ttagag
   
3. >b        chr1:127114-127121
   
4. CAGCGTGG
   
5. >c        chr2:13016-13025
   
6. ttgtaaatat
   
7. >d        chr3:138115-138125
   
8. ATAACTTGAGG
   
9. >e        chr4:128065-128071
   
10. cggatac
    
11. >f        chr5:128116-128120
    
12. tTCCC
    

复制代码

我的程序如下，希望能改的高效些

1. import time                  
   
2. print time.ctime()
   
3. 
   
4. def get_oneseq(chr, start, end):
   
5.         """
   
6.         this function used to get one sequence from fasta files
   
7.         """
   
8.        
   
9.         f = open(chr+'.fa', 'r')
   
10.         sequence = ''.join(f.readlines()[1:]).replace('\n', '')
    
11.         start_end_seq = sequence[start-1:end]
    
12.         return start_end_seq
    
13.         f.close()
    
14. 
    
15. def get_seqs(fragment):
    
16.         """
    
17.         retrieve fragment from a file
    
18.         """
    
19.        
    
20.         f = open(fragment, 'r')
    
21.         f2 = open(fragment+'.fa', 'w')
    
22.         lists = f.readlines()[1:]       
    
23.         for line in lists:
    
24.                 temp = line.replace('\n', '').split('\t')
    
25.                 chr = temp[0]
    
26.                 start = temp[1]
    
27.                 end = temp[2]
    
28.                 name = temp[3]
    
29.                 seq = get_oneseq(chr, int(start), int(end))
    
30.                 print >>f2, '>'+name+'\t'+chr+':'+str(start)+'-'+str(end)+'\n'+seq
    
31.         f.close()
    
32.         f2.close()
    
33. 
    
34. if __name__ == '__main__':
    
35.         get_seqs('fragment_file.txt')         
    
36. print time.ctime()
    

复制代码

谢谢。

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不是很大可以考虑全部打开读入，不要反复打开关闭。

blackjimmy

chr*.fa文件共有近3G，就是压缩了也有近1G。故没法上传阿，整个读取的话，也不是很现实阿。
luffy.deng 的改进方法我实验了一下，基本上没有太大效果，十分感谢先。
而我看有的程序，提取10000多条就1分多钟，而我的程序差不多1秒才1条。
用gzip打开的话，就更慢了。

另外：人家的程序是c写的，并且预先将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来写。

flw

BIO Perl

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 编辑 ]

reiase

Cython