Index: acs/dbdm/dbdm_4/parse-fasta.py
===================================================================
--- /liacs/dbdm/dbdm_4/parse-fasta.py	(revision 58)
+++ 	(revision )
@@ -1,105 +1,0 @@
-#!/usr/bin/env python
-#
-# Parse 2 FASTA files and print statistics
-# BSDLicence
-# Rick van der Zwet - 0433373 - <info@rickvaderzwet.nl>
-from Bio import SeqIO,Seq
-from Bio import Alphabet
-from Bio.Alphabet.IUPAC import ambiguous_dna,unambiguous_dna
-import Bio.Data.CodonTable
-from MultiReplace  import MultiReplace
-
-fasta_translate = { 
-    'r' : 'ga', # purine
-    'y' : 'tc', # pyrimide
-    'k' : 'gt', # keto
-    'm' : 'ac', # amino 
-    's' : 'gc', # strong
-    'w' : 'at', # weak
-    'b' : 'gtc',
-    'd' : 'gat',
-    'h' : 'act',
-    'v' : 'gca',
-    }
-    
-fasta = MultiReplace(fasta_translate)
-
-def parse_file(file):
-    handle = open(file, "rU")
-    for seq_record in SeqIO.parse(handle, "fasta",ambiguous_dna):
-        # How to translate damm thing into plain nucleic acid codes
-        # http://en.wikipedia.org/wiki/FASTA_format
-        retval = seq_record.seq.__str__()
-    
-    handle.close()
-    return(retval)
-
-def stats(seq):
-    pdict = {}
-    for n in range(1, len(seq)):
-        protein = seq[n]
-        if not pdict.has_key(protein):
-            pdict[protein] = 1
-        else:
-            pdict[protein] += 1
-    
-    print pdict
-
-
-def concat(head,tail,max_length=1000):
-    "Concat two strings together removing common parts"
-    l_head = len(head)
-    l_tail = len(tail)
-
-    # Start/Stop at the right time
-    start = 1
-    if (l_head < l_tail):
-        stop = l_head + 1
-    else:
-        stop = l_tail + 1
-
-    # Make sure not to run for-ever
-    if (stop > max_length):
-        stop = max_length
-
-    # Find largest common part
-    # XXX: Not very effient (on very large overlap sets
-    for i in reversed(range(start,stop)):
-        #print "tail[0:%i] '%s' == '%s'" % (i, head[-i:], tail[0:i])
-        if head[-i:] == tail[0:i]:
-            return((i,(tail[0:i]),head + tail[i:]))
-
-    # None found return full part
-    return(-1,'',head + tail)
-
-# File names
-file1 = "data/AE005174v2-1.fas"
-file2 = "data/AE005174v2-2.fas"
-
-# Get data
-seq1 = parse_file(file1)
-seq2 = parse_file(file1)
-
-seq1 = fasta.replace(seq1)
-seq2 = fasta.replace(seq2)
-
-# Find overlap
-(retval, common, result) = concat(seq2,seq1)
-print retval, common
-
-print file1
-stats(seq1)
-print file1
-stats(seq1)
-
-# Strictly speaking there is a gap of about 4 kbs (4000 bs) between seq1 and
-# seq2, so lets' put that into the the statistics as well. Due to circular
-# nature, does not matter wether we add it in the beginning or in the end
-print "Total (inc 4 kbs gap (n))"
-result = result + "n" * 4000;
-stats(result)
-
-# Write to file for later further processing
-out = open("full_contig.raw","w")
-out.write(result)
-out.close()
Index: /liacs/dbdm/dbdm_4/parse_fasta.py
===================================================================
--- /liacs/dbdm/dbdm_4/parse_fasta.py	(revision 59)
+++ /liacs/dbdm/dbdm_4/parse_fasta.py	(revision 59)
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+#
+# Parse 2 FASTA files and print statistics
+# BSDLicence
+# Rick van der Zwet - 0433373 - <info@rickvaderzwet.nl>
+from Bio import SeqIO,Seq
+from Bio import Alphabet
+from Bio.Alphabet.IUPAC import ambiguous_dna,unambiguous_dna
+import Bio.Data.CodonTable
+from MultiReplace  import MultiReplace
+from reading_frames import reading_frames
+
+fasta_translate = { 
+    'r' : 'ga', # purine
+    'y' : 'tc', # pyrimide
+    'k' : 'gt', # keto
+    'm' : 'ac', # amino 
+    's' : 'gc', # strong
+    'w' : 'at', # weak
+    'b' : 'gtc',
+    'd' : 'gat',
+    'h' : 'act',
+    'v' : 'gca',
+    }
+    
+fasta = MultiReplace(fasta_translate)
+
+def parse_file(file):
+    handle = open(file, "rU")
+    for seq_record in SeqIO.parse(handle, "fasta",ambiguous_dna):
+        # How to translate damm thing into plain nucleic acid codes
+        # http://en.wikipedia.org/wiki/FASTA_format
+        retval = seq_record.seq.__str__()
+    
+    handle.close()
+    return(retval)
+
+def stats(seq):
+    pdict = {}
+    for n in range(1, len(seq)):
+        protein = seq[n]
+        if not pdict.has_key(protein):
+            pdict[protein] = 1
+        else:
+            pdict[protein] += 1
+    
+    print pdict
+
+
+def concat(head,tail,max_length=1000):
+    "Concat two strings together removing common parts"
+    l_head = len(head)
+    l_tail = len(tail)
+
+    # Start/Stop at the right time
+    start = 1
+    if (l_head < l_tail):
+        stop = l_head + 1
+    else:
+        stop = l_tail + 1
+
+    # Make sure not to run for-ever
+    if (stop > max_length):
+        stop = max_length
+
+    # Find largest common part
+    # XXX: Not very effient (on very large overlap sets
+    for i in reversed(range(start,stop)):
+        #print "tail[0:%i] '%s' == '%s'" % (i, head[-i:], tail[0:i])
+        if head[-i:] == tail[0:i]:
+            return((i,(tail[0:i]),head + tail[i:]))
+
+    # None found return full part
+    return(-1,'',head + tail)
+
+# File names
+file1 = "data/AE005174v2-1.fas"
+file2 = "data/AE005174v2-2.fas"
+
+# Get data
+seq1 = parse_file(file1)
+seq2 = parse_file(file1)
+
+seq1 = fasta.replace(seq1)
+seq2 = fasta.replace(seq2)
+
+# Find overlap
+(retval, common, result) = concat(seq2,seq1)
+print retval, common
+
+print file1
+stats(seq1)
+reading_frames(seq1)
+print file2
+stats(seq2)
+reading_frames(seq2)
+
+
+# Strictly speaking there is a gap of about 4 kbs (4000 bs) between seq1 and
+# seq2, so lets' put that into the the statistics as well. Due to circular
+# nature, does not matter wether we add it in the beginning or in the end
+print "Total (inc 4 kbs gap (n))"
+result = result + "n" * 4000;
+stats(result)
+reading_frames(result)
+
+# Write to file for later further processing
+out = open("full_contig.raw","w")
+out.write(result)
+out.close()
Index: /liacs/dbdm/dbdm_4/reading_frames.py
===================================================================
--- /liacs/dbdm/dbdm_4/reading_frames.py	(revision 59)
+++ /liacs/dbdm/dbdm_4/reading_frames.py	(revision 59)
@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+#
+# Parse FASTA file and print statistics on reading frames
+# BSDLicence
+# Rick van der Zwet - 0433373 - <info@rickvaderzwet.nl>
+import sys
+
+def _frame_stat(seq, start=0):
+    pdict = {}
+    for n in range(start,len(seq),2):
+        codon = seq[n:n+3]
+        if len(codon) < 3:
+            continue
+        if not pdict.has_key(codon):
+            pdict[codon] = 1
+        else:
+            pdict[codon] += 1
+
+    return(pdict)
+
+def reading_frames(seq):
+    '''Parse from left to right and right to left, at position 1,2,3 in
+    in the so called nucleotide triplets 
+    See: http://en.wikipedia.org/wiki/Genetic_code'''
+
+    # Populate all empty
+    final = {}
+    for a in 'acgtn':
+        for b in 'acgtn':
+            for c in 'acgtn': 
+                final[a + b + c] = [0,0,0,0,0,0]
+
+    for start in [0,1,2]:
+        print "Normal; start %i" % (start)
+        retval = _frame_stat(seq,start)
+        for codon,v in retval.iteritems():
+            final[codon][start] += v
+        print "Reverse; start %i" % (start)
+        retval = _frame_stat(seq[::-1],start)
+        for codon,v in retval.iteritems():
+            final[codon][start+3] += v
+
+    print "CODON :  N:0   , N:1  , N:2  ,  R:0 , R:1  , R:2  "
+
+    for codon in sorted(final.keys()):
+        print codon,"  : ", ",".join(["%6i" % x for x in final[codon]])
+
+
+
+if __name__ == "__main__":
+    # Load data
+    try:
+        handle = open(sys.argv[1],"rU")
+    except IndexError:
+        print "Usage %s <data_file>" % (sys.argv[0])
+        sys.exit(64)
+    except IOError:
+        print "Unable to open '%s'" % (sys.argv[1])
+        sys.exit(64)
+    
+    seq = handle.read()
+    handle.close()
+    reading_frames(seq)
