binaries/src/globplot/biopython-1.50/Bio/Seq.py

   1 # Copyright 2000-2002 Brad Chapman.
   2 # Copyright 2004-2005 by M de Hoon.
   3 # Copyright 2007-2009 by Peter Cock.
   4 # All rights reserved.
   5 # This code is part of the Biopython distribution and governed by its
   6 # license.  Please see the LICENSE file that should have been included
   7 # as part of this package.
   8 """Provides objects to represent biological sequences with alphabets.
   9
  10 See also U{http://biopython.org/wiki/Seq} and the chapter in our tutorial:
  11  - U{http://biopython.org/DIST/docs/tutorial/Tutorial.html}
  12  - U{http://biopython.org/DIST/docs/tutorial/Tutorial.pdf}
  13 """
  14 __docformat__ ="epytext en" #Don't just use plain text in epydoc API pages!
  15
  16 import string #for maketrans only
  17 import array
  18 import sys
  19
  20 #TODO - Remove this work around once we drop python 2.3 support
  21 try:
  22     set = set
  23 except NameError:
  24     from sets import Set as set
  25
  26 import Alphabet
  27 from Alphabet import IUPAC
  28 from Data.IUPACData import ambiguous_dna_complement, ambiguous_rna_complement
  29 from Bio.Data import CodonTable
  30
  31 def _maketrans(complement_mapping) :
  32     """Makes a python string translation table (PRIVATE).
  33
  34     Arguments:
  35      - complement_mapping - a dictionary such as ambiguous_dna_complement
  36        and ambiguous_rna_complement from Data.IUPACData.
  37
  38     Returns a translation table (a string of length 256) for use with the
  39     python string's translate method to use in a (reverse) complement.
  40
  41     Compatible with lower case and upper case sequences.
  42
  43     For internal use only.
  44     """
  45     before = ''.join(complement_mapping.keys())
  46     after  = ''.join(complement_mapping.values())
  47     before = before + before.lower()
  48     after  = after + after.lower()
  49     return string.maketrans(before, after)
  50
  51 _dna_complement_table = _maketrans(ambiguous_dna_complement)
  52 _rna_complement_table = _maketrans(ambiguous_rna_complement)
  53
  54 class Seq(object):
  55     """A read-only sequence object (essentially a string with an alphabet).
  56
  57     Like normal python strings, our basic sequence object is immutable.
  58     This prevents you from doing my_seq[5] = "A" for example, but does allow
  59     Seq objects to be used as dictionary keys.
  60
  61     The Seq object provides a number of string like methods (such as count,
  62     find, split and strip), which are alphabet aware where appropriate.
  63
  64     The Seq object also provides some biological methods, such as complement,
  65     reverse_complement, transcribe, back_transcribe and translate (which are
  66     not applicable to sequences with a protein alphabet).
  67     """
  68     def __init__(self, data, alphabet = Alphabet.generic_alphabet):
  69         """Create a Seq object.
  70
  71         Arguments:
  72          - seq      - Sequence, required (string)
  73          - alphabet - Optional argument, an Alphabet object from Bio.Alphabet
  74
  75         You will typically use Bio.SeqIO to read in sequences from files as
  76         SeqRecord objects, whose sequence will be exposed as a Seq object via
  77         the seq property.
  78
  79         However, will often want to create your own Seq objects directly:
  80
  81         >>> from Bio.Seq import Seq
  82         >>> from Bio.Alphabet import IUPAC
  83         >>> my_seq = Seq("MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF",
  84         ...              IUPAC.protein)
  85         >>> my_seq
  86         Seq('MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF', IUPACProtein())
  87         >>> print my_seq
  88         MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF
  89         """
  90         # Enforce string storage
  91         assert (type(data) == type("") or # must use a string
  92                 type(data) == type(u""))  # but can be a unicode string
  93         self._data = data
  94         self.alphabet = alphabet  # Seq API requirement
  95
  96     # A data property is/was a Seq API requirement
  97     def _set_data(self, value) :
  98         #TODO - In the next release, actually raise an exception?
  99         #The Seq object is like a python string, it should be read only!
 100         import warnings
 101         warnings.warn("Writing to the Seq object's .data propery is deprecated.",
 102                       DeprecationWarning)
 103         self._data = value
 104     data = property(fget= lambda self : str(self),
 105                     fset=_set_data,
 106                     doc="Sequence as a string (DEPRECATED)")
 107
 108     def __repr__(self):
 109         """Returns a (truncated) representation of the sequence for debugging."""
 110         if len(self) > 60 :
 111             #Shows the last three letters as it is often useful to see if there
 112             #is a stop codon at the end of a sequence.
 113             #Note total length is 54+3+3=60
 114             return "%s('%s...%s', %s)" % (self.__class__.__name__,
 115                                    str(self)[:54], str(self)[-3:],
 116                                    repr(self.alphabet))
 117         else :
 118             return "%s(%s, %s)" % (self.__class__.__name__,
 119                                   repr(self.data),
 120                                    repr(self.alphabet))
 121     def __str__(self):
 122         """Returns the full sequence as a python string.
 123
 124         Note that Biopython 1.44 and earlier would give a truncated
 125         version of repr(my_seq) for str(my_seq).  If you are writing code
 126         which need to be backwards compatible with old Biopython, you
 127         should continue to use my_seq.tostring() rather than str(my_seq).
 128         """
 129         return self._data
 130
 131     """
 132     TODO - Work out why this breaks test_Restriction.py
 133     (Comparing Seq objects would be nice to have.  May need to think about
 134     hashes and the in operator for when have list/dictionary of Seq objects...)
 135     def __cmp__(self, other):
 136         if hasattr(other, "alphabet") :
 137             #other should be a Seq or a MutableSeq
 138             if not Alphabet._check_type_compatible([self.alphabet,
 139                                                     other.alphabet]) :
 140                 raise TypeError("Incompatable alphabets %s and %s" \
 141                                 % (repr(self.alphabet), repr(other.alphabet)))
 142             #They should be the same sequence type (or one of them is generic)
 143             return cmp(str(self), str(other))
 144         elif isinstance(other, basestring) :
 145             return cmp(str(self), other)
 146         else :
 147             raise TypeError
 148     """
 149
 150     def __len__(self): return len(self._data)       # Seq API requirement
 151
 152     def __getitem__(self, index) :                 # Seq API requirement
 153         #Note since Python 2.0, __getslice__ is deprecated
 154         #and __getitem__ is used instead.
 155         #See http://docs.python.org/ref/sequence-methods.html
 156         if isinstance(index, int) :
 157             #Return a single letter as a string
 158             return self._data[index]
 159         else :
 160             #Return the (sub)sequence as another Seq object
 161             return Seq(self._data[index], self.alphabet)
 162
 163     def __add__(self, other):
 164         """Add another sequence or string to this sequence."""
 165         if hasattr(other, "alphabet") :
 166             #other should be a Seq or a MutableSeq
 167             if not Alphabet._check_type_compatible([self.alphabet,
 168                                                     other.alphabet]) :
 169                 raise TypeError("Incompatable alphabets %s and %s" \
 170                                 % (repr(self.alphabet), repr(other.alphabet)))
 171             #They should be the same sequence type (or one of them is generic)
 172             a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
 173             return self.__class__(str(self) + str(other), a)
 174         elif isinstance(other, basestring) :
 175             #other is a plain string - use the current alphabet
 176             return self.__class__(str(self) + other, self.alphabet)
 177         else :
 178             raise TypeError
 179
 180     def __radd__(self, other):
 181         if hasattr(other, "alphabet") :
 182             #other should be a Seq or a MutableSeq
 183             if not Alphabet._check_type_compatible([self.alphabet,
 184                                                     other.alphabet]) :
 185                 raise TypeError("Incompatable alphabets %s and %s" \
 186                                 % (repr(self.alphabet), repr(other.alphabet)))
 187             #They should be the same sequence type (or one of them is generic)
 188             a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
 189             return self.__class__(str(other) + str(self), a)
 190         elif isinstance(other, basestring) :
 191             #other is a plain string - use the current alphabet
 192             return self.__class__(other + str(self), self.alphabet)
 193         else :
 194             raise TypeError
 195
 196     def tostring(self):                            # Seq API requirement
 197         """Returns the full sequence as a python string.
 198
 199         Although not formally deprecated, you are now encouraged to use
 200         str(my_seq) instead of my_seq.tostring()."""
 201         return str(self)
 202
 203     def tomutable(self):   # Needed?  Or use a function?
 204         """Returns the full sequence as a MutableSeq object.
 205
 206         >>> from Bio.Seq import Seq
 207         >>> from Bio.Alphabet import IUPAC
 208         >>> my_seq = Seq("MKQHKAMIVALIVICITAVVAAL",
 209         ...              IUPAC.protein)
 210         >>> my_seq
 211         Seq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
 212         >>> my_seq.tomutable()
 213         MutableSeq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
 214
 215         Note that the alphabet is preserved.
 216         """
 217         return MutableSeq(str(self), self.alphabet)
 218
 219     def _get_seq_str_and_check_alphabet(self, other_sequence) :
 220         """string/Seq/MutableSeq to string, checking alphabet (PRIVATE).
 221
 222         For a string argument, returns the string.
 223
 224         For a Seq or MutableSeq, it checks the alphabet is compatible
 225         (raising an exception if it isn't), and then returns a string.
 226         """
 227         try :
 228             other_alpha = other_sequence.alphabet
 229         except AttributeError :
 230             #Assume other_sequence is a string
 231             return other_sequence
 232
 233         #Other should be a Seq or a MutableSeq
 234         if not Alphabet._check_type_compatible([self.alphabet, other_alpha]) :
 235             raise TypeError("Incompatable alphabets %s and %s" \
 236                             % (repr(self.alphabet), repr(other_alpha)))
 237         #Return as a string
 238         return str(other_sequence)
 239
 240     def count(self, sub, start=0, end=sys.maxint):
 241         """Non-overlapping count method, like that of a python string.
 242
 243         This behaves like the python string method of the same name,
 244         which does a non-overlapping count!
 245
 246         Returns an integer, the number of occurrences of substring
 247         argument sub in the (sub)sequence given by [start:end].
 248         Optional arguments start and end are interpreted as in slice
 249         notation.
 250
 251         Arguments:
 252          - sub - a string or another Seq object to look for
 253          - start - optional integer, slice start
 254          - end - optional integer, slice end
 255
 256         e.g.
 257
 258         >>> from Bio.Seq import Seq
 259         >>> my_seq = Seq("AAAATGA")
 260         >>> print my_seq.count("A")
 261         5
 262         >>> print my_seq.count("ATG")
 263         1
 264         >>> print my_seq.count(Seq("AT"))
 265         1
 266         >>> print my_seq.count("AT", 2, -1)
 267         1
 268
 269         HOWEVER, please note because python strings and Seq objects (and
 270         MutableSeq objects) do a non-overlapping search, this may not give
 271         the answer you expect:
 272
 273         >>> "AAAA".count("AA")
 274         2
 275         >>> print Seq("AAAA").count("AA")
 276         2
 277
 278         A non-overlapping search would give the answer as three!
 279         """
 280         #If it has one, check the alphabet:
 281         sub_str = self._get_seq_str_and_check_alphabet(sub)
 282         return str(self).count(sub_str, start, end)
 283
 284     def find(self, sub, start=0, end=sys.maxint):
 285         """Find method, like that of a python string.
 286
 287         This behaves like the python string method of the same name.
 288
 289         Returns an integer, the index of the first occurrence of substring
 290         argument sub in the (sub)sequence given by [start:end].
 291
 292         Arguments:
 293          - sub - a string or another Seq object to look for
 294          - start - optional integer, slice start
 295          - end - optional integer, slice end
 296
 297         Returns -1 if the subsequence is NOT found.
 298
 299         e.g. Locating the first typical start codon, AUG, in an RNA sequence:
 300
 301         >>> from Bio.Seq import Seq
 302         >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
 303         >>> my_rna.find("AUG")
 304         3
 305         """
 306         #If it has one, check the alphabet:
 307         sub_str = self._get_seq_str_and_check_alphabet(sub)
 308         return str(self).find(sub_str, start, end)
 309
 310     def rfind(self, sub, start=0, end=sys.maxint):
 311         """Find from right method, like that of a python string.
 312
 313         This behaves like the python string method of the same name.
 314
 315         Returns an integer, the index of the last (right most) occurrence of
 316         substring argument sub in the (sub)sequence given by [start:end].
 317
 318         Arguments:
 319          - sub - a string or another Seq object to look for
 320          - start - optional integer, slice start
 321          - end - optional integer, slice end
 322
 323         Returns -1 if the subsequence is NOT found.
 324
 325         e.g. Locating the last typical start codon, AUG, in an RNA sequence:
 326
 327         >>> from Bio.Seq import Seq
 328         >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
 329         >>> my_rna.rfind("AUG")
 330         15
 331         """
 332         #If it has one, check the alphabet:
 333         sub_str = self._get_seq_str_and_check_alphabet(sub)
 334         return str(self).rfind(sub_str, start, end)
 335
 336     def startswith(self, prefix, start=0, end=sys.maxint) :
 337         """Does the Seq start with the given prefix?  Returns True/False.
 338
 339         This behaves like the python string method of the same name.
 340
 341         Return True if the sequence starts with the specified prefix
 342         (a string or another Seq object), False otherwise.
 343         With optional start, test sequence beginning at that position.
 344         With optional end, stop comparing sequence at that position.
 345         prefix can also be a tuple of strings to try.  e.g.
 346
 347         >>> from Bio.Seq import Seq
 348         >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
 349         >>> my_rna.startswith("GUC")
 350         True
 351         >>> my_rna.startswith("AUG")
 352         False
 353         >>> my_rna.startswith("AUG", 3)
 354         True
 355         >>> my_rna.startswith(("UCC","UCA","UCG"),1)
 356         True
 357         """
 358         #If it has one, check the alphabet:
 359         if isinstance(prefix, tuple) :
 360             #TODO - Once we drop support for Python 2.4, instead of this
 361             #loop offload to the string method (requires Python 2.5+).
 362             #Check all the alphabets first...
 363             prefix_strings = [self._get_seq_str_and_check_alphabet(p) \
 364                               for p in prefix]
 365             for prefix_str in prefix_strings :
 366                 if str(self).startswith(prefix_str, start, end) :
 367                     return True
 368             return False
 369         else :
 370             prefix_str = self._get_seq_str_and_check_alphabet(prefix)
 371             return str(self).startswith(prefix_str, start, end)
 372
 373     def endswith(self, suffix, start=0, end=sys.maxint) :
 374         """Does the Seq end with the given suffix?  Returns True/False.
 375
 376         This behaves like the python string method of the same name.
 377
 378         Return True if the sequence ends with the specified suffix
 379         (a string or another Seq object), False otherwise.
 380         With optional start, test sequence beginning at that position.
 381         With optional end, stop comparing sequence at that position.
 382         suffix can also be a tuple of strings to try.  e.g.
 383
 384         >>> from Bio.Seq import Seq
 385         >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
 386         >>> my_rna.endswith("UUG")
 387         True
 388         >>> my_rna.endswith("AUG")
 389         False
 390         >>> my_rna.endswith("AUG", 0, 18)
 391         True
 392         >>> my_rna.endswith(("UCC","UCA","UUG"))
 393         True
 394         """
 395         #If it has one, check the alphabet:
 396         if isinstance(suffix, tuple) :
 397             #TODO - Once we drop support for Python 2.4, instead of this
 398             #loop offload to the string method (requires Python 2.5+).
 399             #Check all the alphabets first...
 400             suffix_strings = [self._get_seq_str_and_check_alphabet(p) \
 401                               for p in suffix]
 402             for suffix_str in suffix_strings :
 403                 if str(self).endswith(suffix_str, start, end) :
 404                     return True
 405             return False
 406         else :
 407             suffix_str = self._get_seq_str_and_check_alphabet(suffix)
 408             return str(self).endswith(suffix_str, start, end)
 409
 410
 411     def split(self, sep=None, maxsplit=-1) :
 412         """Split method, like that of a python string.
 413
 414         This behaves like the python string method of the same name.
 415
 416         Return a list of the 'words' in the string (as Seq objects),
 417         using sep as the delimiter string.  If maxsplit is given, at
 418         most maxsplit splits are done.  If maxsplit is ommited, all
 419         splits are made.
 420
 421         Following the python string method, sep will by default be any
 422         white space (tabs, spaces, newlines) but this is unlikely to
 423         apply to biological sequences.
 424
 425         e.g.
 426
 427         >>> from Bio.Seq import Seq
 428         >>> my_rna = Seq("GUCAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUUG")
 429         >>> my_aa = my_rna.translate()
 430         >>> my_aa
 431         Seq('VMAIVMGR*KGAR*L', HasStopCodon(ExtendedIUPACProtein(), '*'))
 432         >>> my_aa.split("*")
 433         [Seq('VMAIVMGR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('KGAR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('L', HasStopCodon(ExtendedIUPACProtein(), '*'))]
 434         >>> my_aa.split("*",1)
 435         [Seq('VMAIVMGR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('KGAR*L', HasStopCodon(ExtendedIUPACProtein(), '*'))]
 436
 437         See also the rsplit method:
 438
 439         >>> my_aa.rsplit("*",1)
 440         [Seq('VMAIVMGR*KGAR', HasStopCodon(ExtendedIUPACProtein(), '*')), Seq('L', HasStopCodon(ExtendedIUPACProtein(), '*'))]
 441         """
 442         #If it has one, check the alphabet:
 443         sep_str = self._get_seq_str_and_check_alphabet(sep)
 444         #TODO - If the sep is the defined stop symbol, or gap char,
 445         #should we adjust the alphabet?
 446         return [Seq(part, self.alphabet) \
 447                 for part in str(self).split(sep_str, maxsplit)]
 448
 449     def rsplit(self, sep=None, maxsplit=-1) :
 450         """Right split method, like that of a python string.
 451
 452         This behaves like the python string method of the same name.
 453
 454         Return a list of the 'words' in the string (as Seq objects),
 455         using sep as the delimiter string.  If maxsplit is given, at
 456         most maxsplit splits are done COUNTING FROM THE RIGHT.
 457         If maxsplit is ommited, all splits are made.
 458
 459         Following the python string method, sep will by default be any
 460         white space (tabs, spaces, newlines) but this is unlikely to
 461         apply to biological sequences.
 462
 463         e.g. print my_seq.rsplit("*",1)
 464
 465         See also the split method.
 466         """
 467         #If it has one, check the alphabet:
 468         sep_str = self._get_seq_str_and_check_alphabet(sep)
 469         try :
 470             return [Seq(part, self.alphabet) \
 471                     for part in str(self).rsplit(sep_str, maxsplit)]
 472         except AttributeError :
 473             #Python 2.3 doesn't have a string rsplit method, which we can
 474             #word around by reversing the sequence, using (left) split,
 475             #and then reversing the answer. Not very efficient!
 476             words = [Seq(word[::-1], self.alphabet) for word \
 477                      in str(self)[::-1].split(sep_str[::-1], maxsplit)]
 478             words.reverse()
 479             return words
 480
 481     def strip(self, chars=None) :
 482         """Returns a new Seq object with leading and trailing ends stripped.
 483
 484         This behaves like the python string method of the same name.
 485
 486         Optional argument chars defines which characters to remove.  If
 487         ommitted or None (default) then as for the python string method,
 488         this defaults to removing any white space.
 489
 490         e.g. print my_seq.strip("-")
 491
 492         See also the lstrip and rstrip methods.
 493         """
 494         #If it has one, check the alphabet:
 495         strip_str = self._get_seq_str_and_check_alphabet(chars)
 496         return Seq(str(self).strip(strip_str), self.alphabet)
 497
 498     def lstrip(self, chars=None) :
 499         """Returns a new Seq object with leading (left) end stripped.
 500
 501         This behaves like the python string method of the same name.
 502
 503         Optional argument chars defines which characters to remove.  If
 504         ommitted or None (default) then as for the python string method,
 505         this defaults to removing any white space.
 506
 507         e.g. print my_seq.lstrip("-")
 508
 509         See also the strip and rstrip methods.
 510         """
 511         #If it has one, check the alphabet:
 512         strip_str = self._get_seq_str_and_check_alphabet(chars)
 513         return Seq(str(self).lstrip(strip_str), self.alphabet)
 514
 515     def rstrip(self, chars=None) :
 516         """Returns a new Seq object with trailing (right) end stripped.
 517
 518         This behaves like the python string method of the same name.
 519
 520         Optional argument chars defines which characters to remove.  If
 521         ommitted or None (default) then as for the python string method,
 522         this defaults to removing any white space.
 523
 524         e.g. Removing a nucleotide sequence's polyadenylation (poly-A tail):
 525
 526         >>> from Bio.Alphabet import IUPAC
 527         >>> from Bio.Seq import Seq
 528         >>> my_seq = Seq("CGGTACGCTTATGTCACGTAGAAAAAA", IUPAC.unambiguous_dna)
 529         >>> my_seq
 530         Seq('CGGTACGCTTATGTCACGTAGAAAAAA', IUPACUnambiguousDNA())
 531         >>> my_seq.rstrip("A")
 532         Seq('CGGTACGCTTATGTCACGTAG', IUPACUnambiguousDNA())
 533
 534         See also the strip and lstrip methods.
 535         """
 536         #If it has one, check the alphabet:
 537         strip_str = self._get_seq_str_and_check_alphabet(chars)
 538         return Seq(str(self).rstrip(strip_str), self.alphabet)
 539
 540     def complement(self):
 541         """Returns the complement sequence. New Seq object.
 542
 543         >>> from Bio.Seq import Seq
 544         >>> from Bio.Alphabet import IUPAC
 545         >>> my_dna = Seq("CCCCCGATAG", IUPAC.unambiguous_dna)
 546         >>> my_dna
 547         Seq('CCCCCGATAG', IUPACUnambiguousDNA())
 548         >>> my_dna.complement()
 549         Seq('GGGGGCTATC', IUPACUnambiguousDNA())
 550
 551         You can of course used mixed case sequences,
 552
 553         >>> from Bio.Seq import Seq
 554         >>> from Bio.Alphabet import generic_dna
 555         >>> my_dna = Seq("CCCCCgatA-GD", generic_dna)
 556         >>> my_dna
 557         Seq('CCCCCgatA-GD', DNAAlphabet())
 558         >>> my_dna.complement()
 559         Seq('GGGGGctaT-CH', DNAAlphabet())
 560
 561         Note in the above example, ambiguous character D denotes
 562         G, A or T so its complement is H (for C, T or A).
 563
 564         Trying to complement a protein sequence raises an exception.
 565
 566         >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
 567         >>> my_protein.complement()
 568         Traceback (most recent call last):
 569            ...
 570         ValueError: Proteins do not have complements!
 571         """
 572         base = Alphabet._get_base_alphabet(self.alphabet)
 573         if isinstance(base, Alphabet.ProteinAlphabet) :
 574             raise ValueError("Proteins do not have complements!")
 575         if isinstance(base, Alphabet.DNAAlphabet) :
 576             ttable = _dna_complement_table
 577         elif isinstance(base, Alphabet.RNAAlphabet) :
 578             ttable = _rna_complement_table
 579         elif ('U' in self._data or 'u' in self._data) \
 580         and ('T' in self._data or 't' in self._data):
 581             #TODO - Handle this cleanly?
 582             raise ValueError("Mixed RNA/DNA found")
 583         elif 'U' in self._data or 'u' in self._data:
 584             ttable = _rna_complement_table
 585         else:
 586             ttable = _dna_complement_table
 587         #Much faster on really long sequences than the previous loop based one.
 588         #thx to Michael Palmer, University of Waterloo
 589         return Seq(str(self).translate(ttable), self.alphabet)
 590
 591     def reverse_complement(self):
 592         """Returns the reverse complement sequence. New Seq object.
 593
 594         >>> from Bio.Seq import Seq
 595         >>> from Bio.Alphabet import IUPAC
 596         >>> my_dna = Seq("CCCCCGATAGNR", IUPAC.ambiguous_dna)
 597         >>> my_dna
 598         Seq('CCCCCGATAGNR', IUPACAmbiguousDNA())
 599         >>> my_dna.reverse_complement()
 600         Seq('YNCTATCGGGGG', IUPACAmbiguousDNA())
 601
 602         Note in the above example, since R = G or A, its complement
 603         is Y (which denotes  C or T).
 604
 605         You can of course used mixed case sequences,
 606
 607         >>> from Bio.Seq import Seq
 608         >>> from Bio.Alphabet import generic_dna
 609         >>> my_dna = Seq("CCCCCgatA-G", generic_dna)
 610         >>> my_dna
 611         Seq('CCCCCgatA-G', DNAAlphabet())
 612         >>> my_dna.reverse_complement()
 613         Seq('C-TatcGGGGG', DNAAlphabet())
 614
 615         Trying to complement a protein sequence raises an exception:
 616
 617         >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
 618         >>> my_protein.reverse_complement()
 619         Traceback (most recent call last):
 620            ...
 621         ValueError: Proteins do not have complements!
 622         """
 623         #Use -1 stride/step to reverse the complement
 624         return self.complement()[::-1]
 625
 626     def transcribe(self):
 627         """Returns the RNA sequence from a DNA sequence. New Seq object.
 628
 629         >>> from Bio.Seq import Seq
 630         >>> from Bio.Alphabet import IUPAC
 631         >>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG",
 632         ...                  IUPAC.unambiguous_dna)
 633         >>> coding_dna
 634         Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())
 635         >>> coding_dna.transcribe()
 636         Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())
 637
 638         Trying to transcribe a protein or RNA sequence raises an exception:
 639
 640         >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
 641         >>> my_protein.transcribe()
 642         Traceback (most recent call last):
 643            ...
 644         ValueError: Proteins cannot be transcribed!
 645         """
 646         base = Alphabet._get_base_alphabet(self.alphabet)
 647         if isinstance(base, Alphabet.ProteinAlphabet) :
 648             raise ValueError("Proteins cannot be transcribed!")
 649         if isinstance(base, Alphabet.RNAAlphabet) :
 650             raise ValueError("RNA cannot be transcribed!")
 651
 652         if self.alphabet==IUPAC.unambiguous_dna:
 653             alphabet = IUPAC.unambiguous_rna
 654         elif self.alphabet==IUPAC.ambiguous_dna:
 655             alphabet = IUPAC.ambiguous_rna
 656         else:
 657             alphabet = Alphabet.generic_rna
 658         return Seq(str(self).replace('T','U').replace('t','u'), alphabet)
 659
 660     def back_transcribe(self):
 661         """Returns the DNA sequence from an RNA sequence. New Seq object.
 662
 663         >>> from Bio.Seq import Seq
 664         >>> from Bio.Alphabet import IUPAC
 665         >>> messenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG",
 666         ...                     IUPAC.unambiguous_rna)
 667         >>> messenger_rna
 668         Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())
 669         >>> messenger_rna.back_transcribe()
 670         Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())
 671
 672         Trying to back-transcribe a protein or DNA sequence raises an
 673         exception:
 674
 675         >>> my_protein = Seq("MAIVMGR", IUPAC.protein)
 676         >>> my_protein.back_transcribe()
 677         Traceback (most recent call last):
 678            ...
 679         ValueError: Proteins cannot be back transcribed!
 680         """
 681         base = Alphabet._get_base_alphabet(self.alphabet)
 682         if isinstance(base, Alphabet.ProteinAlphabet) :
 683             raise ValueError("Proteins cannot be back transcribed!")
 684         if isinstance(base, Alphabet.DNAAlphabet) :
 685             raise ValueError("DNA cannot be back transcribed!")
 686
 687         if self.alphabet==IUPAC.unambiguous_rna:
 688             alphabet = IUPAC.unambiguous_dna
 689         elif self.alphabet==IUPAC.ambiguous_rna:
 690             alphabet = IUPAC.ambiguous_dna
 691         else:
 692             alphabet = Alphabet.generic_dna
 693         return Seq(str(self).replace("U", "T").replace("u", "t"), alphabet)
 694
 695     def translate(self, table="Standard", stop_symbol="*", to_stop=False):
 696         """Turns a nucleotide sequence into a protein sequence. New Seq object.
 697
 698         This method will translate DNA or RNA sequences, and those with a
 699         nucleotide or generic alphabet.  Trying to translate a protein
 700         sequence raises an exception.
 701
 702         Arguments:
 703          - table - Which codon table to use?  This can be either a name
 704                    (string) or an NCBI identifier (integer).  This defaults
 705                    to the "Standard" table.
 706          - stop_symbol - Single character string, what to use for terminators.
 707                          This defaults to the asterisk, "*".
 708          - to_stop - Boolean, defaults to False meaning do a full translation
 709                      continuing on past any stop codons (translated as the
 710                      specified stop_symbol).  If True, translation is
 711                      terminated at the first in frame stop codon (and the
 712                      stop_symbol is not appended to the returned protein
 713                      sequence).
 714
 715         e.g. Using the standard table:
 716
 717         >>> coding_dna = Seq("GTGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
 718         >>> coding_dna.translate()
 719         Seq('VAIVMGR*KGAR*', HasStopCodon(ExtendedIUPACProtein(), '*'))
 720         >>> coding_dna.translate(stop_symbol="@")
 721         Seq('VAIVMGR@KGAR@', HasStopCodon(ExtendedIUPACProtein(), '@'))
 722         >>> coding_dna.translate(to_stop=True)
 723         Seq('VAIVMGR', ExtendedIUPACProtein())
 724
 725         Now using NCBI table 2, where TGA is not a stop codon:
 726
 727         >>> coding_dna.translate(table=2)
 728         Seq('VAIVMGRWKGAR*', HasStopCodon(ExtendedIUPACProtein(), '*'))
 729         >>> coding_dna.translate(table=2, to_stop=True)
 730         Seq('VAIVMGRWKGAR', ExtendedIUPACProtein())
 731
 732         If the sequence has no in-frame stop codon, then the to_stop argument
 733         has no effect:
 734
 735         >>> coding_dna2 = Seq("TTGGCCATTGTAATGGGCCGC")
 736         >>> coding_dna2.translate()
 737         Seq('LAIVMGR', ExtendedIUPACProtein())
 738         >>> coding_dna2.translate(to_stop=True)
 739         Seq('LAIVMGR', ExtendedIUPACProtein())
 740
 741         NOTE - Ambiguous codons like "TAN" or "NNN" could be an amino acid
 742         or a stop codon.  These are translated as "X".  Any invalid codon
 743         (e.g. "TA?" or "T-A") will throw a TranslationError.
 744
 745         NOTE - Does NOT support gapped sequences.
 746
 747         NOTE - This does NOT behave like the python string's translate
 748         method.  For that use str(my_seq).translate(...) instead.
 749         """
 750         try:
 751             table_id = int(table)
 752         except ValueError:
 753             table_id = None
 754         if isinstance(table, str) and len(table)==256 :
 755             raise ValueError("The Seq object translate method DOES NOT take " \
 756                              + "a 256 character string mapping table like " \
 757                              + "the python string object's translate method. " \
 758                              + "Use str(my_seq).translate(...) instead.")
 759         if isinstance(Alphabet._get_base_alphabet(self.alphabet),
 760                       Alphabet.ProteinAlphabet) :
 761             raise ValueError("Proteins cannot be translated!")
 762         if self.alphabet==IUPAC.unambiguous_dna:
 763             #Will use standard IUPAC protein alphabet, no need for X
 764             if table_id is None:
 765                 codon_table = CodonTable.unambiguous_dna_by_name[table]
 766             else:
 767                 codon_table = CodonTable.unambiguous_dna_by_id[table_id]
 768         #elif self.alphabet==IUPAC.ambiguous_dna:
 769         #    if table_id is None:
 770         #        codon_table = CodonTable.ambiguous_dna_by_name[table]
 771         #    else:
 772         #        codon_table = CodonTable.ambiguous_dna_by_id[table_id]
 773         elif self.alphabet==IUPAC.unambiguous_rna:
 774             #Will use standard IUPAC protein alphabet, no need for X
 775             if table_id is None:
 776                 codon_table = CodonTable.unambiguous_rna_by_name[table]
 777             else:
 778                 codon_table = CodonTable.unambiguous_rna_by_id[table_id]
 779         #elif self.alphabet==IUPAC.ambiguous_rna:
 780         #    if table_id is None:
 781         #        codon_table = CodonTable.ambiguous_rna_by_name[table]
 782         #    else:
 783         #        codon_table = CodonTable.ambiguous_rna_by_id[table_id]
 784         else:
 785             #This will use the extend IUPAC protein alphabet with X etc.
 786             #The same table can be used for RNA or DNA (we use this for
 787             #translating strings).
 788             if table_id is None:
 789                 codon_table = CodonTable.ambiguous_generic_by_name[table]
 790             else:
 791                 codon_table = CodonTable.ambiguous_generic_by_id[table_id]
 792         protein = _translate_str(str(self), codon_table, stop_symbol, to_stop)
 793         if stop_symbol in protein :
 794             alphabet = Alphabet.HasStopCodon(codon_table.protein_alphabet,
 795                                              stop_symbol = stop_symbol)
 796         else :
 797             alphabet = codon_table.protein_alphabet
 798         return Seq(protein, alphabet)
 799
 800 class UnknownSeq(Seq):
 801     """A read-only sequence object of known length but unknown contents.
 802
 803     If you have an unknown sequence, you can represent this with a normal
 804     Seq object, for example:
 805
 806     >>> my_seq = Seq("N"*5)
 807     >>> my_seq
 808     Seq('NNNNN', Alphabet())
 809     >>> len(my_seq)
 810     5
 811     >>> print my_seq
 812     NNNNN
 813
 814     However, this is rather wasteful of memory (especially for large
 815     sequences), which is where this class is most usefull:
 816
 817     >>> unk_five = UnknownSeq(5)
 818     >>> unk_five
 819     UnknownSeq(5, alphabet = Alphabet(), character = '?')
 820     >>> len(unk_five)
 821     5
 822     >>> print(unk_five)
 823     ?????
 824
 825     You can add unknown sequence together, provided their alphabets and
 826     characters are compatible, and get another memory saving UnknownSeq:
 827
 828     >>> unk_four = UnknownSeq(4)
 829     >>> unk_four
 830     UnknownSeq(4, alphabet = Alphabet(), character = '?')
 831     >>> unk_four + unk_five
 832     UnknownSeq(9, alphabet = Alphabet(), character = '?')
 833
 834     If the alphabet or characters don't match up, the addition gives an
 835     ordinary Seq object:
 836
 837     >>> unk_nnnn = UnknownSeq(4, character = "N")
 838     >>> unk_nnnn
 839     UnknownSeq(4, alphabet = Alphabet(), character = 'N')
 840     >>> unk_nnnn + unk_four
 841     Seq('NNNN????', Alphabet())
 842
 843     Combining with a real Seq gives a new Seq object:
 844
 845     >>> known_seq = Seq("ACGT")
 846     >>> unk_four + known_seq
 847     Seq('????ACGT', Alphabet())
 848     >>> known_seq + unk_four
 849     Seq('ACGT????', Alphabet())
 850     """
 851     def __init__(self, length, alphabet = Alphabet.generic_alphabet, character = None) :
 852         """Create a new UnknownSeq object.
 853
 854         If character is ommited, it is determed from the alphabet, "N" for
 855         nucleotides, "X" for proteins, and "?" otherwise.
 856         """
 857         self._length = int(length)
 858         if self._length < 0 :
 859             #TODO - Block zero length UnknownSeq?  You can just use a Seq!
 860             raise ValueError("Length must not be negative.")
 861         self.alphabet = alphabet
 862         if character :
 863             if len(character) != 1 :
 864                 raise ValueError("character argument should be a single letter string.")
 865             self._character = character
 866         else :
 867             base = Alphabet._get_base_alphabet(alphabet)
 868             #TODO? Check the case of the letters in the alphabet?
 869             #We may have to use "n" instead of "N" etc.
 870             if isinstance(base, Alphabet.NucleotideAlphabet) :
 871                 self._character = "N"
 872             elif isinstance(base, Alphabet.ProteinAlphabet) :
 873                 self._character = "X"
 874             else :
 875                 self._character = "?"
 876
 877     def __len__(self) :
 878         """Returns the stated length of the unknown sequence."""
 879         return self._length
 880
 881     def __str__(self) :
 882         """Returns the unknown sequence as full string of the given length."""
 883         return self._character * self._length
 884
 885     def __repr__(self):
 886         return "UnknownSeq(%i, alphabet = %s, character = %s)" \
 887                % (self._length, repr(self.alphabet), repr(self._character))
 888
 889     def __add__(self, other) :
 890         if isinstance(other, UnknownSeq) \
 891         and other._character == self._character :
 892             #TODO - Check the alphabets match
 893             return UnknownSeq(len(self)+len(other),
 894                               self.alphabet, self._character)
 895         #Offload to the base class...
 896         return Seq(str(self), self.alphabet) + other
 897
 898     def __radd__(self, other) :
 899         if isinstance(other, UnknownSeq) \
 900         and other._character == self._character :
 901             #TODO - Check the alphabets match
 902             return UnknownSeq(len(self)+len(other),
 903                               self.alphabet, self._character)
 904         #Offload to the base class...
 905         return other + Seq(str(self), self.alphabet)
 906
 907     def __getitem__(self, index):
 908         if isinstance(index, int) :
 909             #TODO - Check the bounds without wasting memory
 910             return str(self)[index]
 911         else :
 912             #TODO - Work out the length without wasting memory
 913             return UnknownSeq(len(("#"*self._length)[index]),
 914                               self.alphabet, self._character)
 915
 916     def count(self, sub, start=0, end=sys.maxint):
 917         """Non-overlapping count method, like that of a python string.
 918
 919         This behaves like the python string (and Seq object) method of the
 920         same name, which does a non-overlapping count!
 921
 922         Returns an integer, the number of occurrences of substring
 923         argument sub in the (sub)sequence given by [start:end].
 924         Optional arguments start and end are interpreted as in slice
 925         notation.
 926
 927         Arguments:
 928          - sub - a string or another Seq object to look for
 929          - start - optional integer, slice start
 930          - end - optional integer, slice end
 931
 932         >>> "NNNN".count("N")
 933         4
 934         >>> Seq("NNNN").count("N")
 935         4
 936         >>> UnknownSeq(4, character="N").count("N")
 937         4
 938         >>> UnknownSeq(4, character="N").count("A")
 939         0
 940         >>> UnknownSeq(4, character="N").count("AA")
 941         0
 942
 943         HOWEVER, please note because that python strings and Seq objects (and
 944         MutableSeq objects) do a non-overlapping search, this may not give
 945         the answer you expect:
 946
 947         >>> UnknownSeq(4, character="N").count("NN")
 948         2
 949         >>> UnknownSeq(4, character="N").count("NNN")
 950         1
 951         """
 952         sub_str = self._get_seq_str_and_check_alphabet(sub)
 953         if len(sub_str) == 1 :
 954             if str(sub_str) == self._character :
 955                 if start==0 and end >= self._length :
 956                     return self._length
 957                 else :
 958                     #This could be done more cleverly...
 959                     return str(self).count(sub_str, start, end)
 960             else :
 961                 return 0
 962         else :
 963             if set(sub_str) == set(self._character) :
 964                 if start==0 and end >= self._length :
 965                     return self._length // len(sub_str)
 966                 else :
 967                     #This could be done more cleverly...
 968                     return str(self).count(sub_str, start, end)
 969             else :
 970                 return 0
 971
 972     def complement(self) :
 973         """The complement of an unknown nucleotide equals itself.
 974
 975         >>> my_nuc = UnknownSeq(8)
 976         >>> my_nuc
 977         UnknownSeq(8, alphabet = Alphabet(), character = '?')
 978         >>> print my_nuc
 979         ????????
 980         >>> my_nuc.complement()
 981         UnknownSeq(8, alphabet = Alphabet(), character = '?')
 982         >>> print my_nuc.complement()
 983         ????????
 984         """
 985         if isinstance(Alphabet._get_base_alphabet(self.alphabet),
 986                       Alphabet.ProteinAlphabet) :
 987             raise ValueError("Proteins do not have complements!")
 988         return self
 989
 990     def reverse_complement(self) :
 991         """The reverse complement of an unknown nucleotide equals itself.
 992
 993         >>> my_nuc = UnknownSeq(10)
 994         >>> my_nuc
 995         UnknownSeq(10, alphabet = Alphabet(), character = '?')
 996         >>> print my_nuc
 997         ??????????
 998         >>> my_nuc.reverse_complement()
 999         UnknownSeq(10, alphabet = Alphabet(), character = '?')
1000         >>> print my_nuc.reverse_complement()
1001         ??????????
1002         """
1003         if isinstance(Alphabet._get_base_alphabet(self.alphabet),
1004                       Alphabet.ProteinAlphabet) :
1005             raise ValueError("Proteins do not have complements!")
1006         return self
1007
1008     def transcribe(self) :
1009         """Returns unknown RNA sequence from an unknown DNA sequence.
1010
1011         >>> my_dna = UnknownSeq(10, character="N")
1012         >>> my_dna
1013         UnknownSeq(10, alphabet = Alphabet(), character = 'N')
1014         >>> print my_dna
1015         NNNNNNNNNN
1016         >>> my_rna = my_dna.transcribe()
1017         >>> my_rna
1018         UnknownSeq(10, alphabet = RNAAlphabet(), character = 'N')
1019         >>> print my_rna
1020         NNNNNNNNNN
1021         """
1022         #Offload the alphabet stuff
1023         s = Seq(self._character, self.alphabet).transcribe()
1024         return UnknownSeq(self._length, s.alphabet, self._character)
1025
1026     def back_transcribe(self) :
1027         """Returns unknown DNA sequence from an unknown RNA sequence.
1028
1029         >>> my_rna = UnknownSeq(20, character="N")
1030         >>> my_rna
1031         UnknownSeq(20, alphabet = Alphabet(), character = 'N')
1032         >>> print my_rna
1033         NNNNNNNNNNNNNNNNNNNN
1034         >>> my_dna = my_rna.back_transcribe()
1035         >>> my_dna
1036         UnknownSeq(20, alphabet = DNAAlphabet(), character = 'N')
1037         >>> print my_dna
1038         NNNNNNNNNNNNNNNNNNNN
1039         """
1040         #Offload the alphabet stuff
1041         s = Seq(self._character, self.alphabet).back_transcribe()
1042         return UnknownSeq(self._length, s.alphabet, self._character)
1043
1044     def translate(self, **kwargs) :
1045         """Translate an unknown nucleotide sequence into an unknown protein.
1046
1047         e.g.
1048
1049         >>> my_seq = UnknownSeq(11, character="N")
1050         >>> print my_seq
1051         NNNNNNNNNNN
1052         >>> my_protein = my_seq.translate()
1053         >>> my_protein
1054         UnknownSeq(3, alphabet = ProteinAlphabet(), character = 'X')
1055         >>> print my_protein
1056         XXX
1057
1058         In comparison, using a normal Seq object:
1059
1060         >>> my_seq = Seq("NNNNNNNNNNN")
1061         >>> print my_seq
1062         NNNNNNNNNNN
1063         >>> my_protein = my_seq.translate()
1064         >>> my_protein
1065         Seq('XXX', ExtendedIUPACProtein())
1066         >>> print my_protein
1067         XXX
1068
1069         """
1070         if isinstance(Alphabet._get_base_alphabet(self.alphabet),
1071                       Alphabet.ProteinAlphabet) :
1072             raise ValueError("Proteins cannot be translated!")
1073         return UnknownSeq(self._length//3, Alphabet.generic_protein, "X")
1074
1075
1076 class MutableSeq(object):
1077     """An editable sequence object (with an alphabet).
1078
1079     Unlike normal python strings and our basic sequence object (the Seq class)
1080     which are immuatable, the MutableSeq lets you edit the sequence in place.
1081     However, this means you cannot use a MutableSeq object as a dictionary key.
1082
1083     >>> from Bio.Seq import MutableSeq
1084     >>> from Bio.Alphabet import generic_dna
1085     >>> my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
1086     >>> my_seq
1087     MutableSeq('ACTCGTCGTCG', DNAAlphabet())
1088     >>> my_seq[5]
1089     'T'
1090     >>> my_seq[5] = "A"
1091     >>> my_seq
1092     MutableSeq('ACTCGACGTCG', DNAAlphabet())
1093     >>> my_seq[5]
1094     'A'
1095     >>> my_seq[5:8] = "NNN"
1096     >>> my_seq
1097     MutableSeq('ACTCGNNNTCG', DNAAlphabet())
1098     >>> len(my_seq)
1099     11
1100
1101     Note that the MutableSeq object does not support as many string-like
1102     or biological methods as the Seq object.
1103     """
1104     def __init__(self, data, alphabet = Alphabet.generic_alphabet):
1105         if type(data) == type(""):
1106             self.data = array.array("c", data)
1107         else:
1108             self.data = data   # assumes the input is an array
1109         self.alphabet = alphabet
1110
1111     def __repr__(self):
1112         """Returns a (truncated) representation of the sequence for debugging."""
1113         if len(self) > 60 :
1114             #Shows the last three letters as it is often useful to see if there
1115             #is a stop codon at the end of a sequence.
1116             #Note total length is 54+3+3=60
1117             return "%s('%s...%s', %s)" % (self.__class__.__name__,
1118                                    str(self[:54]), str(self[-3:]),
1119                                    repr(self.alphabet))
1120         else :
1121             return "%s('%s', %s)" % (self.__class__.__name__,
1122                                    str(self),
1123                                    repr(self.alphabet))
1124
1125     def __str__(self):
1126         """Returns the full sequence as a python string.
1127
1128         Note that Biopython 1.44 and earlier would give a truncated
1129         version of repr(my_seq) for str(my_seq).  If you are writing code
1130         which needs to be backwards compatible with old Biopython, you
1131         should continue to use my_seq.tostring() rather than str(my_seq).
1132         """
1133         #See test_GAQueens.py for an historic usage of a non-string alphabet!
1134         return "".join(self.data)
1135
1136     def __cmp__(self, other):
1137         """Compare the sequence for to another sequence or a string.
1138
1139         If compared to another sequence the alphabets must be compatible.
1140         Comparing DNA to RNA, or Nucleotide to Protein will raise an
1141         exception.
1142
1143         Otherwise only the sequence itself is compared, not the precise
1144         alphabet.
1145
1146         This method indirectly supports ==, < , etc."""
1147         if hasattr(other, "alphabet") :
1148             #other should be a Seq or a MutableSeq
1149             if not Alphabet._check_type_compatible([self.alphabet,
1150                                                     other.alphabet]) :
1151                 raise TypeError("Incompatable alphabets %s and %s" \
1152                                 % (repr(self.alphabet), repr(other.alphabet)))
1153             #They should be the same sequence type (or one of them is generic)
1154             if isinstance(other, MutableSeq):
1155                 #See test_GAQueens.py for an historic usage of a non-string
1156                 #alphabet!  Comparing the arrays supports this.
1157                 return cmp(self.data, other.data)
1158             else :
1159                 return cmp(str(self), str(other))
1160         elif isinstance(other, basestring) :
1161             return cmp(str(self), other)
1162         else :
1163             raise TypeError
1164
1165     def __len__(self): return len(self.data)
1166
1167     def __getitem__(self, index) :
1168         #Note since Python 2.0, __getslice__ is deprecated
1169         #and __getitem__ is used instead.
1170         #See http://docs.python.org/ref/sequence-methods.html
1171         if isinstance(index, int) :
1172             #Return a single letter as a string
1173             return self.data[index]
1174         else :
1175             #Return the (sub)sequence as another Seq object
1176             return MutableSeq(self.data[index], self.alphabet)
1177
1178     def __setitem__(self, index, value):
1179         #Note since Python 2.0, __setslice__ is deprecated
1180         #and __setitem__ is used instead.
1181         #See http://docs.python.org/ref/sequence-methods.html
1182         if isinstance(index, int) :
1183             #Replacing a single letter with a new string
1184             self.data[index] = value
1185         else :
1186             #Replacing a sub-sequence
1187             if isinstance(value, MutableSeq):
1188                 self.data[index] = value.data
1189             elif isinstance(value, type(self.data)):
1190                 self.data[index] = value
1191             else:
1192                 self.data[index] = array.array("c", str(value))
1193
1194     def __delitem__(self, index):
1195         #Note since Python 2.0, __delslice__ is deprecated
1196         #and __delitem__ is used instead.
1197         #See http://docs.python.org/ref/sequence-methods.html
1198
1199         #Could be deleting a single letter, or a slice
1200         del self.data[index]
1201
1202     def __add__(self, other):
1203         """Add another sequence or string to this sequence.
1204
1205         Returns a new MutableSeq object."""
1206         if hasattr(other, "alphabet") :
1207             #other should be a Seq or a MutableSeq
1208             if not Alphabet._check_type_compatible([self.alphabet,
1209                                                     other.alphabet]) :
1210                 raise TypeError("Incompatable alphabets %s and %s" \
1211                                 % (repr(self.alphabet), repr(other.alphabet)))
1212             #They should be the same sequence type (or one of them is generic)
1213             a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
1214             if isinstance(other, MutableSeq):
1215                 #See test_GAQueens.py for an historic usage of a non-string
1216                 #alphabet!  Adding the arrays should support this.
1217                 return self.__class__(self.data + other.data, a)
1218             else :
1219                 return self.__class__(str(self) + str(other), a)
1220         elif isinstance(other, basestring) :
1221             #other is a plain string - use the current alphabet
1222             return self.__class__(str(self) + str(other), self.alphabet)
1223         else :
1224             raise TypeError
1225
1226     def __radd__(self, other):
1227         if hasattr(other, "alphabet") :
1228             #other should be a Seq or a MutableSeq
1229             if not Alphabet._check_type_compatible([self.alphabet,
1230                                                     other.alphabet]) :
1231                 raise TypeError("Incompatable alphabets %s and %s" \
1232                                 % (repr(self.alphabet), repr(other.alphabet)))
1233             #They should be the same sequence type (or one of them is generic)
1234             a = Alphabet._consensus_alphabet([self.alphabet, other.alphabet])
1235             if isinstance(other, MutableSeq):
1236                 #See test_GAQueens.py for an historic usage of a non-string
1237                 #alphabet!  Adding the arrays should support this.
1238                 return self.__class__(other.data + self.data, a)
1239             else :
1240                 return self.__class__(str(other) + str(self), a)
1241         elif isinstance(other, basestring) :
1242             #other is a plain string - use the current alphabet
1243             return self.__class__(str(other) + str(self), self.alphabet)
1244         else :
1245             raise TypeError
1246
1247     def append(self, c):
1248         self.data.append(c)
1249
1250     def insert(self, i, c):
1251         self.data.insert(i, c)
1252
1253     def pop(self, i = (-1)):
1254         c = self.data[i]
1255         del self.data[i]
1256         return c
1257
1258     def remove(self, item):
1259         for i in range(len(self.data)):
1260             if self.data[i] == item:
1261                 del self.data[i]
1262                 return
1263         raise ValueError("MutableSeq.remove(x): x not in list")
1264
1265     def count(self, sub, start=0, end=sys.maxint):
1266         """Non-overlapping count method, like that of a python string.
1267
1268         This behaves like the python string method of the same name,
1269         which does a non-overlapping count!
1270
1271         Returns an integer, the number of occurrences of substring
1272         argument sub in the (sub)sequence given by [start:end].
1273         Optional arguments start and end are interpreted as in slice
1274         notation.
1275
1276         Arguments:
1277          - sub - a string or another Seq object to look for
1278          - start - optional integer, slice start
1279          - end - optional integer, slice end
1280
1281         e.g.
1282
1283         >>> from Bio.Seq import MutableSeq
1284         >>> my_mseq = MutableSeq("AAAATGA")
1285         >>> print my_mseq.count("A")
1286         5
1287         >>> print my_mseq.count("ATG")
1288         1
1289         >>> print my_mseq.count(Seq("AT"))
1290         1
1291         >>> print my_mseq.count("AT", 2, -1)
1292         1
1293
1294         HOWEVER, please note because that python strings, Seq objects and
1295         MutableSeq objects do a non-overlapping search, this may not give
1296         the answer you expect:
1297
1298         >>> "AAAA".count("AA")
1299         2
1300         >>> print MutableSeq("AAAA").count("AA")
1301         2
1302
1303         A non-overlapping search would give the answer as three!
1304         """
1305         try :
1306             #TODO - Should we check the alphabet?
1307             search = sub.tostring()
1308         except AttributeError :
1309             search = sub
1310
1311         if not isinstance(search, basestring) :
1312             raise TypeError("expected a string, Seq or MutableSeq")
1313
1314         if len(search) == 1 :
1315             #Try and be efficient and work directly from the array.
1316             count = 0
1317             for c in self.data[start:end]:
1318                 if c == search: count += 1
1319             return count
1320         else :
1321             #TODO - Can we do this more efficiently?
1322             return self.tostring().count(search, start, end)
1323
1324     def index(self, item):
1325         for i in range(len(self.data)):
1326             if self.data[i] == item:
1327                 return i
1328         raise ValueError("MutableSeq.index(x): x not in list")
1329
1330     def reverse(self):
1331         """Modify the mutable sequence to reverse itself.
1332
1333         No return value.
1334         """
1335         self.data.reverse()
1336
1337     def complement(self):
1338         """Modify the mutable sequence to take on its complement.
1339
1340         Trying to complement a protein sequence raises an exception.
1341
1342         No return value.
1343         """
1344         if isinstance(Alphabet._get_base_alphabet(self.alphabet),
1345                       Alphabet.ProteinAlphabet) :
1346             raise ValueError("Proteins do not have complements!")
1347         if self.alphabet in (IUPAC.ambiguous_dna, IUPAC.unambiguous_dna):
1348             d = ambiguous_dna_complement
1349         elif self.alphabet in (IUPAC.ambiguous_rna, IUPAC.unambiguous_rna):
1350             d = ambiguous_rna_complement
1351         elif 'U' in self.data and 'T' in self.data :
1352             #TODO - Handle this cleanly?
1353             raise ValueError("Mixed RNA/DNA found")
1354         elif 'U' in self.data:
1355             d = ambiguous_rna_complement
1356         else:
1357             d = ambiguous_dna_complement
1358         c = dict([(x.lower(), y.lower()) for x,y in d.iteritems()])
1359         d.update(c)
1360         self.data = map(lambda c: d[c], self.data)
1361         self.data = array.array('c', self.data)
1362
1363     def reverse_complement(self):
1364         """Modify the mutable sequence to take on its reverse complement.
1365
1366         Trying to reverse complement a protein sequence raises an exception.
1367
1368         No return value.
1369         """
1370         self.complement()
1371         self.data.reverse()
1372
1373     ## Sorting a sequence makes no sense.
1374     # def sort(self, *args): self.data.sort(*args)
1375
1376     def extend(self, other):
1377         if isinstance(other, MutableSeq):
1378             for c in other.data:
1379                 self.data.append(c)
1380         else:
1381             for c in other:
1382                 self.data.append(c)
1383
1384     def tostring(self):
1385         """Returns the full sequence as a python string.
1386
1387         Although not formally deprecated, you are now encouraged to use
1388         str(my_seq) instead of my_seq.tostring().
1389
1390         Because str(my_seq) will give you the full sequence as a python string,
1391         there is often no need to make an explicit conversion.  For example,
1392
1393         print "ID={%s}, sequence={%s}" % (my_name, my_seq)
1394
1395         On Biopython 1.44 or older you would have to have done this:
1396
1397         print "ID={%s}, sequence={%s}" % (my_name, my_seq.tostring())
1398         """
1399         return "".join(self.data)
1400
1401     def toseq(self):
1402         """Returns the full sequence as a new immutable Seq object.
1403
1404         >>> from Bio.Seq import Seq
1405         >>> from Bio.Alphabet import IUPAC
1406         >>> my_mseq = MutableSeq("MKQHKAMIVALIVICITAVVAAL", \
1407                                  IUPAC.protein)
1408         >>> my_mseq
1409         MutableSeq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
1410         >>> my_mseq.toseq()
1411         Seq('MKQHKAMIVALIVICITAVVAAL', IUPACProtein())
1412
1413         Note that the alphabet is preserved.
1414         """
1415         return Seq("".join(self.data), self.alphabet)
1416
1417 # The transcribe, backward_transcribe, and translate functions are
1418 # user-friendly versions of the corresponding functions in Bio.Transcribe
1419 # and Bio.Translate. The functions work both on Seq objects, and on strings.
1420
1421 def transcribe(dna):
1422     """Transcribes a DNA sequence into RNA.
1423
1424     If given a string, returns a new string object.
1425
1426     Given a Seq or MutableSeq, returns a new Seq object with an RNA alphabet.
1427
1428     Trying to transcribe a protein or RNA sequence raises an exception.
1429
1430     e.g.
1431
1432     >>> transcribe("ACTGN")
1433     'ACUGN'
1434     """
1435     if isinstance(dna, Seq) :
1436         return dna.transcribe()
1437     elif isinstance(dna, MutableSeq):
1438         return dna.toseq().transcribe()
1439     else:
1440         return dna.replace('T','U').replace('t','u')
1441
1442 def back_transcribe(rna):
1443     """Back-transcribes an RNA sequence into DNA.
1444
1445     If given a string, returns a new string object.
1446
1447     Given a Seq or MutableSeq, returns a new Seq object with an RNA alphabet.
1448
1449     Trying to transcribe a protein or DNA sequence raises an exception.
1450
1451     e.g.
1452
1453     >>> back_transcribe("ACUGN")
1454     'ACTGN'
1455     """
1456     if isinstance(rna, Seq) :
1457         return rna.back_transcribe()
1458     elif isinstance(rna, MutableSeq):
1459         return rna.toseq().back_transcribe()
1460     else:
1461         return rna.replace('U','T').replace('u','t')
1462
1463 def _translate_str(sequence, table, stop_symbol="*",
1464                    to_stop=False, pos_stop="X") :
1465     """Helper function to translate a nucleotide string (PRIVATE).
1466
1467     Arguments:
1468      - sequence    - a string
1469      - table       - a CodonTable object (NOT a table name or id number)
1470      - stop_symbol - a single character string, what to use for terminators.
1471      - to_stop     - boolean, should translation terminate at the first
1472                      in frame stop codon?  If there is no in-frame stop codon
1473                      then translation continues to the end.
1474      - pos_stop    - a single character string for a possible stop codon
1475                      (e.g. TAN or NNN)
1476
1477     Returns a string.
1478
1479     e.g.
1480
1481     >>> from Bio.Data import CodonTable
1482     >>> table = CodonTable.ambiguous_dna_by_id[1]
1483     >>> _translate_str("AAA", table)
1484     'K'
1485     >>> _translate_str("TAR", table)
1486     '*'
1487     >>> _translate_str("TAN", table)
1488     'X'
1489     >>> _translate_str("TAN", table, pos_stop="@")
1490     '@'
1491     >>> _translate_str("TA?", table)
1492     Traceback (most recent call last):
1493        ...
1494     TranslationError: Codon 'TA?' is invalid
1495     """
1496     sequence = sequence.upper()
1497     amino_acids = []
1498     forward_table = table.forward_table
1499     stop_codons = table.stop_codons
1500     if table.nucleotide_alphabet.letters is not None :
1501         valid_letters = set(table.nucleotide_alphabet.letters.upper())
1502     else :
1503         #Assume the worst case, ambiguous DNA or RNA:
1504         valid_letters = set(IUPAC.ambiguous_dna.letters.upper() + \
1505                             IUPAC.ambiguous_rna.letters.upper())
1506
1507     n = len(sequence)
1508     for i in xrange(0,n-n%3,3) :
1509         codon = sequence[i:i+3]
1510         try :
1511             amino_acids.append(forward_table[codon])
1512         except (KeyError, CodonTable.TranslationError) :
1513             #Todo? Treat "---" as a special case (gapped translation)
1514             if codon in table.stop_codons :
1515                 if to_stop : break
1516                 amino_acids.append(stop_symbol)
1517             elif valid_letters.issuperset(set(codon)) :
1518                 #Possible stop codon (e.g. NNN or TAN)
1519                 amino_acids.append(pos_stop)
1520             else :
1521                 raise CodonTable.TranslationError(\
1522                     "Codon '%s' is invalid" % codon)
1523     return "".join(amino_acids)
1524
1525 def translate(sequence, table="Standard", stop_symbol="*", to_stop=False):
1526     """Translate a nucleotide sequence into amino acids.
1527
1528     If given a string, returns a new string object. Given a Seq or
1529     MutableSeq, returns a Seq object with a protein alphabet.
1530
1531     Arguments:
1532      - table - Which codon table to use?  This can be either a name
1533                (string) or an NCBI identifier (integer).  Defaults
1534                to the "Standard" table.
1535      - stop_symbol - Single character string, what to use for any
1536                      terminators, defaults to the asterisk, "*".
1537      - to_stop - Boolean, defaults to False meaning do a full
1538                  translation continuing on past any stop codons
1539                  (translated as the specified stop_symbol).  If
1540                  True, translation is terminated at the first in
1541                  frame stop codon (and the stop_symbol is not
1542                  appended to the returned protein sequence).
1543
1544     A simple string example using the default (standard) genetic code:
1545
1546     >>> coding_dna = "GTGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"
1547     >>> translate(coding_dna)
1548     'VAIVMGR*KGAR*'
1549     >>> translate(coding_dna, stop_symbol="@")
1550     'VAIVMGR@KGAR@'
1551     >>> translate(coding_dna, to_stop=True)
1552     'VAIVMGR'
1553
1554     Now using NCBI table 2, where TGA is not a stop codon:
1555
1556     >>> translate(coding_dna, table=2)
1557     'VAIVMGRWKGAR*'
1558     >>> translate(coding_dna, table=2, to_stop=True)
1559     'VAIVMGRWKGAR'
1560
1561     Note that if the sequence has no in-frame stop codon, then the to_stop
1562     argument has no effect:
1563
1564     >>> coding_dna2 = "GTGGCCATTGTAATGGGCCGC"
1565     >>> translate(coding_dna2)
1566     'VAIVMGR'
1567     >>> translate(coding_dna2, to_stop=True)
1568     'VAIVMGR'
1569
1570     NOTE - Ambiguous codons like "TAN" or "NNN" could be an amino acid
1571     or a stop codon.  These are translated as "X".  Any invalid codon
1572     (e.g. "TA?" or "T-A") will throw a TranslationError.
1573
1574     NOTE - Does NOT support gapped sequences.
1575
1576     It will however translate either DNA or RNA.
1577     """
1578     if isinstance(sequence, Seq) :
1579         return sequence.translate(table, stop_symbol, to_stop)
1580     elif isinstance(sequence, MutableSeq):
1581         #Return a Seq object
1582         return sequence.toseq().translate(table, stop_symbol, to_stop)
1583     else:
1584         #Assume its a string, return a string
1585         try :
1586             codon_table = CodonTable.ambiguous_generic_by_id[int(table)]
1587         except ValueError :
1588             codon_table = CodonTable.ambiguous_generic_by_name[table]
1589         return _translate_str(sequence, codon_table, stop_symbol, to_stop)
1590
1591 def reverse_complement(sequence):
1592     """Returns the reverse complement sequence of a nucleotide string.
1593
1594     If given a string, returns a new string object.
1595     Given a Seq or a MutableSeq, returns a new Seq object with the same alphabet.
1596
1597     Supports unambiguous and ambiguous nucleotide sequences.
1598
1599     e.g.
1600
1601     >>> reverse_complement("ACTG-NH")
1602     'DN-CAGT'
1603     """
1604     if isinstance(sequence, Seq) :
1605         #Return a Seq
1606         return sequence.reverse_complement()
1607     elif isinstance(sequence, MutableSeq) :
1608         #Return a Seq
1609         #Don't use the MutableSeq reverse_complement method as it is 'in place'.
1610         return sequence.toseq().reverse_complement()
1611
1612     #Assume its a string.
1613     #In order to avoid some code duplication, the old code would turn the string
1614     #into a Seq, use the reverse_complement method, and convert back to a string.
1615     #This worked, but is over five times slower on short sequences!
1616     if ('U' in sequence or 'u' in sequence) \
1617     and ('T' in sequence or 't' in sequence):
1618         raise ValueError("Mixed RNA/DNA found")
1619     elif 'U' in sequence or 'u' in sequence:
1620         ttable = _rna_complement_table
1621     else:
1622         ttable = _dna_complement_table
1623     return sequence.translate(ttable)[::-1]
1624
1625 def _test():
1626     """Run the Bio.Seq module's doctests."""
1627     print "Runing doctests..."
1628     import doctest
1629     doctest.testmod()
1630     print "Done"
1631
1632 if __name__ == "__main__":
1633     _test()