0

我正在尝试对我之前的问题做类似的事情。
我的目的是加入所有相等的序列。但这一次我没有字母，而是数字。

对齐文件可以在这里找到 - phylip 文件

问题是当我尝试这样做时：

records = list(SeqIO.parse(file(filename),'phylip'))

我收到此错误：

ValueError: Sequence 1 length 49, expected length 1001000000100000100000001000000000000000

我不明白为什么，因为这是我正在创建的第二个文件，而第一个文件运行良好..

以下是用于构建对齐文件的代码：

fl.write('\t')
fl.write(str(161))
fl.write('\t')
fl.write(str(size))
fl.write('\n')

for i in info_plex:
    if 'ref' in i[0]:
        i[0] = 'H37Rv'
    fl.write(str(i[0]))
    num = 10 - len(i[0])
    fl.write(' ' * num)
    for x in i[1:]:
        fl.write(str(x))
    fl.write('\n')

所以它不应该将 1001000000100000100000001000000000000000 解释为一个数字，因为它是一个字符串。

有任何想法吗？

谢谢！

4

2 回答 2

2

您的 PHYLIP 文件已损坏。标题说 161 个序列，但有 166 个。修复当前版本的 Biopython 似乎可以很好地加载您的文件。在创建标题行时可能使用 len(info_plex) 。

PS 在您的问题中包含 Biopython 的版本是个好主意。

于 2011-09-06T11:15:12.327 回答
1

Kevin Jacobs 在您之前的问题中的代码使用了BiopythonSeq ，它使用的类型序列

«本质上是像AGTACACTGGT这样的字母串，这看起来很自然，因为这是在生物文件格式中看到序列的最常见方式。»

Seq«对象和标准 Python 字符串之间有两个重要的区别。(...)

首先，他们有不同的方法。(...)

其次，Seq 对象有一个重要的属性，alphabet它是一个对象，它描述了组成序列字符串的各个字符的“含义”，以及应该如何解释它们。例如，AGTACACTGGT 是 DNA 序列，还是恰好富含丙氨酸、甘氨酸、半胱氨酸和苏氨酸的蛋白质序列？

字母表对象可能是使Seq
对象不仅仅是一个字符串的重要因素。Biopython 当前可用的字母在 Bio.Alphabet 模块中定义。»

http://biopython.org/DIST/docs/tutorial/Tutorial.html

您的问题的原因很简单，就是SeqIO.parse()无法Seq从包含没有 alphabet属性能够管理它们的字符的文件中创建对象。

.

所以，你必须使用另一种方法。不要试图在不同的问题上采用不适应的方法。

这是我的方式：

from itertools import groupby
from operator import itemgetter
import re

regx = re.compile('^(\d+)[ \t]+([01]+)',re.MULTILINE)

with open('pastie-2486250.rb') as f:
    records = regx.findall(f.read())

records.sort(key=itemgetter(1))

print 'len(records) == %s\n' % len(records)

n = 0
for seq,equal in groupby(records, itemgetter(1)):
    ids = tuple(x[0] for x in equal)
    if len(ids)>1:
        print '>%s :\n%s' % (','.join(ids), seq)
    else:
        n+=1

print '\nNumber of unique occurences : %s' % n

结果

len(records) == 165

>154995,168481 :
0000000000001000000010000100000001000000000000000
>123031,74772 :
0000000000001111000101100011100000100010000000000
>176816,178586,80016 :
0100000000000010010010000010110011100000000000000
>129575,45329 :
0100000000101101100000101110001000000100000000000

Number of unique occurences : 156

.

编辑

我已经理解了我的问题：我在我的代码中使用了“fasta”而不是“phylip”。

'phylip' 是属性的有效值alphabet，它可以正常工作

records = list(SeqIO.parse(file('pastie-2486250.rb'),'phylip'))

def seq_getter(s):  return str(s.seq)

records.sort(key=seq_getter)

ecr = []

for seq,equal in groupby(records, seq_getter):
    ids = tuple(s.id for s in equal)
    if len(ids)>1:
        ecr.append(  '>%s\n%s' % (','.join(ids),seq) )                   

print '\n'.join(ecr)

生产

,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

>154995,168481
0000000000001000000010000100000001000000000000000
>123031,74772
0000000000001111000101100011100000100010000000000
>176816,178586,80016
0100000000000010010010000010110011100000000000000
>129575,45329
0100000000101101100000101110001000000100000000000

在有趣的数据之前有大量的字符,,,,,,,,,,,,,,,,，我想知道它是什么。

.

但我的代码不是没用的。看：

from time import clock
from itertools import groupby
from operator import itemgetter
import re
from Bio       import SeqIO


def seq_getter(s):  return str(s.seq)
t0 = clock()
with open('pastie-2486250.rb') as f:
    records = list(SeqIO.parse(f,'phylip'))
records.sort(key=seq_getter)
print clock()-t0,'seconds'



t0 = clock()
regx = re.compile('^(\d+)[ \t]+([01]+)',re.MULTILINE)
with open('pastie-2486250.rb') as f:
    records = regx.findall(f.read())
records.sort(key=itemgetter(1))
print clock()-t0,'seconds'

结果

12.4826178327 seconds
0.228640588399 seconds

比率 = 55！

于 2011-09-06T17:54:35.010 回答