binaries/src/disembl/DisEMBL.py

   1 #!/usr/bin/python
   2 # Copyright (C) 2004 Rune Linding & Lars Juhl Jensen - EMBL
   3 # The DisEMBL is licensed under the GPL license
   4 # (http://www.opensource.org/licenses/gpl-license.php)
   5 # DisEMBL pipeline
   6
   7
   8 from string import *
   9 import sys,re
  10 from sys import argv
  11 import os
  12 from os import system,popen3
  13
  14 relpath = re.sub("/DisEMBL.py$","",argv[0])
  15 # cwd = re.sub("/$","", os.getcwd())
  16 newpath = relpath + "/biopython-1.50"
  17 sys.path.append(newpath)
  18
  19 from Bio import File
  20 from Bio import Fasta
  21
  22 import fpformat
  23 import tempfile
  24
  25
  26 # change these to the correct paths
  27 NN_bin = relpath + "/" + 'disembl'
  28 SG_bin = relpath + "/" + 'sav_gol'
  29
  30 def JensenNet(sequence):
  31     outFile = tempfile.mktemp()
  32     inFile= tempfile.mktemp()
  33     open(inFile,'w').write(sequence+'\n')
  34     system(NN_bin + '< ' + inFile +' > ' + outFile)
  35     REM465 = []
  36     COILS = []
  37     HOTLOOPS = []
  38     resultsFile = open(outFile,'r')
  39     results = resultsFile.readlines()
  40     resultsFile.close()
  41     for result in results:
  42         coil = float(fpformat.fix(split(result)[0],6))
  43         COILS.append(coil)
  44         hotloop = float(fpformat.fix(split(result)[1],6))
  45         HOTLOOPS.append(hotloop)
  46         rem465 = float(fpformat.fix(split(result)[2],6))
  47         REM465.append(rem465)
  48     os.remove(inFile)
  49     os.remove(outFile)
  50     return COILS, HOTLOOPS, REM465
  51
  52
  53 def SavitzkyGolay(window,derivative,datalist):
  54     if len(datalist) < 2*window:
  55         window = len(datalist)/2
  56     elif window == 0:
  57         window = 1
  58     stdin, stdout, stderr = popen3(SG_bin + ' -V0 -D' + str(derivative) + ' -n' + str(window)+','+str(window))
  59     for data in datalist:
  60         stdin.write(`data`+'\n')
  61     try:
  62         stdin.close()
  63     except:
  64         print stderr.readlines()
  65     results = stdout.readlines()
  66     stdout.close()
  67     SG_results = []
  68     for result in results:
  69         f = float(fpformat.fix(result,6))
  70         if f < 0:
  71             SG_results.append(0)
  72         else:
  73             SG_results.append(f)
  74     return SG_results
  75
  76 def getSlices(NNdata, fold, join_frame, peak_frame, expect_val):
  77     slices = []
  78     inSlice = 0
  79     for i in range(len(NNdata)):
  80         if inSlice:
  81             if NNdata[i] < expect_val:
  82                 if maxSlice >= fold*expect_val:
  83                     slices.append([beginSlice, endSlice])
  84                 inSlice = 0
  85             else:
  86                 endSlice += 1
  87                 if NNdata[i] > maxSlice:
  88                     maxSlice = NNdata[i]
  89         elif NNdata[i] >= expect_val:
  90             beginSlice = i
  91             endSlice = i
  92             inSlice = 1
  93             maxSlice = NNdata[i]
  94     if inSlice and maxSlice >= fold*expect_val:
  95         slices.append([beginSlice, endSlice])
  96
  97     i = 0
  98     while i < len(slices):
  99         if i+1 < len(slices) and slices[i+1][0]-slices[i][1] <= join_frame:
 100             slices[i] = [ slices[i][0], slices[i+1][1] ]
 101             del slices[i+1]
 102         elif slices[i][1]-slices[i][0]+1 < peak_frame:
 103             del slices[i]
 104         else:
 105             i += 1
 106     return slices
 107
 108
 109 def reportSlicesTXT(slices, sequence):
 110     if slices == []:
 111         s = lower(sequence)
 112     else:
 113         if slices[0][0] > 0:
 114             s = lower(sequence[0:slices[0][0]])
 115         else:
 116             s = ''
 117         for i in range(len(slices)):
 118             if i > 0:
 119                 sys.stdout.write(', ')
 120             sys.stdout.write( str(slices[i][0]+1) + '-' + str(slices[i][1]+1) )
 121
 122     print ''
 123
 124
 125
 126 def runDisEMBLpipeline():
 127     try:
 128         smooth_frame = 8
 129         peak_frame = 8
 130         join_frame = 4
 131         fold_coils = 1.2
 132         fold_hotloops = 1.4
 133         fold_rem465 = 1.2
 134         mode = 'scores'
 135         try:
 136             file = open(sys.argv[1],'r')
 137         except:
 138             mode = 'default'
 139     except:
 140         print '\nDisEMBL.py sequence_file \n'
 141         print 'A default run would be: ./DisEMBL.py fasta_file'
 142         raise SystemExit
 143     #db = sys.stdin
 144     parser = Fasta.RecordParser()
 145     iterator = Fasta.Iterator(file,parser)
 146     while 1:
 147         try:
 148             cur_record = iterator.next()
 149             sequence = upper(cur_record.sequence)
 150             # Run NN
 151             COILS_raw, HOTLOOPS_raw, REM465_raw = JensenNet(sequence)
 152             # Run Savitzky-Golay
 153             REM465_smooth = SavitzkyGolay(smooth_frame,0,REM465_raw)
 154             COILS_smooth = SavitzkyGolay(smooth_frame,0,COILS_raw)
 155             HOTLOOPS_smooth = SavitzkyGolay(smooth_frame,0,HOTLOOPS_raw)
 156
 157             sys.stdout.write('> '+cur_record.title+'\n')
 158             sys.stdout.write('# COILS ')
 159             reportSlicesTXT( getSlices(COILS_smooth, fold_coils, join_frame, peak_frame, 0.43), sequence )
 160             sys.stdout.write('# REM465 ')
 161             reportSlicesTXT( getSlices(REM465_smooth, fold_rem465, join_frame, peak_frame, 0.50), sequence )
 162             sys.stdout.write('# HOTLOOPS ')
 163             reportSlicesTXT( getSlices(HOTLOOPS_smooth, fold_hotloops, join_frame, peak_frame, 0.086), sequence )
 164             sys.stdout.write('# RESIDUE COILS REM465 HOTLOOPS\n')
 165             for i in range(len(REM465_smooth)):
 166                 sys.stdout.write(sequence[i]+'\t'+fpformat.fix(COILS_smooth[i],5)+'\t'+fpformat.fix(REM465_smooth[i],5)+'\t'+fpformat.fix(HOTLOOPS_smooth[i],5)+'\n')
 167         except AttributeError:
 168             break
 169     file.close()
 170     return
 171
 172 runDisEMBLpipeline()