Changeset 62 for liacs

Timestamp:

Dec 23, 2009, 2:35:48 PM (15 years ago)

Author:

Rick van der Zwet

Message:

First results on shifting tests

Location:

liacs/dbdm/dbdm_4

Files:

• ecoli_hmm.py (modified) (2 diffs)
• full_contig.txt.gz (added)
• shifts.txt (added)

liacs/dbdm/dbdm_4/ecoli_hmm.py

@@ -10 +10 @@
 import csv
 import string
+
+dna_flip = string.maketrans('atgc','tacg')
+def get_codon(seq, position, order):
+    codon = seq[position:position+3]
+    if not order:
+        # When living on the other side of the string make sure to flip 
+        # around before comparison
+        codon = codon[::-1].translate(dna_flip)
+    return(codon)
+
 
 
@@ -50 +60 @@
     # rend = Right End
 
-    # Training
+    # Training and testing
     # A -- T, G -- C
-    start_condons = ['atg', 'gtg']
-    stop_condons = ['taa', 'tga', 'tag']
-    dna_flip = string.maketrans('atgc','tacg')
+    start_codons = ['atg', 'gtg']
+    stop_codons = ['taa', 'tga', 'tag']
     edl_file = 'data/edl_genes.csv'
+
+    # Counter limit of how many 'hard' errors are allowed before bailing out
+    stop_limit = 5
+
+    # Current shifting offset
+    base_shift = 0
+    stop_counter = 0
+
     reader = csv.reader(open(edl_file,"rU"))
     reader.next()
     try:
         for row in reader:
-            (featureType, zNumber, contig, lend, rend, orientation, 
-            segmentType, oIslandNumber, geneName, note, function, product, 
-            translationNotes) = row
-            lend = int(lend) - 1
-            rend = int(rend)
-            # Make a forward orientation as positive
-            if orientation == '>':
-                order = True
-            else:
-                order = False
+            finished = False
+            while not finished:
+                (featureType, zNumber, contig, lend, rend, orientation, 
+                segmentType, oIslandNumber, geneName, note, function, product, 
+                translationNotes) = row
+                lend = int(lend)
+                rend = int(rend)
+                # Make a forward orientation mark as positive
+                if orientation == '>':
+                    order = True
+                else:
+                    order = False
 
-            if order:
-                start_codon = seq[lend:lend+3]
-                stop_codon = seq[rend-3:rend]
-            else:
-                start_codon = seq[rend-3:rend].translate(dna_flip)[::-1]
-                stop_codon = seq[lend:lend+3].translate(dna_flip)[::-1]
+                # Living the world upside down to order is off
+                if order:
+                    start_pos = base_shift + lend - 1
+                    stop_pos = base_shift + rend - 3
+                else:
+                    start_pos = base_shift + rend -3
+                    stop_pos = base_shift + lend - 1
+                    
+                start_codon = get_codon(seq,start_pos,order)
+                stop_codon = get_codon(seq,stop_pos,order)
+                start_match = start_codon in start_codons
+                stop_match = stop_codon in stop_codons
 
-            print "%6s, %s, %s, %s, %s, %s," % (geneName,lend,rend, start_codon, stop_codon, orientation),
+                print "%6s, %s, %s, %s, %s, %s," % (geneName,lend,rend, start_codon, stop_codon, orientation),
+                print "%s, %s" % (start_match, stop_match)
+    
+                # Check for fucked up offsets, shifts
+                start_shift = None
+                stop_shift = None 
+                # Technically speaking should the cope only the reading frames,
+                # but what the heck looking if 'free'
+                search_range = 10
+                if not start_match:
+                    # Somewhere else in the reading frame?
+                    for r in range(-search_range,search_range+1):
+                        l = start_pos + r
+                        t = get_codon(seq,l,order)
+                        new_match = t in start_codons
+                        if new_match:
+                            start_shift = r
+                            print "# Start codon reading frame match - shift %i: " % r,
+                            print "%i, %s: %s" % (l,t,new_match)
+                if not stop_match:
+                    # Somewhere else in the reading frame?
+                    for r in range(-search_range,search_range+1):
+                        l = stop_pos + r
+                        t = get_codon(seq,l,order)
+                        new_match = t in stop_codons
+                        if new_match:
+                            stop_shift = r
+                            print "# Stop codon reading frame match - shift %i: " % r,
+                            print "%i, %s: %s" % (l,t,new_match)
 
-            print "%s, %s" % (start_codon in start_condons, stop_codon in stop_condons)
-            if not start_codon in start_condons or not stop_codon in stop_condons:
-                return
+
+                # Both wrong is something screwy in the data, check offset fix
+                if not start_match and not stop_match:
+                    if stop_shift != None and start_shift == stop_shift:
+                        print "# Matching shift %i" % start_shift
+                        base_shift += start_shift
+                        finished = False
+                        continue
+                    # Exercise left to the reader
+                    if (stop_counter > stop_limit):
+                        return
+                    else:
+                        stop_counter += 1
+                        finished = True
+                else:
+                    stop_counter = 0
+                    finished = True
+
     except csv.Error, e:
         sys.exit('file %s, line %d: %s' % (filename, reader.line_num, e))