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

blackjimmy

UID: 30347
帖子: 50
积分: 114
在线时间: 3 小时

1^# blackjimmy 发表于 2009-02-15 12:54

用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()

谢谢。

luffy.deng

UID: 17648
帖子: 1
积分: 2
在线时间: 10 分钟

2^# luffy.deng 发表于 2009-02-15 15:30

QUOTE:

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

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

luffy.deng

UID: 17648
帖子: 1
积分: 2
在线时间: 10 分钟

3^# luffy.deng 发表于 2009-02-15 15:37

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:
      ....

另外能不能把测试文本传上来，光看描述不直观。

luffy.deng

UID: 17648
帖子: 1
积分: 2
在线时间: 10 分钟

4^# luffy.deng 发表于 2009-02-15 15:59

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

blackjimmy

UID: 30347
帖子: 50
积分: 114
在线时间: 3 小时

5^# blackjimmy 发表于 2009-02-15 19:12

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

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

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

luffy.deng

UID: 17648
帖子: 1
积分: 2
在线时间: 10 分钟

6^# luffy.deng 发表于 2009-02-15 19:54

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]快。

luffy.deng

UID: 17648
帖子: 1
积分: 2
在线时间: 10 分钟

7^# luffy.deng 发表于 2009-02-15 19:55

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

flw

UID: 21204
帖子: 144
积分: 331
在线时间: 1 天 4 小时

8^# flw 发表于 2009-02-15 19:57

BIO Perl

luffy.deng

UID: 17648
帖子: 1
积分: 2
在线时间: 10 分钟

9^# luffy.deng 发表于 2009-02-15 20:09

即使是二进制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试了一下还可以，不知道大牛们什么意见？

reiase

UID: 31880
帖子: 179
积分: 411
在线时间: 1 天 19 小时

10^# reiase 发表于 2009-02-15 21:28

Cython