--- /dev/null
+# Copyright 2000-2002 Brad Chapman.
+# Copyright 2004-2005 by M de Hoon.
+# Copyright 2007-2009 by Peter Cock.
+# All rights reserved.
+# This code is part of the Biopython distribution and governed by its
+# license. Please see the LICENSE file that should have been included
+# as part of this package.
+"""Provides objects to represent biological sequences with alphabets.
+
+See also U{http://biopython.org/wiki/Seq} and the chapter in our tutorial:
+ - U{http://biopython.org/DIST/docs/tutorial/Tutorial.html}
+ - U{http://biopython.org/DIST/docs/tutorial/Tutorial.pdf}
+"""
+__docformat__ ="epytext en" #Don't just use plain text in epydoc API pages!
+
+import string #for maketrans only
+import array
+import sys
+
+#TODO - Remove this work around once we drop python 2.3 support
+try:
+ set = set
+except NameError:
+ from sets import Set as set
+
+import Alphabet
+from Alphabet import IUPAC
+from Data.IUPACData import ambiguous_dna_complement, ambiguous_rna_complement
+from Bio.Data import CodonTable
+
+def _maketrans(complement_mapping) :
+ """Makes a python string translation table (PRIVATE).
+
+ Arguments:
+ - complement_mapping - a dictionary such as ambiguous_dna_complement
+ and ambiguous_rna_complement from Data.IUPACData.
+
+ Returns a translation table (a string of length 256) for use with the
+ python string's translate method to use in a (reverse) complement.
+
+ Compatible with lower case and upper case sequences.
+
+ For internal use only.
+ """
+ before = ''.join(complement_mapping.keys())
+ after = ''.join(complement_mapping.values())
+ before = before + before.lower()
+ after = after + after.lower()
+ return string.maketrans(before, after)
+
+_dna_complement_table = _maketrans(ambiguous_dna_complement)
+_rna_complement_table = _maketrans(ambiguous_rna_complement)
+
+class Seq(object):
+ """A read-only sequence object (essentially a string with an alphabet).
+
+ Like normal python strings, our basic sequence object is immutable.
+ This prevents you from doing my_seq[5] = "A" for example, but does allow
+ Seq objects to be used as dictionary keys.
+
+ The Seq object provides a number of string like methods (such as count,
+ find, split and strip), which are alphabet aware where appropriate.
+
+ The Seq object also provides some biological methods, such as complement,
+ reverse_complement, transcribe, back_transcribe and translate (which are
+ not applicable to sequences with a protein alphabet).
+ """
+ def __init__(self, data, alphabet = Alphabet.generic_alphabet):
+ """Create a Seq object.
+
+ Arguments:
+ - seq - Sequence, required (string)
+ - alphabet - Optional argument, an Alphabet object from Bio.Alphabet
+
+ You will typically use Bio.SeqIO to read in sequences from files as
+ SeqRecord objects, whose sequence will be exposed as a Seq object via
+ the seq property.
+
+ However, will often want to create your own Seq objects directly:
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> my_seq = Seq("MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF",
+ ... IUPAC.protein)
+ >>> my_seq
+ Seq('MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF', IUPACProtein())
+ >>> print my_seq
+ MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF
+ """
+ # Enforce string storage
+ assert (type(data) == type("") or # must use a string
+ type(data) == type(u"")) # but can be a unicode string
+ self._data = data
+ self.alphabet = alphabet # Seq API requirement
+
+ # A data property is/was a Seq API requirement
+ def _set_data(self, value) :
+ #TODO - In the next release, actually raise an exception?
+ #The Seq object is like a python string, it should be read only!
+ import warnings
+ warnings.warn("Writing to the Seq object's .data propery is deprecated.",
+ DeprecationWarning)
+ self._data = value
+ data = property(fget= lambda self : str(self),
+ fset=_set_data,
+ doc="Sequence as a string (DEPRECATED)")
+
+ def __repr__(self):
+ """Returns a (truncated) representation of the sequence for debugging."""
+ if len(self) > 60 :
+ #Shows the last three letters as it is often useful to see if there
+ #is a stop codon at the end of a sequence.
+ #Note total length is 54+3+3=60
+ return "%s('%s...%s', %s)" % (self.__class__.__name__,
+ str(self)[:54], str(self)[-3:],
+ repr(self.alphabet))
+ else :
+ return "%s(%s, %s)" % (self.__class__.__name__,
+ repr(self.data),
+ repr(self.alphabet))
+ def __str__(self):
+ """Returns the full sequence as a python string.
+
+ Note that Biopython 1.44 and earlier would give a truncated
+ version of repr(my_seq) for str(my_seq). If you are writing code
+ which need to be backwards compatible with old Biopython, you
+ should continue to use my_seq.tostring() rather than str(my_seq).
+ """
+ return self._data
+
+ """
+ TODO - Work out why this breaks test_Restriction.py
+ (Comparing Seq objects would be nice to have. May need to think about
+ hashes and the in operator for when have list/dictionary of Seq objects...)
+ def __cmp__(self, other):
+ if hasattr(other, "alphabet") :
+ #other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet,
+ other.alphabet]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other.alphabet)))
+ #They should be the same sequence type (or one of them is generic)
+ return cmp(str(self), str(other))
+ elif isinstance(other, basestring) :
+ return cmp(str(self), other)
+ else :
+ raise TypeError
+ """
+
+ def __len__(self): return len(self._data) # Seq API requirement
+
+ def __getitem__(self, index) : # Seq API requirement
+ #Note since Python 2.0, __getslice__ is deprecated
+ #and __getitem__ is used instead.
+ #See http://docs.python.org/ref/sequence-methods.html
+ if isinstance(index, int) :
+ #Return a single letter as a string
+ return self._data[index]
+ else :
+ #Return the (sub)sequence as another Seq object
+ return Seq(self._data[index], self.alphabet)
+
+ def __add__(self, other):
+ """Add another sequence or string to this sequence."""
+ if hasattr(other, "alphabet") :
+ #other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet,
+ other.alphabet]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other.alphabet)))
+ #They should be the same sequence type (or one of them is generic)
+ a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
+ return self.__class__(str(self) + str(other), a)
+ elif isinstance(other, basestring) :
+ #other is a plain string - use the current alphabet
+ return self.__class__(str(self) + other, self.alphabet)
+ else :
+ raise TypeError
+
+ def __radd__(self, other):
+ if hasattr(other, "alphabet") :
+ #other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet,
+ other.alphabet]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other.alphabet)))
+ #They should be the same sequence type (or one of them is generic)
+ a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
+ return self.__class__(str(other) + str(self), a)
+ elif isinstance(other, basestring) :
+ #other is a plain string - use the current alphabet
+ return self.__class__(other + str(self), self.alphabet)
+ else :
+ raise TypeError
+
+ def tostring(self): # Seq API requirement
+ """Returns the full sequence as a python string.
+
+ Although not formally deprecated, you are now encouraged to use
+ str(my_seq) instead of my_seq.tostring()."""
+ return str(self)
+
+ def tomutable(self): # Needed? Or use a function?
+ """Returns the full sequence as a MutableSeq object.
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> my_seq = Seq("MKQHKAMIVALIVICITAVVAAL",
+ ... IUPAC.protein)
+ >>> my_seq
+ Seq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
+ >>> my_seq.tomutable()
+ MutableSeq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
+
+ Note that the alphabet is preserved.
+ """
+ return MutableSeq(str(self), self.alphabet)
+
+ def _get_seq_str_and_check_alphabet(self, other_sequence) :
+ """string/Seq/MutableSeq to string, checking alphabet (PRIVATE).
+
+ For a string argument, returns the string.
+
+ For a Seq or MutableSeq, it checks the alphabet is compatible
+ (raising an exception if it isn't), and then returns a string.
+ """
+ try :
+ other_alpha = other_sequence.alphabet
+ except AttributeError :
+ #Assume other_sequence is a string
+ return other_sequence
+
+ #Other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet, other_alpha]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other_alpha)))
+ #Return as a string
+ return str(other_sequence)
+
+ def count(self, sub, start=0, end=sys.maxint):
+ """Non-overlapping count method, like that of a python string.
+
+ This behaves like the python string method of the same name,
+ which does a non-overlapping count!
+
+ Returns an integer, the number of occurrences of substring
+ argument sub in the (sub)sequence given by [start:end].
+ Optional arguments start and end are interpreted as in slice
+ notation.
+
+ Arguments:
+ - sub - a string or another Seq object to look for
+ - start - optional integer, slice start
+ - end - optional integer, slice end
+
+ e.g.
+
+ >>> from Bio.Seq import Seq
+ >>> my_seq = Seq("AAAATGA")
+ >>> print my_seq.count("A")
+ 5
+ >>> print my_seq.count("ATG")
+ 1
+ >>> print my_seq.count(Seq("AT"))
+ 1
+ >>> print my_seq.count("AT", 2, -1)
+ 1
+
+ HOWEVER, please note because python strings and Seq objects (and
+ MutableSeq objects) do a non-overlapping search, this may not give
+ the answer you expect:
+
+ >>> "AAAA".count("AA")
+ 2
+ >>> print Seq("AAAA").count("AA")
+ 2
+
+ A non-overlapping search would give the answer as three!
+ """
+ #If it has one, check the alphabet:
+ sub_str = self._get_seq_str_and_check_alphabet(sub)
+ return str(self).count(sub_str, start, end)
+
+ def find(self, sub, start=0, end=sys.maxint):
+ """Find method, like that of a python string.
+
+ This behaves like the python string method of the same name.
+
+ Returns an integer, the index of the first occurrence of substring
+ argument sub in the (sub)sequence given by [start:end].
+
+ Arguments:
+ - sub - a string or another Seq object to look for
+ - start - optional integer, slice start
+ - end - optional integer, slice end
+
+ Returns -1 if the subsequence is NOT found.
+
+ e.g. Locating the first typical start codon, AUG, in an RNA sequence:
+
+ >>> from Bio.Seq import Seq
+ >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
+ >>> my_rna.find("AUG")
+ 3
+ """
+ #If it has one, check the alphabet:
+ sub_str = self._get_seq_str_and_check_alphabet(sub)
+ return str(self).find(sub_str, start, end)
+
+ def rfind(self, sub, start=0, end=sys.maxint):
+ """Find from right method, like that of a python string.
+
+ This behaves like the python string method of the same name.
+
+ Returns an integer, the index of the last (right most) occurrence of
+ substring argument sub in the (sub)sequence given by [start:end].
+
+ Arguments:
+ - sub - a string or another Seq object to look for
+ - start - optional integer, slice start
+ - end - optional integer, slice end
+
+ Returns -1 if the subsequence is NOT found.
+
+ e.g. Locating the last typical start codon, AUG, in an RNA sequence:
+
+ >>> from Bio.Seq import Seq
+ >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
+ >>> my_rna.rfind("AUG")
+ 15
+ """
+ #If it has one, check the alphabet:
+ sub_str = self._get_seq_str_and_check_alphabet(sub)
+ return str(self).rfind(sub_str, start, end)
+
+ def startswith(self, prefix, start=0, end=sys.maxint) :
+ """Does the Seq start with the given prefix? Returns True/False.
+
+ This behaves like the python string method of the same name.
+
+ Return True if the sequence starts with the specified prefix
+ (a string or another Seq object), False otherwise.
+ With optional start, test sequence beginning at that position.
+ With optional end, stop comparing sequence at that position.
+ prefix can also be a tuple of strings to try. e.g.
+
+ >>> from Bio.Seq import Seq
+ >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
+ >>> my_rna.startswith("GUC")
+ True
+ >>> my_rna.startswith("AUG")
+ False
+ >>> my_rna.startswith("AUG", 3)
+ True
+ >>> my_rna.startswith(("UCC","UCA","UCG"),1)
+ True
+ """
+ #If it has one, check the alphabet:
+ if isinstance(prefix, tuple) :
+ #TODO - Once we drop support for Python 2.4, instead of this
+ #loop offload to the string method (requires Python 2.5+).
+ #Check all the alphabets first...
+ prefix_strings = [self._get_seq_str_and_check_alphabet(p) \
+ for p in prefix]
+ for prefix_str in prefix_strings :
+ if str(self).startswith(prefix_str, start, end) :
+ return True
+ return False
+ else :
+ prefix_str = self._get_seq_str_and_check_alphabet(prefix)
+ return str(self).startswith(prefix_str, start, end)
+
+ def endswith(self, suffix, start=0, end=sys.maxint) :
+ """Does the Seq end with the given suffix? Returns True/False.
+
+ This behaves like the python string method of the same name.
+
+ Return True if the sequence ends with the specified suffix
+ (a string or another Seq object), False otherwise.
+ With optional start, test sequence beginning at that position.
+ With optional end, stop comparing sequence at that position.
+ suffix can also be a tuple of strings to try. e.g.
+
+ >>> from Bio.Seq import Seq
+ >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
+ >>> my_rna.endswith("UUG")
+ True
+ >>> my_rna.endswith("AUG")
+ False
+ >>> my_rna.endswith("AUG", 0, 18)
+ True
+ >>> my_rna.endswith(("UCC","UCA","UUG"))
+ True
+ """
+ #If it has one, check the alphabet:
+ if isinstance(suffix, tuple) :
+ #TODO - Once we drop support for Python 2.4, instead of this
+ #loop offload to the string method (requires Python 2.5+).
+ #Check all the alphabets first...
+ suffix_strings = [self._get_seq_str_and_check_alphabet(p) \
+ for p in suffix]
+ for suffix_str in suffix_strings :
+ if str(self).endswith(suffix_str, start, end) :
+ return True
+ return False
+ else :
+ suffix_str = self._get_seq_str_and_check_alphabet(suffix)
+ return str(self).endswith(suffix_str, start, end)
+
+
+ def split(self, sep=None, maxsplit=-1) :
+ """Split method, like that of a python string.
+
+ This behaves like the python string method of the same name.
+
+ Return a list of the 'words' in the string (as Seq objects),
+ using sep as the delimiter string. If maxsplit is given, at
+ most maxsplit splits are done. If maxsplit is ommited, all
+ splits are made.
+
+ Following the python string method, sep will by default be any
+ white space (tabs, spaces, newlines) but this is unlikely to
+ apply to biological sequences.
+
+ e.g.
+
+ >>> from Bio.Seq import Seq
+ >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
+ >>> my_aa = my_rna.translate()
+ >>> my_aa
+ Seq('VMAIVMGR*KGAR*L', HasStopCodon(ExtendedIUPACProtein(), '*'))
+ >>> my_aa.split("*")
+ [Seq('VMAIVMGR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('KGAR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('L', HasStopCodon(ExtendedIUPACProtein(), '*'))]
+ >>> my_aa.split("*",1)
+ [Seq('VMAIVMGR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('KGAR*L', HasStopCodon(ExtendedIUPACProtein(), '*'))]
+
+ See also the rsplit method:
+
+ >>> my_aa.rsplit("*",1)
+ [Seq('VMAIVMGR*KGAR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('L', HasStopCodon(ExtendedIUPACProtein(), '*'))]
+ """
+ #If it has one, check the alphabet:
+ sep_str = self._get_seq_str_and_check_alphabet(sep)
+ #TODO - If the sep is the defined stop symbol, or gap char,
+ #should we adjust the alphabet?
+ return [Seq(part, self.alphabet) \
+ for part in str(self).split(sep_str, maxsplit)]
+
+ def rsplit(self, sep=None, maxsplit=-1) :
+ """Right split method, like that of a python string.
+
+ This behaves like the python string method of the same name.
+
+ Return a list of the 'words' in the string (as Seq objects),
+ using sep as the delimiter string. If maxsplit is given, at
+ most maxsplit splits are done COUNTING FROM THE RIGHT.
+ If maxsplit is ommited, all splits are made.
+
+ Following the python string method, sep will by default be any
+ white space (tabs, spaces, newlines) but this is unlikely to
+ apply to biological sequences.
+
+ e.g. print my_seq.rsplit("*",1)
+
+ See also the split method.
+ """
+ #If it has one, check the alphabet:
+ sep_str = self._get_seq_str_and_check_alphabet(sep)
+ try :
+ return [Seq(part, self.alphabet) \
+ for part in str(self).rsplit(sep_str, maxsplit)]
+ except AttributeError :
+ #Python 2.3 doesn't have a string rsplit method, which we can
+ #word around by reversing the sequence, using (left) split,
+ #and then reversing the answer. Not very efficient!
+ words = [Seq(word[::-1], self.alphabet) for word \
+ in str(self)[::-1].split(sep_str[::-1], maxsplit)]
+ words.reverse()
+ return words
+
+ def strip(self, chars=None) :
+ """Returns a new Seq object with leading and trailing ends stripped.
+
+ This behaves like the python string method of the same name.
+
+ Optional argument chars defines which characters to remove. If
+ ommitted or None (default) then as for the python string method,
+ this defaults to removing any white space.
+
+ e.g. print my_seq.strip("-")
+
+ See also the lstrip and rstrip methods.
+ """
+ #If it has one, check the alphabet:
+ strip_str = self._get_seq_str_and_check_alphabet(chars)
+ return Seq(str(self).strip(strip_str), self.alphabet)
+
+ def lstrip(self, chars=None) :
+ """Returns a new Seq object with leading (left) end stripped.
+
+ This behaves like the python string method of the same name.
+
+ Optional argument chars defines which characters to remove. If
+ ommitted or None (default) then as for the python string method,
+ this defaults to removing any white space.
+
+ e.g. print my_seq.lstrip("-")
+
+ See also the strip and rstrip methods.
+ """
+ #If it has one, check the alphabet:
+ strip_str = self._get_seq_str_and_check_alphabet(chars)
+ return Seq(str(self).lstrip(strip_str), self.alphabet)
+
+ def rstrip(self, chars=None) :
+ """Returns a new Seq object with trailing (right) end stripped.
+
+ This behaves like the python string method of the same name.
+
+ Optional argument chars defines which characters to remove. If
+ ommitted or None (default) then as for the python string method,
+ this defaults to removing any white space.
+
+ e.g. Removing a nucleotide sequence's polyadenylation (poly-A tail):
+
+ >>> from Bio.Alphabet import IUPAC
+ >>> from Bio.Seq import Seq
+ >>> my_seq = Seq("CGGTACGCTTATGTCACGTAGAAAAAA", IUPAC.unambiguous_dna)
+ >>> my_seq
+ Seq('CGGTACGCTTATGTCACGTAGAAAAAA', IUPACUnambiguousDNA())
+ >>> my_seq.rstrip("A")
+ Seq('CGGTACGCTTATGTCACGTAG', IUPACUnambiguousDNA())
+
+ See also the strip and lstrip methods.
+ """
+ #If it has one, check the alphabet:
+ strip_str = self._get_seq_str_and_check_alphabet(chars)
+ return Seq(str(self).rstrip(strip_str), self.alphabet)
+
+ def complement(self):
+ """Returns the complement sequence. New Seq object.
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> my_dna = Seq("CCCCCGATAG", IUPAC.unambiguous_dna)
+ >>> my_dna
+ Seq('CCCCCGATAG', IUPACUnambiguousDNA())
+ >>> my_dna.complement()
+ Seq('GGGGGCTATC', IUPACUnambiguousDNA())
+
+ You can of course used mixed case sequences,
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import generic_dna
+ >>> my_dna = Seq("CCCCCgatA-GD", generic_dna)
+ >>> my_dna
+ Seq('CCCCCgatA-GD', DNAAlphabet())
+ >>> my_dna.complement()
+ Seq('GGGGGctaT-CH', DNAAlphabet())
+
+ Note in the above example, ambiguous character D denotes
+ G, A or T so its complement is H (for C, T or A).
+
+ Trying to complement a protein sequence raises an exception.
+
+ >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
+ >>> my_protein.complement()
+ Traceback (most recent call last):
+ ...
+ ValueError: Proteins do not have complements!
+ """
+ base = Alphabet._get_base_alphabet(self.alphabet)
+ if isinstance(base, Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins do not have complements!")
+ if isinstance(base, Alphabet.DNAAlphabet) :
+ ttable = _dna_complement_table
+ elif isinstance(base, Alphabet.RNAAlphabet) :
+ ttable = _rna_complement_table
+ elif ('U' in self._data or 'u' in self._data) \
+ and ('T' in self._data or 't' in self._data):
+ #TODO - Handle this cleanly?
+ raise ValueError("Mixed RNA/DNA found")
+ elif 'U' in self._data or 'u' in self._data:
+ ttable = _rna_complement_table
+ else:
+ ttable = _dna_complement_table
+ #Much faster on really long sequences than the previous loop based one.
+ #thx to Michael Palmer, University of Waterloo
+ return Seq(str(self).translate(ttable), self.alphabet)
+
+ def reverse_complement(self):
+ """Returns the reverse complement sequence. New Seq object.
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> my_dna = Seq("CCCCCGATAGNR", IUPAC.ambiguous_dna)
+ >>> my_dna
+ Seq('CCCCCGATAGNR', IUPACAmbiguousDNA())
+ >>> my_dna.reverse_complement()
+ Seq('YNCTATCGGGGG', IUPACAmbiguousDNA())
+
+ Note in the above example, since R = G or A, its complement
+ is Y (which denotes C or T).
+
+ You can of course used mixed case sequences,
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import generic_dna
+ >>> my_dna = Seq("CCCCCgatA-G", generic_dna)
+ >>> my_dna
+ Seq('CCCCCgatA-G', DNAAlphabet())
+ >>> my_dna.reverse_complement()
+ Seq('C-TatcGGGGG', DNAAlphabet())
+
+ Trying to complement a protein sequence raises an exception:
+
+ >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
+ >>> my_protein.reverse_complement()
+ Traceback (most recent call last):
+ ...
+ ValueError: Proteins do not have complements!
+ """
+ #Use -1 stride/step to reverse the complement
+ return self.complement()[::-1]
+
+ def transcribe(self):
+ """Returns the RNA sequence from a DNA sequence. New Seq object.
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG",
+ ... IUPAC.unambiguous_dna)
+ >>> coding_dna
+ Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())
+ >>> coding_dna.transcribe()
+ Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())
+
+ Trying to transcribe a protein or RNA sequence raises an exception:
+
+ >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
+ >>> my_protein.transcribe()
+ Traceback (most recent call last):
+ ...
+ ValueError: Proteins cannot be transcribed!
+ """
+ base = Alphabet._get_base_alphabet(self.alphabet)
+ if isinstance(base, Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins cannot be transcribed!")
+ if isinstance(base, Alphabet.RNAAlphabet) :
+ raise ValueError("RNA cannot be transcribed!")
+
+ if self.alphabet==IUPAC.unambiguous_dna:
+ alphabet = IUPAC.unambiguous_rna
+ elif self.alphabet==IUPAC.ambiguous_dna:
+ alphabet = IUPAC.ambiguous_rna
+ else:
+ alphabet = Alphabet.generic_rna
+ return Seq(str(self).replace('T','U').replace('t','u'), alphabet)
+
+ def back_transcribe(self):
+ """Returns the DNA sequence from an RNA sequence. New Seq object.
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> messenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG",
+ ... IUPAC.unambiguous_rna)
+ >>> messenger_rna
+ Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())
+ >>> messenger_rna.back_transcribe()
+ Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())
+
+ Trying to back-transcribe a protein or DNA sequence raises an
+ exception:
+
+ >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
+ >>> my_protein.back_transcribe()
+ Traceback (most recent call last):
+ ...
+ ValueError: Proteins cannot be back transcribed!
+ """
+ base = Alphabet._get_base_alphabet(self.alphabet)
+ if isinstance(base, Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins cannot be back transcribed!")
+ if isinstance(base, Alphabet.DNAAlphabet) :
+ raise ValueError("DNA cannot be back transcribed!")
+
+ if self.alphabet==IUPAC.unambiguous_rna:
+ alphabet = IUPAC.unambiguous_dna
+ elif self.alphabet==IUPAC.ambiguous_rna:
+ alphabet = IUPAC.ambiguous_dna
+ else:
+ alphabet = Alphabet.generic_dna
+ return Seq(str(self).replace("U", "T").replace("u", "t"), alphabet)
+
+ def translate(self, table="Standard", stop_symbol="*", to_stop=False):
+ """Turns a nucleotide sequence into a protein sequence. New Seq object.
+
+ This method will translate DNA or RNA sequences, and those with a
+ nucleotide or generic alphabet. Trying to translate a protein
+ sequence raises an exception.
+
+ Arguments:
+ - table - Which codon table to use? This can be either a name
+ (string) or an NCBI identifier (integer). This defaults
+ to the "Standard" table.
+ - stop_symbol - Single character string, what to use for terminators.
+ This defaults to the asterisk, "*".
+ - to_stop - Boolean, defaults to False meaning do a full translation
+ continuing on past any stop codons (translated as the
+ specified stop_symbol). If True, translation is
+ terminated at the first in frame stop codon (and the
+ stop_symbol is not appended to the returned protein
+ sequence).
+
+ e.g. Using the standard table:
+
+ >>> coding_dna = Seq("GTGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
+ >>> coding_dna.translate()
+ Seq('VAIVMGR*KGAR*', HasStopCodon(ExtendedIUPACProtein(), '*'))
+ >>> coding_dna.translate(stop_symbol="@")
+ Seq('VAIVMGR@KGAR@', HasStopCodon(ExtendedIUPACProtein(), '@'))
+ >>> coding_dna.translate(to_stop=True)
+ Seq('VAIVMGR', ExtendedIUPACProtein())
+
+ Now using NCBI table 2, where TGA is not a stop codon:
+
+ >>> coding_dna.translate(table=2)
+ Seq('VAIVMGRWKGAR*', HasStopCodon(ExtendedIUPACProtein(), '*'))
+ >>> coding_dna.translate(table=2, to_stop=True)
+ Seq('VAIVMGRWKGAR', ExtendedIUPACProtein())
+
+ If the sequence has no in-frame stop codon, then the to_stop argument
+ has no effect:
+
+ >>> coding_dna2 = Seq("TTGGCCATTGTAATGGGCCGC")
+ >>> coding_dna2.translate()
+ Seq('LAIVMGR', ExtendedIUPACProtein())
+ >>> coding_dna2.translate(to_stop=True)
+ Seq('LAIVMGR', ExtendedIUPACProtein())
+
+ NOTE - Ambiguous codons like "TAN" or "NNN" could be an amino acid
+ or a stop codon. These are translated as "X". Any invalid codon
+ (e.g. "TA?" or "T-A") will throw a TranslationError.
+
+ NOTE - Does NOT support gapped sequences.
+
+ NOTE - This does NOT behave like the python string's translate
+ method. For that use str(my_seq).translate(...) instead.
+ """
+ try:
+ table_id = int(table)
+ except ValueError:
+ table_id = None
+ if isinstance(table, str) and len(table)==256 :
+ raise ValueError("The Seq object translate method DOES NOT take " \
+ + "a 256 character string mapping table like " \
+ + "the python string object's translate method. " \
+ + "Use str(my_seq).translate(...) instead.")
+ if isinstance(Alphabet._get_base_alphabet(self.alphabet),
+ Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins cannot be translated!")
+ if self.alphabet==IUPAC.unambiguous_dna:
+ #Will use standard IUPAC protein alphabet, no need for X
+ if table_id is None:
+ codon_table = CodonTable.unambiguous_dna_by_name[table]
+ else:
+ codon_table = CodonTable.unambiguous_dna_by_id[table_id]
+ #elif self.alphabet==IUPAC.ambiguous_dna:
+ # if table_id is None:
+ # codon_table = CodonTable.ambiguous_dna_by_name[table]
+ # else:
+ # codon_table = CodonTable.ambiguous_dna_by_id[table_id]
+ elif self.alphabet==IUPAC.unambiguous_rna:
+ #Will use standard IUPAC protein alphabet, no need for X
+ if table_id is None:
+ codon_table = CodonTable.unambiguous_rna_by_name[table]
+ else:
+ codon_table = CodonTable.unambiguous_rna_by_id[table_id]
+ #elif self.alphabet==IUPAC.ambiguous_rna:
+ # if table_id is None:
+ # codon_table = CodonTable.ambiguous_rna_by_name[table]
+ # else:
+ # codon_table = CodonTable.ambiguous_rna_by_id[table_id]
+ else:
+ #This will use the extend IUPAC protein alphabet with X etc.
+ #The same table can be used for RNA or DNA (we use this for
+ #translating strings).
+ if table_id is None:
+ codon_table = CodonTable.ambiguous_generic_by_name[table]
+ else:
+ codon_table = CodonTable.ambiguous_generic_by_id[table_id]
+ protein = _translate_str(str(self), codon_table, stop_symbol, to_stop)
+ if stop_symbol in protein :
+ alphabet = Alphabet.HasStopCodon(codon_table.protein_alphabet,
+ stop_symbol = stop_symbol)
+ else :
+ alphabet = codon_table.protein_alphabet
+ return Seq(protein, alphabet)
+
+class UnknownSeq(Seq):
+ """A read-only sequence object of known length but unknown contents.
+
+ If you have an unknown sequence, you can represent this with a normal
+ Seq object, for example:
+
+ >>> my_seq = Seq("N"*5)
+ >>> my_seq
+ Seq('NNNNN', Alphabet())
+ >>> len(my_seq)
+ 5
+ >>> print my_seq
+ NNNNN
+
+ However, this is rather wasteful of memory (especially for large
+ sequences), which is where this class is most usefull:
+
+ >>> unk_five = UnknownSeq(5)
+ >>> unk_five
+ UnknownSeq(5, alphabet = Alphabet(), character = '?')
+ >>> len(unk_five)
+ 5
+ >>> print(unk_five)
+ ?????
+
+ You can add unknown sequence together, provided their alphabets and
+ characters are compatible, and get another memory saving UnknownSeq:
+
+ >>> unk_four = UnknownSeq(4)
+ >>> unk_four
+ UnknownSeq(4, alphabet = Alphabet(), character = '?')
+ >>> unk_four + unk_five
+ UnknownSeq(9, alphabet = Alphabet(), character = '?')
+
+ If the alphabet or characters don't match up, the addition gives an
+ ordinary Seq object:
+
+ >>> unk_nnnn = UnknownSeq(4, character = "N")
+ >>> unk_nnnn
+ UnknownSeq(4, alphabet = Alphabet(), character = 'N')
+ >>> unk_nnnn + unk_four
+ Seq('NNNN????', Alphabet())
+
+ Combining with a real Seq gives a new Seq object:
+
+ >>> known_seq = Seq("ACGT")
+ >>> unk_four + known_seq
+ Seq('????ACGT', Alphabet())
+ >>> known_seq + unk_four
+ Seq('ACGT????', Alphabet())
+ """
+ def __init__(self, length, alphabet = Alphabet.generic_alphabet, character = None) :
+ """Create a new UnknownSeq object.
+
+ If character is ommited, it is determed from the alphabet, "N" for
+ nucleotides, "X" for proteins, and "?" otherwise.
+ """
+ self._length = int(length)
+ if self._length < 0 :
+ #TODO - Block zero length UnknownSeq? You can just use a Seq!
+ raise ValueError("Length must not be negative.")
+ self.alphabet = alphabet
+ if character :
+ if len(character) != 1 :
+ raise ValueError("character argument should be a single letter string.")
+ self._character = character
+ else :
+ base = Alphabet._get_base_alphabet(alphabet)
+ #TODO? Check the case of the letters in the alphabet?
+ #We may have to use "n" instead of "N" etc.
+ if isinstance(base, Alphabet.NucleotideAlphabet) :
+ self._character = "N"
+ elif isinstance(base, Alphabet.ProteinAlphabet) :
+ self._character = "X"
+ else :
+ self._character = "?"
+
+ def __len__(self) :
+ """Returns the stated length of the unknown sequence."""
+ return self._length
+
+ def __str__(self) :
+ """Returns the unknown sequence as full string of the given length."""
+ return self._character * self._length
+
+ def __repr__(self):
+ return "UnknownSeq(%i, alphabet = %s, character = %s)" \
+ % (self._length, repr(self.alphabet), repr(self._character))
+
+ def __add__(self, other) :
+ if isinstance(other, UnknownSeq) \
+ and other._character == self._character :
+ #TODO - Check the alphabets match
+ return UnknownSeq(len(self)+len(other),
+ self.alphabet, self._character)
+ #Offload to the base class...
+ return Seq(str(self), self.alphabet) + other
+
+ def __radd__(self, other) :
+ if isinstance(other, UnknownSeq) \
+ and other._character == self._character :
+ #TODO - Check the alphabets match
+ return UnknownSeq(len(self)+len(other),
+ self.alphabet, self._character)
+ #Offload to the base class...
+ return other + Seq(str(self), self.alphabet)
+
+ def __getitem__(self, index):
+ if isinstance(index, int) :
+ #TODO - Check the bounds without wasting memory
+ return str(self)[index]
+ else :
+ #TODO - Work out the length without wasting memory
+ return UnknownSeq(len(("#"*self._length)[index]),
+ self.alphabet, self._character)
+
+ def count(self, sub, start=0, end=sys.maxint):
+ """Non-overlapping count method, like that of a python string.
+
+ This behaves like the python string (and Seq object) method of the
+ same name, which does a non-overlapping count!
+
+ Returns an integer, the number of occurrences of substring
+ argument sub in the (sub)sequence given by [start:end].
+ Optional arguments start and end are interpreted as in slice
+ notation.
+
+ Arguments:
+ - sub - a string or another Seq object to look for
+ - start - optional integer, slice start
+ - end - optional integer, slice end
+
+ >>> "NNNN".count("N")
+ 4
+ >>> Seq("NNNN").count("N")
+ 4
+ >>> UnknownSeq(4, character="N").count("N")
+ 4
+ >>> UnknownSeq(4, character="N").count("A")
+ 0
+ >>> UnknownSeq(4, character="N").count("AA")
+ 0
+
+ HOWEVER, please note because that python strings and Seq objects (and
+ MutableSeq objects) do a non-overlapping search, this may not give
+ the answer you expect:
+
+ >>> UnknownSeq(4, character="N").count("NN")
+ 2
+ >>> UnknownSeq(4, character="N").count("NNN")
+ 1
+ """
+ sub_str = self._get_seq_str_and_check_alphabet(sub)
+ if len(sub_str) == 1 :
+ if str(sub_str) == self._character :
+ if start==0 and end >= self._length :
+ return self._length
+ else :
+ #This could be done more cleverly...
+ return str(self).count(sub_str, start, end)
+ else :
+ return 0
+ else :
+ if set(sub_str) == set(self._character) :
+ if start==0 and end >= self._length :
+ return self._length // len(sub_str)
+ else :
+ #This could be done more cleverly...
+ return str(self).count(sub_str, start, end)
+ else :
+ return 0
+
+ def complement(self) :
+ """The complement of an unknown nucleotide equals itself.
+
+ >>> my_nuc = UnknownSeq(8)
+ >>> my_nuc
+ UnknownSeq(8, alphabet = Alphabet(), character = '?')
+ >>> print my_nuc
+ ????????
+ >>> my_nuc.complement()
+ UnknownSeq(8, alphabet = Alphabet(), character = '?')
+ >>> print my_nuc.complement()
+ ????????
+ """
+ if isinstance(Alphabet._get_base_alphabet(self.alphabet),
+ Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins do not have complements!")
+ return self
+
+ def reverse_complement(self) :
+ """The reverse complement of an unknown nucleotide equals itself.
+
+ >>> my_nuc = UnknownSeq(10)
+ >>> my_nuc
+ UnknownSeq(10, alphabet = Alphabet(), character = '?')
+ >>> print my_nuc
+ ??????????
+ >>> my_nuc.reverse_complement()
+ UnknownSeq(10, alphabet = Alphabet(), character = '?')
+ >>> print my_nuc.reverse_complement()
+ ??????????
+ """
+ if isinstance(Alphabet._get_base_alphabet(self.alphabet),
+ Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins do not have complements!")
+ return self
+
+ def transcribe(self) :
+ """Returns unknown RNA sequence from an unknown DNA sequence.
+
+ >>> my_dna = UnknownSeq(10, character="N")
+ >>> my_dna
+ UnknownSeq(10, alphabet = Alphabet(), character = 'N')
+ >>> print my_dna
+ NNNNNNNNNN
+ >>> my_rna = my_dna.transcribe()
+ >>> my_rna
+ UnknownSeq(10, alphabet = RNAAlphabet(), character = 'N')
+ >>> print my_rna
+ NNNNNNNNNN
+ """
+ #Offload the alphabet stuff
+ s = Seq(self._character, self.alphabet).transcribe()
+ return UnknownSeq(self._length, s.alphabet, self._character)
+
+ def back_transcribe(self) :
+ """Returns unknown DNA sequence from an unknown RNA sequence.
+
+ >>> my_rna = UnknownSeq(20, character="N")
+ >>> my_rna
+ UnknownSeq(20, alphabet = Alphabet(), character = 'N')
+ >>> print my_rna
+ NNNNNNNNNNNNNNNNNNNN
+ >>> my_dna = my_rna.back_transcribe()
+ >>> my_dna
+ UnknownSeq(20, alphabet = DNAAlphabet(), character = 'N')
+ >>> print my_dna
+ NNNNNNNNNNNNNNNNNNNN
+ """
+ #Offload the alphabet stuff
+ s = Seq(self._character, self.alphabet).back_transcribe()
+ return UnknownSeq(self._length, s.alphabet, self._character)
+
+ def translate(self, **kwargs) :
+ """Translate an unknown nucleotide sequence into an unknown protein.
+
+ e.g.
+
+ >>> my_seq = UnknownSeq(11, character="N")
+ >>> print my_seq
+ NNNNNNNNNNN
+ >>> my_protein = my_seq.translate()
+ >>> my_protein
+ UnknownSeq(3, alphabet = ProteinAlphabet(), character = 'X')
+ >>> print my_protein
+ XXX
+
+ In comparison, using a normal Seq object:
+
+ >>> my_seq = Seq("NNNNNNNNNNN")
+ >>> print my_seq
+ NNNNNNNNNNN
+ >>> my_protein = my_seq.translate()
+ >>> my_protein
+ Seq('XXX', ExtendedIUPACProtein())
+ >>> print my_protein
+ XXX
+
+ """
+ if isinstance(Alphabet._get_base_alphabet(self.alphabet),
+ Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins cannot be translated!")
+ return UnknownSeq(self._length//3, Alphabet.generic_protein, "X")
+
+
+class MutableSeq(object):
+ """An editable sequence object (with an alphabet).
+
+ Unlike normal python strings and our basic sequence object (the Seq class)
+ which are immuatable, the MutableSeq lets you edit the sequence in place.
+ However, this means you cannot use a MutableSeq object as a dictionary key.
+
+ >>> from Bio.Seq import MutableSeq
+ >>> from Bio.Alphabet import generic_dna
+ >>> my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
+ >>> my_seq
+ MutableSeq('ACTCGTCGTCG', DNAAlphabet())
+ >>> my_seq[5]
+ 'T'
+ >>> my_seq[5] = "A"
+ >>> my_seq
+ MutableSeq('ACTCGACGTCG', DNAAlphabet())
+ >>> my_seq[5]
+ 'A'
+ >>> my_seq[5:8] = "NNN"
+ >>> my_seq
+ MutableSeq('ACTCGNNNTCG', DNAAlphabet())
+ >>> len(my_seq)
+ 11
+
+ Note that the MutableSeq object does not support as many string-like
+ or biological methods as the Seq object.
+ """
+ def __init__(self, data, alphabet = Alphabet.generic_alphabet):
+ if type(data) == type(""):
+ self.data = array.array("c", data)
+ else:
+ self.data = data # assumes the input is an array
+ self.alphabet = alphabet
+
+ def __repr__(self):
+ """Returns a (truncated) representation of the sequence for debugging."""
+ if len(self) > 60 :
+ #Shows the last three letters as it is often useful to see if there
+ #is a stop codon at the end of a sequence.
+ #Note total length is 54+3+3=60
+ return "%s('%s...%s', %s)" % (self.__class__.__name__,
+ str(self[:54]), str(self[-3:]),
+ repr(self.alphabet))
+ else :
+ return "%s('%s', %s)" % (self.__class__.__name__,
+ str(self),
+ repr(self.alphabet))
+
+ def __str__(self):
+ """Returns the full sequence as a python string.
+
+ Note that Biopython 1.44 and earlier would give a truncated
+ version of repr(my_seq) for str(my_seq). If you are writing code
+ which needs to be backwards compatible with old Biopython, you
+ should continue to use my_seq.tostring() rather than str(my_seq).
+ """
+ #See test_GAQueens.py for an historic usage of a non-string alphabet!
+ return "".join(self.data)
+
+ def __cmp__(self, other):
+ """Compare the sequence for to another sequence or a string.
+
+ If compared to another sequence the alphabets must be compatible.
+ Comparing DNA to RNA, or Nucleotide to Protein will raise an
+ exception.
+
+ Otherwise only the sequence itself is compared, not the precise
+ alphabet.
+
+ This method indirectly supports ==, < , etc."""
+ if hasattr(other, "alphabet") :
+ #other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet,
+ other.alphabet]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other.alphabet)))
+ #They should be the same sequence type (or one of them is generic)
+ if isinstance(other, MutableSeq):
+ #See test_GAQueens.py for an historic usage of a non-string
+ #alphabet! Comparing the arrays supports this.
+ return cmp(self.data, other.data)
+ else :
+ return cmp(str(self), str(other))
+ elif isinstance(other, basestring) :
+ return cmp(str(self), other)
+ else :
+ raise TypeError
+
+ def __len__(self): return len(self.data)
+
+ def __getitem__(self, index) :
+ #Note since Python 2.0, __getslice__ is deprecated
+ #and __getitem__ is used instead.
+ #See http://docs.python.org/ref/sequence-methods.html
+ if isinstance(index, int) :
+ #Return a single letter as a string
+ return self.data[index]
+ else :
+ #Return the (sub)sequence as another Seq object
+ return MutableSeq(self.data[index], self.alphabet)
+
+ def __setitem__(self, index, value):
+ #Note since Python 2.0, __setslice__ is deprecated
+ #and __setitem__ is used instead.
+ #See http://docs.python.org/ref/sequence-methods.html
+ if isinstance(index, int) :
+ #Replacing a single letter with a new string
+ self.data[index] = value
+ else :
+ #Replacing a sub-sequence
+ if isinstance(value, MutableSeq):
+ self.data[index] = value.data
+ elif isinstance(value, type(self.data)):
+ self.data[index] = value
+ else:
+ self.data[index] = array.array("c", str(value))
+
+ def __delitem__(self, index):
+ #Note since Python 2.0, __delslice__ is deprecated
+ #and __delitem__ is used instead.
+ #See http://docs.python.org/ref/sequence-methods.html
+
+ #Could be deleting a single letter, or a slice
+ del self.data[index]
+
+ def __add__(self, other):
+ """Add another sequence or string to this sequence.
+
+ Returns a new MutableSeq object."""
+ if hasattr(other, "alphabet") :
+ #other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet,
+ other.alphabet]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other.alphabet)))
+ #They should be the same sequence type (or one of them is generic)
+ a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
+ if isinstance(other, MutableSeq):
+ #See test_GAQueens.py for an historic usage of a non-string
+ #alphabet! Adding the arrays should support this.
+ return self.__class__(self.data + other.data, a)
+ else :
+ return self.__class__(str(self) + str(other), a)
+ elif isinstance(other, basestring) :
+ #other is a plain string - use the current alphabet
+ return self.__class__(str(self) + str(other), self.alphabet)
+ else :
+ raise TypeError
+
+ def __radd__(self, other):
+ if hasattr(other, "alphabet") :
+ #other should be a Seq or a MutableSeq
+ if not Alphabet._check_type_compatible([self.alphabet,
+ other.alphabet]) :
+ raise TypeError("Incompatable alphabets %s and %s" \
+ % (repr(self.alphabet), repr(other.alphabet)))
+ #They should be the same sequence type (or one of them is generic)
+ a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
+ if isinstance(other, MutableSeq):
+ #See test_GAQueens.py for an historic usage of a non-string
+ #alphabet! Adding the arrays should support this.
+ return self.__class__(other.data + self.data, a)
+ else :
+ return self.__class__(str(other) + str(self), a)
+ elif isinstance(other, basestring) :
+ #other is a plain string - use the current alphabet
+ return self.__class__(str(other) + str(self), self.alphabet)
+ else :
+ raise TypeError
+
+ def append(self, c):
+ self.data.append(c)
+
+ def insert(self, i, c):
+ self.data.insert(i, c)
+
+ def pop(self, i = (-1)):
+ c = self.data[i]
+ del self.data[i]
+ return c
+
+ def remove(self, item):
+ for i in range(len(self.data)):
+ if self.data[i] == item:
+ del self.data[i]
+ return
+ raise ValueError("MutableSeq.remove(x): x not in list")
+
+ def count(self, sub, start=0, end=sys.maxint):
+ """Non-overlapping count method, like that of a python string.
+
+ This behaves like the python string method of the same name,
+ which does a non-overlapping count!
+
+ Returns an integer, the number of occurrences of substring
+ argument sub in the (sub)sequence given by [start:end].
+ Optional arguments start and end are interpreted as in slice
+ notation.
+
+ Arguments:
+ - sub - a string or another Seq object to look for
+ - start - optional integer, slice start
+ - end - optional integer, slice end
+
+ e.g.
+
+ >>> from Bio.Seq import MutableSeq
+ >>> my_mseq = MutableSeq("AAAATGA")
+ >>> print my_mseq.count("A")
+ 5
+ >>> print my_mseq.count("ATG")
+ 1
+ >>> print my_mseq.count(Seq("AT"))
+ 1
+ >>> print my_mseq.count("AT", 2, -1)
+ 1
+
+ HOWEVER, please note because that python strings, Seq objects and
+ MutableSeq objects do a non-overlapping search, this may not give
+ the answer you expect:
+
+ >>> "AAAA".count("AA")
+ 2
+ >>> print MutableSeq("AAAA").count("AA")
+ 2
+
+ A non-overlapping search would give the answer as three!
+ """
+ try :
+ #TODO - Should we check the alphabet?
+ search = sub.tostring()
+ except AttributeError :
+ search = sub
+
+ if not isinstance(search, basestring) :
+ raise TypeError("expected a string, Seq or MutableSeq")
+
+ if len(search) == 1 :
+ #Try and be efficient and work directly from the array.
+ count = 0
+ for c in self.data[start:end]:
+ if c == search: count += 1
+ return count
+ else :
+ #TODO - Can we do this more efficiently?
+ return self.tostring().count(search, start, end)
+
+ def index(self, item):
+ for i in range(len(self.data)):
+ if self.data[i] == item:
+ return i
+ raise ValueError("MutableSeq.index(x): x not in list")
+
+ def reverse(self):
+ """Modify the mutable sequence to reverse itself.
+
+ No return value.
+ """
+ self.data.reverse()
+
+ def complement(self):
+ """Modify the mutable sequence to take on its complement.
+
+ Trying to complement a protein sequence raises an exception.
+
+ No return value.
+ """
+ if isinstance(Alphabet._get_base_alphabet(self.alphabet),
+ Alphabet.ProteinAlphabet) :
+ raise ValueError("Proteins do not have complements!")
+ if self.alphabet in (IUPAC.ambiguous_dna, IUPAC.unambiguous_dna):
+ d = ambiguous_dna_complement
+ elif self.alphabet in (IUPAC.ambiguous_rna, IUPAC.unambiguous_rna):
+ d = ambiguous_rna_complement
+ elif 'U' in self.data and 'T' in self.data :
+ #TODO - Handle this cleanly?
+ raise ValueError("Mixed RNA/DNA found")
+ elif 'U' in self.data:
+ d = ambiguous_rna_complement
+ else:
+ d = ambiguous_dna_complement
+ c = dict([(x.lower(), y.lower()) for x,y in d.iteritems()])
+ d.update(c)
+ self.data = map(lambda c: d[c], self.data)
+ self.data = array.array('c', self.data)
+
+ def reverse_complement(self):
+ """Modify the mutable sequence to take on its reverse complement.
+
+ Trying to reverse complement a protein sequence raises an exception.
+
+ No return value.
+ """
+ self.complement()
+ self.data.reverse()
+
+ ## Sorting a sequence makes no sense.
+ # def sort(self, *args): self.data.sort(*args)
+
+ def extend(self, other):
+ if isinstance(other, MutableSeq):
+ for c in other.data:
+ self.data.append(c)
+ else:
+ for c in other:
+ self.data.append(c)
+
+ def tostring(self):
+ """Returns the full sequence as a python string.
+
+ Although not formally deprecated, you are now encouraged to use
+ str(my_seq) instead of my_seq.tostring().
+
+ Because str(my_seq) will give you the full sequence as a python string,
+ there is often no need to make an explicit conversion. For example,
+
+ print "ID={%s}, sequence={%s}" % (my_name, my_seq)
+
+ On Biopython 1.44 or older you would have to have done this:
+
+ print "ID={%s}, sequence={%s}" % (my_name, my_seq.tostring())
+ """
+ return "".join(self.data)
+
+ def toseq(self):
+ """Returns the full sequence as a new immutable Seq object.
+
+ >>> from Bio.Seq import Seq
+ >>> from Bio.Alphabet import IUPAC
+ >>> my_mseq = MutableSeq("MKQHKAMIVALIVICITAVVAAL", \
+ IUPAC.protein)
+ >>> my_mseq
+ MutableSeq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
+ >>> my_mseq.toseq()
+ Seq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
+
+ Note that the alphabet is preserved.
+ """
+ return Seq("".join(self.data), self.alphabet)
+
+# The transcribe, backward_transcribe, and translate functions are
+# user-friendly versions of the corresponding functions in Bio.Transcribe
+# and Bio.Translate. The functions work both on Seq objects, and on strings.
+
+def transcribe(dna):
+ """Transcribes a DNA sequence into RNA.
+
+ If given a string, returns a new string object.
+
+ Given a Seq or MutableSeq, returns a new Seq object with an RNA alphabet.
+
+ Trying to transcribe a protein or RNA sequence raises an exception.
+
+ e.g.
+
+ >>> transcribe("ACTGN")
+ 'ACUGN'
+ """
+ if isinstance(dna, Seq) :
+ return dna.transcribe()
+ elif isinstance(dna, MutableSeq):
+ return dna.toseq().transcribe()
+ else:
+ return dna.replace('T','U').replace('t','u')
+
+def back_transcribe(rna):
+ """Back-transcribes an RNA sequence into DNA.
+
+ If given a string, returns a new string object.
+
+ Given a Seq or MutableSeq, returns a new Seq object with an RNA alphabet.
+
+ Trying to transcribe a protein or DNA sequence raises an exception.
+
+ e.g.
+
+ >>> back_transcribe("ACUGN")
+ 'ACTGN'
+ """
+ if isinstance(rna, Seq) :
+ return rna.back_transcribe()
+ elif isinstance(rna, MutableSeq):
+ return rna.toseq().back_transcribe()
+ else:
+ return rna.replace('U','T').replace('u','t')
+
+def _translate_str(sequence, table, stop_symbol="*",
+ to_stop=False, pos_stop="X") :
+ """Helper function to translate a nucleotide string (PRIVATE).
+
+ Arguments:
+ - sequence - a string
+ - table - a CodonTable object (NOT a table name or id number)
+ - stop_symbol - a single character string, what to use for terminators.
+ - to_stop - boolean, should translation terminate at the first
+ in frame stop codon? If there is no in-frame stop codon
+ then translation continues to the end.
+ - pos_stop - a single character string for a possible stop codon
+ (e.g. TAN or NNN)
+
+ Returns a string.
+
+ e.g.
+
+ >>> from Bio.Data import CodonTable
+ >>> table = CodonTable.ambiguous_dna_by_id[1]
+ >>> _translate_str("AAA", table)
+ 'K'
+ >>> _translate_str("TAR", table)
+ '*'
+ >>> _translate_str("TAN", table)
+ 'X'
+ >>> _translate_str("TAN", table, pos_stop="@")
+ '@'
+ >>> _translate_str("TA?", table)
+ Traceback (most recent call last):
+ ...
+ TranslationError: Codon 'TA?' is invalid
+ """
+ sequence = sequence.upper()
+ amino_acids = []
+ forward_table = table.forward_table
+ stop_codons = table.stop_codons
+ if table.nucleotide_alphabet.letters is not None :
+ valid_letters = set(table.nucleotide_alphabet.letters.upper())
+ else :
+ #Assume the worst case, ambiguous DNA or RNA:
+ valid_letters = set(IUPAC.ambiguous_dna.letters.upper() + \
+ IUPAC.ambiguous_rna.letters.upper())
+
+ n = len(sequence)
+ for i in xrange(0,n-n%3,3) :
+ codon = sequence[i:i+3]
+ try :
+ amino_acids.append(forward_table[codon])
+ except (KeyError, CodonTable.TranslationError) :
+ #Todo? Treat "---" as a special case (gapped translation)
+ if codon in table.stop_codons :
+ if to_stop : break
+ amino_acids.append(stop_symbol)
+ elif valid_letters.issuperset(set(codon)) :
+ #Possible stop codon (e.g. NNN or TAN)
+ amino_acids.append(pos_stop)
+ else :
+ raise CodonTable.TranslationError(\
+ "Codon '%s' is invalid" % codon)
+ return "".join(amino_acids)
+
+def translate(sequence, table="Standard", stop_symbol="*", to_stop=False):
+ """Translate a nucleotide sequence into amino acids.
+
+ If given a string, returns a new string object. Given a Seq or
+ MutableSeq, returns a Seq object with a protein alphabet.
+
+ Arguments:
+ - table - Which codon table to use? This can be either a name
+ (string) or an NCBI identifier (integer). Defaults
+ to the "Standard" table.
+ - stop_symbol - Single character string, what to use for any
+ terminators, defaults to the asterisk, "*".
+ - to_stop - Boolean, defaults to False meaning do a full
+ translation continuing on past any stop codons
+ (translated as the specified stop_symbol). If
+ True, translation is terminated at the first in
+ frame stop codon (and the stop_symbol is not
+ appended to the returned protein sequence).
+
+ A simple string example using the default (standard) genetic code:
+
+ >>> coding_dna = "GTGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"
+ >>> translate(coding_dna)
+ 'VAIVMGR*KGAR*'
+ >>> translate(coding_dna, stop_symbol="@")
+ 'VAIVMGR@KGAR@'
+ >>> translate(coding_dna, to_stop=True)
+ 'VAIVMGR'
+
+ Now using NCBI table 2, where TGA is not a stop codon:
+
+ >>> translate(coding_dna, table=2)
+ 'VAIVMGRWKGAR*'
+ >>> translate(coding_dna, table=2, to_stop=True)
+ 'VAIVMGRWKGAR'
+
+ Note that if the sequence has no in-frame stop codon, then the to_stop
+ argument has no effect:
+
+ >>> coding_dna2 = "GTGGCCATTGTAATGGGCCGC"
+ >>> translate(coding_dna2)
+ 'VAIVMGR'
+ >>> translate(coding_dna2, to_stop=True)
+ 'VAIVMGR'
+
+ NOTE - Ambiguous codons like "TAN" or "NNN" could be an amino acid
+ or a stop codon. These are translated as "X". Any invalid codon
+ (e.g. "TA?" or "T-A") will throw a TranslationError.
+
+ NOTE - Does NOT support gapped sequences.
+
+ It will however translate either DNA or RNA.
+ """
+ if isinstance(sequence, Seq) :
+ return sequence.translate(table, stop_symbol, to_stop)
+ elif isinstance(sequence, MutableSeq):
+ #Return a Seq object
+ return sequence.toseq().translate(table, stop_symbol, to_stop)
+ else:
+ #Assume its a string, return a string
+ try :
+ codon_table = CodonTable.ambiguous_generic_by_id[int(table)]
+ except ValueError :
+ codon_table = CodonTable.ambiguous_generic_by_name[table]
+ return _translate_str(sequence, codon_table, stop_symbol, to_stop)
+
+def reverse_complement(sequence):
+ """Returns the reverse complement sequence of a nucleotide string.
+
+ If given a string, returns a new string object.
+ Given a Seq or a MutableSeq, returns a new Seq object with the same alphabet.
+
+ Supports unambiguous and ambiguous nucleotide sequences.
+
+ e.g.
+
+ >>> reverse_complement("ACTG-NH")
+ 'DN-CAGT'
+ """
+ if isinstance(sequence, Seq) :
+ #Return a Seq
+ return sequence.reverse_complement()
+ elif isinstance(sequence, MutableSeq) :
+ #Return a Seq
+ #Don't use the MutableSeq reverse_complement method as it is 'in place'.
+ return sequence.toseq().reverse_complement()
+
+ #Assume its a string.
+ #In order to avoid some code duplication, the old code would turn the string
+ #into a Seq, use the reverse_complement method, and convert back to a string.
+ #This worked, but is over five times slower on short sequences!
+ if ('U' in sequence or 'u' in sequence) \
+ and ('T' in sequence or 't' in sequence):
+ raise ValueError("Mixed RNA/DNA found")
+ elif 'U' in sequence or 'u' in sequence:
+ ttable = _rna_complement_table
+ else:
+ ttable = _dna_complement_table
+ return sequence.translate(ttable)[::-1]
+
+def _test():
+ """Run the Bio.Seq module's doctests."""
+ print "Runing doctests..."
+ import doctest
+ doctest.testmod()
+ print "Done"
+
+if __name__ == "__main__":
+ _test()
git://source.jalview.org / jabaws.git / blobdiff