I read a fasta-formatted genome sequence and try to get its reverse complement:
genomeSeq = FastaIO.FastaIterator(genomeHandle, IUPACUnambiguousDNA).next()
genomeSeq.seq.reverse_complement()
but it doesn't work and I can't understand why:
File "/Users/charodeika/Dropbox/genesGelfand/scripts/genome/src/matrixCount/sigma.py", line 127, in <module>
print genomeSeq.seq[:10].reverse_complement()
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/Bio/Seq.py", line 804, in reverse_complement
return self.complement()[::-1]
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/Bio/Seq.py", line 752, in complement
base = Alphabet._get_base_alphabet(self.alphabet)
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/Bio/Alphabet/__init__.py", line 213, in _get_base_alphabet
"Invalid alphabet found, %s" % repr(a)
AssertionError: Invalid alphabet found, <class 'Bio.Alphabet.IUPAC.IUPACUnambiguousDNA'>
And, of course, I would like to find out how to make it work!
What I think the problem is, that genome sequences contain N's for nucleotides of which they are not sure what nucleotide it is. IUPACUnambiguousDNA only accepts ACTG. So you would want to change your script to ignore the N's (or any other non-dna letter there is in there.
No, certainly not. It doesn't work for slices, that do not contain any other letters as A, T, G,C
Well, when I added
from Bio.Alphabet import IUPAC
and changed
genomeSeq.seq.reverse_complement()
to
genomeSeq = FastaIO.FastaIterator(genomeHandle, IUPAC.unambiguous_dna).next()
it started working nicely. However, I would still be happy to learn why.
Log in to answer this question.
