Changeset 27

Timestamp:

06/10/08 23:46:02 (4 years ago)

Author:

mauro

Message:

 

Location:

biocomp.pscoils/branches/biopython-enabled

Files:

• scripts/pscoils (modified) (3 diffs)
• src/biocomp/pscoils/__init__.py (modified) (5 diffs)
• src/biocomp/pscoils/utils.py (added)
• tests/tests_coils.py (modified) (3 diffs)

biocomp.pscoils/branches/biopython-enabled/scripts/pscoils

@@ -2 +2 @@
 import sys
 import biocomp.pscoils
+from biocomp.pscoils import pred_coil
+from biocomp.pscoils.utils import readFasta, printCoil
 import getopt
 import string
+from Bio.SeqRecord import SeqRecord
 
 def usage(prog):
@@ -31 +34 @@
     for o,a in opts:
         if o == '-f':
-            seqs=biocomp.pscoils.readFasta(a)
+            seqs=readFasta(a)
             modefasta=True
         elif o == '-p':
@@ -55 +58 @@
         #use both 
         for name in seqs.keys(): # WARNING We assume only one sequence!!!
-            gg,gcc,prob,hept_seq,score=biocomp.pscoils.pred_coil2(seqs[name],prof,params,wlambda) 
-            biocomp.pscoils.printCoil(seqs[name],len(seqs[name]),hept_seq,score,prob,gcc,gg,labels)
+                # TODO
+            raise("TODO")
+            gg,gcc,prob,hept_seq,score=biocomp.pscoils.pred_coil2(seqs[name],prof,params,wlambda) 
+            printCoil(seqs[name],len(seqs[name]),hept_seq,score,prob,gcc,gg,labels)
     elif modefasta:
         for name in seqs.keys():
-            gg,gcc,prob,hept_seq,score=biocomp.pscoils.pred_coil(seqs[name],len(seqs[name]),params,biocomp.pscoils.seqScore)
-            biocomp.pscoils.printCoil(seqs[name],len(seqs[name]),hept_seq,score,prob,gcc,gg,labels)
+                seqr = SeqRecord(seqs[name], id=name)
+                seqr = pred_coil(seqr, params)
+                printCoil(seqr,labels)
     elif modeprofile:
+        # TODO
+        raise("TODO")
         gg,gcc,prob,hept_seq,score=biocomp.pscoils.pred_coil(prof,profLen,params,biocomp.pscoils.profScore)
         biocomp.pscoils.printCoil('x'*profLen,profLen,hept_seq,score,prob,gcc,gg,labels)

biocomp.pscoils/branches/biopython-enabled/src/biocomp/pscoils/__init__.py

@@ -35 +35 @@
 import math
 import string
+
 import Params # this must be in the sys.path. Params contains the COILS paramters
+import utils
+
+from Bio.SeqFeature import FeatureLocation, SeqFeature
 
 #
@@ -69 +73 @@
     return prof,pos
     
-#---------------------------------#
-
-
-#----------- readFasta files 
-def readFasta(fname):
-    ''' readFasta(fname) returns a dictionary 
-    '''
-    try:
-        lines=open(fname,'r').readlines()
-    except:
-        sys.stderr.write('cannot open file '+fname+'\n')
-        sys.exit()
-    seqs={}
-    name=''
-    seq=''
-    for line in lines:
-        if line[0]=='>': # assuming fasta file
-           if name != '':
-               seqs[name]=seq
-               seq=''
-           name=line[1:].strip()
-        else:
-           seq+=''.join(line.split())
-    if name != '':
-        seqs[name]=seq
-    return seqs
 #---------------------------------#
 
@@ -134 +112 @@
 
     
-#---------------------------------#
-def pred_coil(seq,seqLen,params,fScore):
+def pred_coil(seqr, params, fScore=None):
     ''' pred_coil(seq,seqLen,params,fScore) returns the coiled coil prediction of sequence seq'''
+    if fScore == None:
+        fScore = seqScore
+    seq = seqr.seq
+    seqLen = len(seqr.seq)
     hept_pos=['a','b','c','d','e','f','g']
     score=[0.0]*seqLen
-    prob=[0.0]*seqLen
-    gcc=[0.0]*seqLen
-    gg=[0.0]*seqLen
     hept_seq=['x']*seqLen
     for i in range(seqLen-params.win+1):
