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

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

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

[Copy to clipboard] [ - ]
CODE:
>chrM
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT
TTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTG
GAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATT
CTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACCTACTA
AAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAAT
GTCTGCACAGCCGCTTTCCACACAGACATCATAACAAAAAATTTCCACCA
AACCCCCCCCTCCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGC
CAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAAT
TTTATCTTTAGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACA
......

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

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

[Copy to clipboard] [ - ]
CODE:
Chr        Mstart        Mends        NAME
chrX        1234        12345        a
chr1        127104        127121        b
chr2        12961        13025        c
chr3        138015        138125        d
chr4        128035        128071        e
chr5        128063        128120        f
chr6        128026        128169        g
chr7        128036        128218        h

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

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

[Copy to clipboard] [ - ]
CODE:
>a        chrX:12340-12345
ttagag
>b        chr1:127114-127121
CAGCGTGG
>c        chr2:13016-13025
ttgtaaatat
>d        chr3:138115-138125
ATAACTTGAGG
>e        chr4:128065-128071
cggatac
>f        chr5:128116-128120
tTCCC

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

[Copy to clipboard] [ - ]
CODE:
import time                  
print time.ctime()

def get_oneseq(chr, start, end):
        """
        this function used to get one sequence from fasta files
        """
       
        f = open(chr+'.fa', 'r')
        sequence = ''.join(f.readlines()[1:]).replace('\n', '')
        start_end_seq = sequence[start-1:end]
        return start_end_seq
        f.close()

def get_seqs(fragment):
        """
        retrieve fragment from a file
        """
       
        f = open(fragment, 'r')
        f2 = open(fragment+'.fa', 'w')
        lists = f.readlines()[1:]       
        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
        f.close()
        f2.close()

if __name__ == '__main__':
        get_seqs('fragment_file.txt')         
print time.ctime()

谢谢。


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

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


QUOTE:
原帖由 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:
        ....
                              
另外能不能把测试文本传上来,光看描述不直观。

fragment如果行数很多。chr*.fa文件被反复打开读取并没有修改,如果chr*.fa不是很大可以考虑全部打开读入,不要反复打开关闭。
chr*.fa文件共有近3G,就是压缩了也有近1G。故没法上传阿,整个读取的话,也不是很现实阿。
luffy.deng 的改进方法我实验了一下,基本上没有太大效果,十分感谢先。
而我看有的程序,提取10000多条就1分多钟,而我的程序差不多1秒才1条。
用gzip打开的话,就更慢了。

另外:人家的程序是c写的,并且预先将chr*.fa转成了二进制。

这个问题还是比较经典的生物学问题。


QUOTE:
原帖由 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]快。

或者考虑io部分直接用c来写。
BIO Perl
即使是二进制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试了一下还可以,不知道大牛们什么意见?

Cython