2

我试图创建一个读取fasta文件“seqs.fa”的python程序,并让程序按名称对序列进行排序。

Fasta 文件如下所示:

>seqA - human
GCTGACGTGGTGAAGTCAC
>seqC - gorilla
GATGACAA
GATGAAGTCAG
>seqB - chimp
GATGACATGGTGAAGTAAC

我的程序如下所示:

import sys

inFile = open(sys.argv[1], 'r')
a = inFile.readlines()
a.sort()
seq = ''.join(a[0:])
seq = seq.replace('\n', "\n")
print seq

预期结果:

>seqA - human
GCTGACGTGGTGAAGTCAC
>seqB - chimp
GATGACATGGTGAAGTAAC
>seqC - gorilla
GATGACAAGATGAAGTCAG

我的结果:

>seqA - human
>seqB - chimp
>seqC - gorilla
GATGACAA
GATGAAGTCAG
GATGACATGGTGAAGTAAC
GCTGACGTGGTGAAGTCAC

最后四行是大猩猩、黑猩猩和人类序列,大猩猩序列分为前两行。

谁能给我一些关于如何排序或解决问题的方法的提示?

4

3 回答 3

5

不要自己实现 FASTA 阅读器!像大多数情况一样,有一些聪明的人已经为你做了这件事。改为使用例如BioPython。像这样:

from Bio import SeqIO
handle = open("seqs.fa", "rU")
l = SeqIO.parse(handle, "fasta")
sortedList = [f for f in sorted(l, key=lambda x : x.id)]
for s in sortedList:
   print s.description
   print str(s.seq)
于 2012-05-11T02:36:09.413 回答
3

您的代码存在一些问题。主要是在readlines()您的描述和序列返回的列表中都是单独的行,因此当您对列表进行排序时,它们是相互分离的。此外,所有描述都在序列之前,因为它们'>'在开头。

二、a[0:]同理a

第三,seq.replace('\n', "\n")什么都不做。单引号和双引号的意思是一样的。您用自己替换换行符。

读取 fasta 文件对于 Python 来说并不是一项非常复杂的任务,但我仍然希望能够原谅我提供使用我正在研究的包 - pyteomics

这是我要使用的代码:

In [1]: from pyteomics import fasta

In [2]: with fasta.read('/tmp/seqs.fa') as f:
   ...:     fasta.write(sorted(f))
   ...:     
>seqA - human
GCTGACGTGGTGAAGTCAC

>seqB - chimp
GATGACATGGTGAAGTAAC

>seqC - gorilla
GATGACAAGATGAAGTCAG

要将其保存到新文件,请将其名称fasta.write作为参数:

fasta.write(sorted(f), 'newfile.fa')

通常,pyteomics.fasta是针对蛋白质序列,而不是 DNA,但它确实可以完成这项工作。也许您可以使用它返回元组中的描述和序列的事实。

于 2012-05-10T22:34:48.050 回答
0
file = open("seqs.fa")    
a = file.readlines()
i = 0
ar = []
while True:
    l1=file.readline()
    l2=file.readline()
    if not (l1 and l2):
        break;
    l = l1.strip('\n') + '////////' + l2
    ar.append(l)
ar = ar.sort()
for l in ar:
    l1 = l.split('////////')[0]+'\n'
    print l1
    l2 = l.split('////////')[1]
    print l2
于 2012-05-10T22:38:20.467 回答