@@ -160 +138 @@
                 hept_seq[i+j]=hept_pos[pos]
     for i in range(seqLen):
-        gg[i],gcc[i],prob[i]= coilProb(score[i],params)
-    return gg,gcc,prob,hept_seq,score
-  
+        gg, gcc, prob = coilProb(score[i],params)
+        seqf = SeqFeature(location=FeatureLocation(i,i), type="pscoils")
+        seqf.qualifiers = {'gg':gg, 
+            'gcc': gcc, 
+            'prob': prob, 
+            'score': score[i],
+            'hept_seq': hept_seq[i]}
+        seqr.features.append(seqf)
+    return seqr  
     
 #---------------------------------#
@@ -195 +179 @@
     return gg,gcc,prob,hept_seq,score
   
-def printCoil(seq,seqLen,hept_seq,score,P,Gcc,Gg,labels='T',io=sys.stdout):
-    ''' '''
-    ccLab="C"
-    loopLab="L"
-    if labels !='T' :
-        for i in range(seqLen):
-            #print "%4d %c %c %7.3f %7.3f (%7.3f %7.3f)" % (i+1,seq[i],hept_seq[i],score[i],P[i],Gcc[i],Gg[i])
-            io.write("%4d %c %c %7.3f %7.3f %7.3f %7.3f\n" % (i+1,seq[i],hept_seq[i],score[i],P[i],Gcc[i],Gg[i]))
-    else:
-        io.write(" Pos A Hep Score   Prob    Gcc     Gg    Pred (Loop=L Coiledcoil=C)\n")
-        for i in range(seqLen):
-            if P[i] > 0.5:
-               clabel=ccLab
-            else:
-               clabel=loopLab
-            io.write("%4d %c %c %7.3f %7.3f %7.3f %7.3f %s\n" % (i+1,seq[i],hept_seq[i],score[i],P[i],Gcc[i],Gg[i],clabel))
             

biocomp.pscoils/branches/biopython-enabled/tests/tests_coils.py

@@ -3 +3 @@
 from unittest import TestCase
 import StringIO
-import biocomp.pscoils
+from biocomp.pscoils import pred_coil
+from biocomp.pscoils.utils import readFasta, printCoil
+from biocomp.pscoils.Params import Params
+from Bio.SeqRecord import SeqRecord
 
 class test_suite(TestCase):
@@ -12 +15 @@
         win=21
         modeprofile=False
-        seqs=biocomp.pscoils.readFasta("tests/ferritin.fasta")
+        seqs=readFasta("tests/ferritin.fasta")
         modefasta=True
         labels='F'
-        params=biocomp.pscoils.Params.Params(weight,win)
+        params=Params(weight,win)
         output = StringIO.StringIO()
         for name in seqs.keys():
-            gg,gcc,prob,hept_seq,score=biocomp.pscoils.pred_coil(seqs[name],len(seqs[name]),params,biocomp.pscoils.seqScore)
-            biocomp.pscoils.printCoil(seqs[name],len(seqs[name]),hept_seq,score,prob,gcc,gg,labels,io=output)
+            seqr = SeqRecord(seqs[name], id=name)
+            seqr = pred_coil(seqr, params)
+            printCoil(seqr,labels,io=output)
         bench_output = open("tests/ferritin.Coils").read()
         self.assertEqual(output.getvalue(),bench_output)        
@@ -29 +33 @@
         win=21
         modeprofile=False
-        seqs=biocomp.pscoils.readFasta("tests/ceba_bovin.fasta")
+        seqs=readFasta("tests/ceba_bovin.fasta")
         modefasta=True
         labels='F'
-        params=biocomp.pscoils.Params.Params(weight,win)
+        params=Params(weight,win)
         output = StringIO.StringIO()
         for name in seqs.keys():
-            gg,gcc,prob,hept_seq,score=biocomp.pscoils.pred_coil(seqs[name],len(seqs[name]),params,biocomp.pscoils.seqScore)
-            biocomp.pscoils.printCoil(seqs[name],len(seqs[name]),hept_seq,score,prob,gcc,gg,labels,io=output)
+            seqr = SeqRecord(seqs[name], id=name)
+            seqr = pred_coil(seqr, params)
+            printCoil(seqr,labels,io=output)
         bench_output = open("tests/ceba_bovin.Coils").read()
         self.assertEqual(output.getvalue(),bench_output)