error comparing sequences - string interpreted as number

Question

I'm trying to do something similar with my previous question.
My purpose is to join all sequences that are equal. But this time instead of letters, I have numbers.

alignment file can be found here - phylip file

the problem is when I try to do this:

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

I get this error:

ValueError: Sequence 1 length 49, expected length 1001000000100000100000001000000000000000

I don't understand why because this is the second file I'm creating and the first one worked perfectly..

Below is the code used to build the alignment file:

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')

So it shouldn't interpret 1001000000100000100000001000000000000000 as a number since its a string..

Any ideas?

Thank you!

Answer 1

Your PHYLIP file is broken. The header says 161 sequences but there are 166. Fixing that the current version of Biopython seems to load your file fine. Maybe use len(info_plex) when creating the header line.

P.S. It would have been a good idea to include the version of Biopython in your question.

Answer 2

The code of Kevin Jacobs in your former question employs Biopython that uses sequences of type Seq that

« are essentially strings of letters like AGTACACTGGT, which seems very natural since this is the most common way that sequences are seen in biological file formats. »

« There are two important differences between Seq objects and standard Python strings. (...)

First of all, they have different methods. (...)

Secondly, the Seq object has an important attribute, alphabet, which is an object describing what the individual characters making up the sequence string “mean”, and how they should be interpreted. For example, is AGTACACTGGT a DNA sequence, or just a protein sequence that happens to be rich in Alanines, Glycines, Cysteines and Threonines?

The alphabet object is perhaps the important thing that makes the Seq
object more than just a string. The currently available alphabets for Biopython are defined in the Bio.Alphabet module. »

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

The reason of your problem is simply that SeqIO.parse() can't create Seq objects from a file containing characters for which there is no alphabet attribute able to manage them.

.

So, you must use another method. Not try to plate an inadapted method on a different problem.

Here's my way:

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

result

len(records) == 165

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

Number of unique occurences : 156

.

Edit

I've understood MY problem: I let 'fasta' instead of 'phylip' in my code.

'phylip' is a valid value for the attribute alphabet, with it it works fine

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)

produces

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

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

There is an incredible amount of characters ,,,,,,,,,,,,,,,, before the interesting data, I wonder what it is.

.

But my code isn't useless. See:

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'

result

12.4826178327 seconds
0.228640588399 seconds

ratio = 55 !