--- /dev/null
+#!/usr/bin/env python
+# Created: Wed May 29 08:07:18 2002
+# thomas@cbs.dtu.dk, Cecilia.Alsmark@ebc.uu.se
+# Copyright 2001 by Thomas Sicheritz-Ponten and Cecilia Alsmark.
+# 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.
+
+"""Miscellaneous functions for dealing with sequences."""
+
+import re, time
+from Bio import SeqIO
+from Bio import Translate
+from Bio.Seq import Seq
+from Bio import Alphabet
+from Bio.Alphabet import IUPAC
+from Bio.Data import IUPACData, CodonTable
+
+
+######################################
+# DNA
+######################
+# {{{
+
+def reverse(seq):
+ """Reverse the sequence. Works on string sequences.
+
+ e.g.
+ >>> reverse("ACGGT")
+ 'TGGCA'
+
+ """
+ r = list(seq)
+ r.reverse()
+ return ''.join(r)
+
+def GC(seq):
+ """Calculates G+C content, returns the percentage (float between 0 and 100).
+
+ Copes mixed case seuqneces, and with the ambiguous nucleotide S (G or C)
+ when counting the G and C content. The percentage is calculated against
+ the full length, e.g.:
+
+ >>> from Bio.SeqUtils import GC
+ >>> GC("ACTGN")
+ 40.0
+
+ Note that this will return zero for an empty sequence.
+ """
+ try :
+ gc = sum(map(seq.count,['G','C','g','c','S','s']))
+ return gc*100.0/len(seq)
+ except ZeroDivisionError :
+ return 0.0
+
+
+def GC123(seq):
+ """Calculates total G+C content plus first, second and third positions.
+
+ Returns a tuple of four floats (percentages between 0 and 100) for the
+ entire sequence, and the three codon positions. e.g.
+
+ >>> from Bio.SeqUtils import GC123
+ >>> GC123("ACTGTN")
+ (40.0, 50.0, 50.0, 0.0)
+
+ Copes with mixed case sequences, but does NOT deal with ambiguous
+ nucleotides.
+ """
+ d= {}
+ for nt in ['A','T','G','C']:
+ d[nt] = [0,0,0]
+
+ for i in range(0,len(seq),3):
+ codon = seq[i:i+3]
+ if len(codon) <3: codon += ' '
+ for pos in range(0,3):
+ for nt in ['A','T','G','C']:
+ if codon[pos] == nt or codon[pos] == nt.lower():
+ d[nt][pos] += 1
+ gc = {}
+ gcall = 0
+ nall = 0
+ for i in range(0,3):
+ try:
+ n = d['G'][i] + d['C'][i] +d['T'][i] + d['A'][i]
+ gc[i] = (d['G'][i] + d['C'][i])*100.0/n
+ except:
+ gc[i] = 0
+
+ gcall = gcall + d['G'][i] + d['C'][i]
+ nall = nall + n
+
+ gcall = 100.0*gcall/nall
+ return gcall, gc[0], gc[1], gc[2]
+
+def GC_skew(seq, window = 100):
+ """Calculates GC skew (G-C)/(G+C) for multuple windows along the sequence.
+
+ Returns a list of ratios (floats), controlled by the length of the sequence
+ and the size of the window.
+
+ Does NOT look at any ambiguous nucleotides.
+ """
+ # 8/19/03: Iddo: added lowercase
+ values = []
+ for i in range(0, len(seq), window):
+ s = seq[i: i + window]
+ g = s.count('G') + s.count('g')
+ c = s.count('C') + s.count('c')
+ skew = (g-c)/float(g+c)
+ values.append(skew)
+ return values
+
+from math import pi, sin, cos, log
+def xGC_skew(seq, window = 1000, zoom = 100,
+ r = 300, px = 100, py = 100):
+ """Calculates and plots normal and accumulated GC skew (GRAPHICS !!!)."""
+ from Tkinter import Scrollbar, Canvas, BOTTOM, BOTH, ALL, \
+ VERTICAL, HORIZONTAL, RIGHT, LEFT, X, Y
+ yscroll = Scrollbar(orient = VERTICAL)
+ xscroll = Scrollbar(orient = HORIZONTAL)
+ canvas = Canvas(yscrollcommand = yscroll.set,
+ xscrollcommand = xscroll.set, background = 'white')
+ win = canvas.winfo_toplevel()
+ win.geometry('700x700')
+
+ yscroll.config(command = canvas.yview)
+ xscroll.config(command = canvas.xview)
+ yscroll.pack(side = RIGHT, fill = Y)
+ xscroll.pack(side = BOTTOM, fill = X)
+ canvas.pack(fill=BOTH, side = LEFT, expand = 1)
+ canvas.update()
+
+ X0, Y0 = r + px, r + py
+ x1, x2, y1, y2 = X0 - r, X0 + r, Y0 -r, Y0 + r
+
+ ty = Y0
+ canvas.create_text(X0, ty, text = '%s...%s (%d nt)' % (seq[:7], seq[-7:], len(seq)))
+ ty +=20
+ canvas.create_text(X0, ty, text = 'GC %3.2f%%' % (GC(seq)))
+ ty +=20
+ canvas.create_text(X0, ty, text = 'GC Skew', fill = 'blue')
+ ty +=20
+ canvas.create_text(X0, ty, text = 'Accumulated GC Skew', fill = 'magenta')
+ ty +=20
+ canvas.create_oval(x1,y1, x2, y2)
+
+ acc = 0
+ start = 0
+ for gc in GC_skew(seq, window):
+ r1 = r
+ acc+=gc
+ # GC skew
+ alpha = pi - (2*pi*start)/len(seq)
+ r2 = r1 - gc*zoom
+ x1 = X0 + r1 * sin(alpha)
+ y1 = Y0 + r1 * cos(alpha)
+ x2 = X0 + r2 * sin(alpha)
+ y2 = Y0 + r2 * cos(alpha)
+ canvas.create_line(x1,y1,x2,y2, fill = 'blue')
+ # accumulated GC skew
+ r1 = r - 50
+ r2 = r1 - acc
+ x1 = X0 + r1 * sin(alpha)
+ y1 = Y0 + r1 * cos(alpha)
+ x2 = X0 + r2 * sin(alpha)
+ y2 = Y0 + r2 * cos(alpha)
+ canvas.create_line(x1,y1,x2,y2, fill = 'magenta')
+
+ canvas.update()
+ start += window
+
+ canvas.configure(scrollregion = canvas.bbox(ALL))
+
+def molecular_weight(seq):
+ """Calculate the molecular weight of a DNA sequence."""
+ if type(seq) == type(''): seq = Seq(seq, IUPAC.unambiguous_dna)
+ weight_table = IUPACData.unambiguous_dna_weights
+ #TODO, use a generator expession once we drop Python 2.3?
+ #e.g. return sum(weight_table[x] for x in seq)
+ total = 0
+ for x in seq:
+ total += weight_table[x]
+ return total
+
+def nt_search(seq, subseq):
+ """Search for a DNA subseq in sequence.
+
+ use ambiguous values (like N = A or T or C or G, R = A or G etc.)
+ searches only on forward strand
+ """
+ pattern = ''
+ for nt in subseq:
+ value = IUPACData.ambiguous_dna_values[nt]
+ if len(value) == 1:
+ pattern += value
+ else:
+ pattern += '[%s]' % value
+
+ pos = -1
+ result = [pattern]
+ l = len(seq)
+ while True:
+ pos+=1
+ s = seq[pos:]
+ m = re.search(pattern, s)
+ if not m: break
+ pos += int(m.start(0))
+ result.append(pos)
+ return result
+
+# }}}
+
+######################################
+# Protein
+######################
+# {{{
+
+# temporary hack for exception free translation of "dirty" DNA
+# should be moved to ???
+
+class ProteinX(Alphabet.ProteinAlphabet):
+ letters = IUPACData.extended_protein_letters + "X"
+
+proteinX = ProteinX()
+
+class MissingTable:
+ def __init__(self, table):
+ self._table = table
+ def get(self, codon, stop_symbol):
+ try:
+ return self._table.get(codon, stop_symbol)
+ except CodonTable.TranslationError:
+ return 'X'
+
+def makeTableX(table):
+ assert table.protein_alphabet == IUPAC.extended_protein
+ return CodonTable.CodonTable(table.nucleotide_alphabet, proteinX,
+ MissingTable(table.forward_table),
+ table.back_table, table.start_codons,
+ table.stop_codons)
+
+# end of hacks
+
+def seq3(seq):
+ """Turn a one letter code protein sequence into one with three letter codes.
+
+ The single input argument 'seq' should be a protein sequence using single
+ letter codes, either as a python string or as a Seq or MutableSeq object.
+
+ This function returns the amino acid sequence as a string using the three
+ letter amino acid codes. Output follows the IUPAC standard (including
+ ambiguous characters B for "Asx", J for "Xle" and X for "Xaa", and also U
+ for "Sel" and O for "Pyl") plus "Ter" for a terminator given as an asterisk. Any unknown
+ character (including possible gap characters), is changed into 'Xaa'.
+
+ e.g.
+ >>> from Bio.SeqUtils import seq3
+ >>> seq3("MAIVMGRWKGAR*")
+ 'MetAlaIleValMetGlyArgTrpLysGlyAlaArgTer'
+
+ This function was inspired by BioPerl's seq3.
+ """
+ threecode = {'A':'Ala', 'B':'Asx', 'C':'Cys', 'D':'Asp',
+ 'E':'Glu', 'F':'Phe', 'G':'Gly', 'H':'His',
+ 'I':'Ile', 'K':'Lys', 'L':'Leu', 'M':'Met',
+ 'N':'Asn', 'P':'Pro', 'Q':'Gln', 'R':'Arg',
+ 'S':'Ser', 'T':'Thr', 'V':'Val', 'W':'Trp',
+ 'Y':'Tyr', 'Z':'Glx', 'X':'Xaa', '*':'Ter',
+ 'U':'Sel', 'O':'Pyl', 'J':'Xle',
+ }
+ #We use a default of 'Xaa' for undefined letters
+ #Note this will map '-' to 'Xaa' which may be undesirable!
+ return ''.join([threecode.get(aa,'Xaa') for aa in seq])
+
+
+# }}}
+
+######################################
+# Mixed ???
+######################
+# {{{
+
+def translate(seq, frame = 1, genetic_code = 1, translator = None):
+ """Translation of DNA in one of the six different reading frames (DEPRECATED).
+
+ Use the Bio.Seq.Translate function, or the Seq object's translate method
+ instead:
+
+ >>> from Bio.Seq import Seq
+ >>> my_seq = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG")
+ >>> my_seq = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAGUA")
+ >>> for frame in [0,1,2] :
+ ... print my_seq[frame:].translate()
+ ...
+ MAIVMGR*KGAR*
+ WPL*WAAERVPDS
+ GHCNGPLKGCPIV
+ >>> for frame in [0,1,2] :
+ ... print my_seq.reverse_complement()[frame:].translate()
+ ...
+ YYRAPFQRPITMA
+ TIGHPFSGPLQWP
+ LSGTLSAAHYNGH
+ """
+ import warnings
+ warnings.warn("Bio.SeqUtils.translate() has been deprecated, and we intend" \
+ +" to remove it in a future release of Biopython. Please use"\
+ +" the method or function in Bio.Seq instead, as described in"\
+ +" the Tutorial.", DeprecationWarning)
+
+ if frame not in [1,2,3,-1,-2,-3]:
+ raise ValueError('invalid frame')
+
+ if not translator:
+ table = makeTableX(CodonTable.ambiguous_dna_by_id[genetic_code])
+ translator = Translate.Translator(table)
+
+ #Does this frame calculation do something sensible? No RC taken!
+ return translator.translate(Seq(seq[frame-1:], IUPAC.ambiguous_dna)).data
+
+def GC_Frame(seq, genetic_code = 1):
+ """Just an alias for six_frame_translations (OBSOLETE).
+
+ Use six_frame_translation directly, as this function may be deprecated
+ in a future release."""
+ return six_frame_translations(seq, genetic_code)
+
+def six_frame_translations(seq, genetic_code = 1):
+ """Formatted string showing the 6 frame translations and GC content.
+
+ nice looking 6 frame translation with GC content - code from xbbtools
+ similar to DNA Striders six-frame translation
+
+ e.g.
+ from Bio.SeqUtils import six_frame_translations
+ print six_frame_translations("AUGGCCAUUGUAAUGGGCCGCUGA")
+ """
+ from Bio.Seq import reverse_complement, translate
+ anti = reverse_complement(seq)
+ comp = anti[::-1]
+ length = len(seq)
+ frames = {}
+ for i in range(0,3):
+ frames[i+1] = translate(seq[i:], genetic_code)
+ frames[-(i+1)] = reverse(translate(anti[i:], genetic_code))
+
+ # create header
+ if length > 20:
+ short = '%s ... %s' % (seq[:10], seq[-10:])
+ else:
+ short = seq
+ #TODO? Remove the date as this would spoil any unit test...
+ date = time.strftime('%y %b %d, %X', time.localtime(time.time()))
+ header = 'GC_Frame: %s, ' % date
+ for nt in ['a','t','g','c']:
+ header += '%s:%d ' % (nt, seq.count(nt.upper()))
+
+ header += '\nSequence: %s, %d nt, %0.2f %%GC\n\n\n' % (short.lower(),length, GC(seq))
+ res = header
+
+ for i in range(0,length,60):
+ subseq = seq[i:i+60]
+ csubseq = comp[i:i+60]
+ p = i/3
+ res = res + '%d/%d\n' % (i+1, i/3+1)
+ res = res + ' ' + ' '.join(map(None,frames[3][p:p+20])) + '\n'
+ res = res + ' ' + ' '.join(map(None,frames[2][p:p+20])) + '\n'
+ res = res + ' '.join(map(None,frames[1][p:p+20])) + '\n'
+ # seq
+ res = res + subseq.lower() + '%5d %%\n' % int(GC(subseq))
+ res = res + csubseq.lower() + '\n'
+ # - frames
+ res = res + ' '.join(map(None,frames[-2][p:p+20])) +' \n'
+ res = res + ' ' + ' '.join(map(None,frames[-1][p:p+20])) + '\n'
+ res = res + ' ' + ' '.join(map(None,frames[-3][p:p+20])) + '\n\n'
+ return res
+
+# }}}
+
+######################################
+# FASTA file utilities
+######################
+# {{{
+
+def fasta_uniqids(file):
+ """Checks and changes the name/ID's to be unique identifiers by adding numbers (OBSOLETE).
+
+ file - a FASTA format filename to read in.
+
+ No return value, the output is written to screen.
+ """
+ dict = {}
+ txt = open(file).read()
+ entries = []
+ for entry in txt.split('>')[1:]:
+ name, seq= entry.split('\n',1)
+ name = name.split()[0].split(',')[0]
+
+ if name in dict:
+ n = 1
+ while 1:
+ n = n + 1
+ _name = name + str(n)
+ if _name not in dict:
+ name = _name
+ break
+
+ dict[name] = seq
+
+ for name, seq in dict.items():
+ print '>%s\n%s' % (name, seq)
+
+def quick_FASTA_reader(file):
+ """Simple FASTA reader, returning a list of string tuples.
+
+ The single argument 'file' should be the filename of a FASTA format file.
+ This function will open and read in the entire file, constructing a list
+ of all the records, each held as a tuple of strings (the sequence name or
+ title, and its sequence).
+
+ This function was originally intended for use on large files, where its
+ low overhead makes it very fast. However, because it returns the data as
+ a single in memory list, this can require a lot of RAM on large files.
+
+ You are generally encouraged to use Bio.SeqIO.parse(handle, "fasta") which
+ allows you to iterate over the records one by one (avoiding having all the
+ records in memory at once). Using Bio.SeqIO also makes it easy to switch
+ between different input file formats. However, please note that rather
+ than simple strings, Bio.SeqIO uses SeqRecord objects for each record.
+ """
+ #Want to split on "\n>" not just ">" in case there are any extra ">"
+ #in the name/description. So, in order to make sure we also split on
+ #the first entry, prepend a "\n" to the start of the file.
+ handle = open(file)
+ txt = "\n" + handle.read()
+ handle.close()
+ entries = []
+ for entry in txt.split('\n>')[1:]:
+ name,seq= entry.split('\n',1)
+ seq = seq.replace('\n','').replace(' ','').upper()
+ entries.append((name, seq))
+ return entries
+
+def apply_on_multi_fasta(file, function, *args):
+ """Apply a function on each sequence in a multiple FASTA file (OBSOLETE).
+
+ file - filename of a FASTA format file
+ function - the function you wish to invoke on each record
+ *args - any extra arguments you want passed to the function
+
+ This function will iterate over each record in a FASTA file as SeqRecord
+ objects, calling your function with the record (and supplied args) as
+ arguments.
+
+ This function returns a list. For those records where your function
+ returns a value, this is taken as a sequence and used to construct a
+ FASTA format string. If your function never has a return value, this
+ means apply_on_multi_fasta will return an empty list.
+ """
+ try:
+ f = globals()[function]
+ except:
+ raise NotImplementedError("%s not implemented" % function)
+
+ handle = open(file, 'r')
+ records = SeqIO.parse(handle, "fasta")
+ results = []
+ for record in records:
+ arguments = [record.sequence]
+ for arg in args: arguments.append(arg)
+ result = f(*arguments)
+ if result:
+ results.append('>%s\n%s' % (record.name, result))
+ handle.close()
+ return results
+
+def quicker_apply_on_multi_fasta(file, function, *args):
+ """Apply a function on each sequence in a multiple FASTA file (OBSOLETE).
+
+ file - filename of a FASTA format file
+ function - the function you wish to invoke on each record
+ *args - any extra arguments you want passed to the function
+
+ This function will use quick_FASTA_reader to load every record in the
+ FASTA file into memory as a list of tuples. For each record, it will
+ call your supplied function with the record as a tuple of the name and
+ sequence as strings (plus any supplied args).
+
+ This function returns a list. For those records where your function
+ returns a value, this is taken as a sequence and used to construct a
+ FASTA format string. If your function never has a return value, this
+ means quicker_apply_on_multi_fasta will return an empty list.
+ """
+ try:
+ f = globals()[function]
+ except:
+ raise NotImplementedError("%s not implemented" % function)
+
+ entries = quick_FASTA_reader(file)
+ results = []
+ for name, seq in entries:
+ arguments = [seq]
+ for arg in args: arguments.append(arg)
+ result = f(*arguments)
+ if result:
+ results.append('>%s\n%s' % (name, result))
+ handle.close()
+ return results
+
+# }}}
+
+######################################
+# Main
+#####################
+# {{{
+
+if __name__ == '__main__':
+ import sys, getopt
+ # crude command line options to use most functions directly on a FASTA file
+ options = {'apply_on_multi_fasta':0,
+ 'quick':0,
+ 'uniq_ids':0,
+ }
+
+ optlist, args = getopt.getopt(sys.argv[1:], '', ['describe', 'apply_on_multi_fasta=',
+ 'help', 'quick', 'uniq_ids', 'search='])
+ for arg in optlist:
+ if arg[0] in ['-h', '--help']:
+ pass
+ elif arg[0] in ['--describe']:
+ # get all new functions from this file
+ mol_funcs = [x[0] for x in locals().items() if type(x[1]) == type(GC)]
+ mol_funcs.sort()
+ print 'available functions:'
+ for f in mol_funcs: print '\t--%s' % f
+ print '\n\ne.g.\n./sequtils.py --apply_on_multi_fasta GC test.fas'
+
+ sys.exit(0)
+ elif arg[0] in ['--apply_on_multi_fasta']:
+ options['apply_on_multi_fasta'] = arg[1]
+ elif arg[0] in ['--search']:
+ options['search'] = arg[1]
+ else:
+ key = re.search('-*(.+)', arg[0]).group(1)
+ options[key] = 1
+
+
+ if options.get('apply_on_multi_fasta'):
+ file = args[0]
+ function = options['apply_on_multi_fasta']
+ arguments = []
+ if options.get('search'):
+ arguments = options['search']
+ if function == 'xGC_skew':
+ arguments = 1000
+ if options.get('quick'):
+ results = quicker_apply_on_multi_fasta(file, function, arguments)
+ else:
+ results = apply_on_multi_fasta(file, function, arguments)
+ for result in results: print result
+
+ elif options.get('uniq_ids'):
+ file = args[0]
+ fasta_uniqids(file)
+
+# }}}
+
git://source.jalview.org / jabaws.git / blobdiff