Changeset 53 for liacs/dbdm

Timestamp:

Dec 20, 2009, 8:56:07 PM (15 years ago)

Author:

Rick van der Zwet

Message:

Temponary commit

File:

• liacs/dbdm/dbdm_4/fasta-hmm.py (modified) (1 diff)

liacs/dbdm/dbdm_4/fasta-hmm.py

@@ -4 +4 @@
 from Bio.Alphabet.IUPAC import ambiguous_dna,unambiguous_dna
 import Bio.Data.CodonTable
+from MultiReplace  import MultiReplace
 
+def parse_file(file):
+    handle = open("data/AE005174v2-1.fas", "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
+        stupid = seq_record.seq.__str__()
+        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',
+            }
+            
+        r = MultiReplace(fasta_translate)
+        stupid = r.replace(stupid)
+    
+        pdict = {}
+        for n in range(1, len(stupid)):
+            protein = stupid[n]
+            if not pdict.has_key(protein):
+                pdict[protein] = 1
+            else:
+                pdict[protein] += 1
+    
+        print pdict
 
-handle = open("data/AE005174v2-1.fas", "rU")
-for seq_record in SeqIO.parse(handle, "fasta",ambiguous_dna):
-    print seq_record.id
-    print repr(seq_record.seq)
-    print seq_record.seq.alphabet
-    print seq_record.letter_annotations
-
-    # How to translate damm thing into plain nucleic acid codes
-    # http://en.wikipedia.org/wiki/FASTA_format
-    stupid = seq_record.seq.to_str().translate({'W' : 'G'})
-
-    pdict = {}
-    for n in range(1, len(stupid)):
-        protein = stupid[n]
-        if not pdict.has_key(protein):
-            pdict[protein] = 1
-        else:
-            pdict[protein] += 1
-
-    print pdict
-
-
+file1 = parse_file("data/AE005174v2-1.fas")
+file2 = parse_file("data/AE005174v2-2.fas")