- 论坛徽章:
- 0
|
问题描述:
我这里有人类基因组的染色体序列共24个文件,分别是chr1.fa, chr2.fa, ..., chr22.fa, chrX.fa, chrY.fa, 他们在当前文件夹下,格式是这样的:
- >chrM
- GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT
- TTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTG
- GAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATT
- CTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACCTACTA
- AAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAAT
- GTCTGCACAGCCGCTTTCCACACAGACATCATAACAAAAAATTTCCACCA
- AACCCCCCCCTCCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGC
- CAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAAT
- TTTATCTTTAGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACA
- ......
复制代码
我需要作的就是按照三个参数来提取序列,即chr*, startnum, endnum。 就是取子序列拉。
另外,我还有一个按照chr, startnum, endnum, name组成的一个fragment片段文件。例如fragment_file.txt, 他的内容如下:
- 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文件将序列提取出来,并输入到另外一个文件中,格式如下:
- >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
复制代码
我的程序如下,希望能改的高效些
- 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()
复制代码
谢谢。 |
|