binaries/src/globplot/GlobPlot.py

   1 #!/usr/bin/env python
   2 # Copyright (C) 2003 Rune Linding - EMBL
   3 # GlobPlot TM
   4 # GlobPlot is licensed under the Academic Free license
   5
   6 from string import *
   7 from sys import argv
   8 import sys,re
   9 import os
  10 from os import system,popen3
  11
  12 relpath = re.sub("/GlobPlot.py$","",argv[0])
  13 # cwd = re.sub("/$","", os.getcwd())
  14 newpath = relpath + "/biopython-1.50"
  15 sys.path.append(newpath)
  16
  17 import Bio
  18 from Bio import File
  19 from Bio import Fasta
  20 import fpformat
  21 import tempfile
  22 import math
  23
  24 # Russell/Linding
  25 RL = {'N':0.229885057471264,'P':0.552316012226663,'Q':-0.187676577424997,'A':-0.261538461538462,'R':-0.176592654077609, \
  26       'S':0.142883029808825,'C':-0.0151515151515152,'T':0.00887797506611258,'D':0.227629796839729,'E':-0.204684629516228, \
  27       'V':-0.386174834235195,'F':-0.225572305974316,'W':-0.243375458622095,'G':0.433225711769886,'H':-0.00121743364986608, \
  28       'Y':-0.20750516775322,'I':-0.422234699606962,'K':-0.100092289621613,'L':-0.337933495925287,'M':-0.225903614457831}
  29
  30 def Sum(seq,par_dict):
  31     sum = 0
  32     results = []
  33     raws = []
  34     sums = []
  35     p = 1
  36     for residue in seq:
  37         try:
  38             parameter = par_dict[residue]
  39         except:
  40             parameter = 0
  41         if p == 1:
  42             sum = parameter
  43         else:
  44             sum = sum + parameter#*math.log10(p)
  45         ssum = float(fpformat.fix(sum,10))
  46         sums.append(ssum)
  47         p +=1
  48     return sums
  49
  50 def getSlices(dydx_data, DOM_join_frame, DOM_peak_frame, DIS_join_frame, DIS_peak_frame):
  51     DOMslices = []
  52     DISslices = []
  53     in_DOMslice = 0
  54     in_DISslice = 0
  55     beginDOMslice = 0
  56     endDOMslice = 0
  57     beginDISslice = 0
  58     endDISslice = 0
  59     for i in range( len(dydx_data) ):
  60     #close dom slice
  61         if in_DOMslice and dydx_data[i] > 0:
  62             DOMslices.append([beginDOMslice, endDOMslice])
  63             in_DOMslice = 0
  64     #close dis slice
  65         elif in_DISslice and dydx_data[i] < 0:
  66             DISslices.append([beginDISslice, endDISslice])
  67             in_DISslice = 0
  68         # elseif inSlice expandslice
  69         elif in_DOMslice:
  70             endDOMslice += 1
  71         elif in_DISslice:
  72             endDISslice += 1
  73     # if not in slice and dydx !== 0 start slice
  74         if dydx_data[i] > 0 and not in_DISslice:
  75             beginDISslice = i
  76             endDISslice = i
  77             in_DISslice = 1
  78         elif dydx_data[i] < 0 and not in_DOMslice:
  79             beginDOMslice = i
  80             endDOMslice = i
  81             in_DOMslice = 1
  82     #last slice
  83     if in_DOMslice:
  84         DOMslices.append([beginDOMslice, endDOMslice])
  85     if in_DISslice:
  86         DISslices.append([beginDISslice,endDISslice])
  87     k = 0
  88     l = 0
  89     while k < len(DOMslices):
  90         if k+1 < len(DOMslices) and DOMslices[k+1][0]-DOMslices[k][1] < DOM_join_frame:
  91             DOMslices[k] = [ DOMslices[k][0], DOMslices[k+1][1] ]
  92             del DOMslices[k+1]
  93         elif DOMslices[k][1]-DOMslices[k][0]+1 < DOM_peak_frame:
  94             del DOMslices[k]
  95         else:
  96             k += 1
  97     while l < len(DISslices):
  98         if l+1 < len(DISslices) and DISslices[l+1][0]-DISslices[l][1] < DIS_join_frame:
  99             DISslices[l] = [ DISslices[l][0], DISslices[l+1][1] ]
 100             del DISslices[l+1]
 101         elif DISslices[l][1]-DISslices[l][0]+1 < DIS_peak_frame:
 102             del DISslices[l]
 103         else:
 104             l += 1
 105     return DOMslices, DISslices
 106
 107
 108 def SavitzkyGolay(window,derivative,datalist):
 109     SG_bin = relpath +'/'+ 'sav_gol -V0 '
 110     stdin, stdout, stderr = popen3(SG_bin + '-D' + str(derivative) + ' -n' + str(window)+','+str(window))
 111     for data in datalist:
 112         stdin.write(`data`+'\n')
 113     try:
 114         stdin.close()
 115     except:
 116         print stderr.readlines()
 117     results = stdout.readlines()
 118     stdout.close()
 119     SG_results = []
 120     for result in results:
 121         SG_results.append(float(fpformat.fix(result,6)))
 122     return SG_results
 123
 124
 125 def reportSlicesTXT(slices, sequence, maskFlag):
 126     if maskFlag == 'DOM':
 127         coordstr = 'GlobDoms '
 128     elif maskFlag == 'DIS':
 129         coordstr = 'Disorder '
 130     else:
 131         raise SystemExit
 132     if slices == []:
 133         #by default the sequence is in uppercase which is our search space
 134         s = sequence
 135     else:
 136         # insert seq before first slide
 137         if slices[0][0] > 0:
 138             s = sequence[0:slices[0][0]]
 139         else:
 140             s = ''
 141         for i in range(len(slices)):
 142             #skip first slice
 143             if i > 0:
 144                 coordstr = coordstr + ', '
 145             coordstr = coordstr + str(slices[i][0]+1) + '-' + str(slices[i][1]+1)
 146             #insert the actual slice
 147             if maskFlag == 'DOM':
 148                 s = s + lower(sequence[slices[i][0]:(slices[i][1]+1)])
 149                 if i < len(slices)-1:
 150                     s = s + upper(sequence[(slices[i][1]+1):(slices[i+1][0])])
 151                 #last slice
 152                 elif slices[i][1] < len(sequence)-1:
 153                     s = s + lower(sequence[(slices[i][1]+1):(len(sequence))])
 154             elif maskFlag == 'DIS':
 155                 s = s + upper(sequence[slices[i][0]:(slices[i][1]+1)])
 156                 #insert untouched seq between disorder segments, 2-run labelling
 157                 if i < len(slices)-1:
 158                     s = s + sequence[(slices[i][1]+1):(slices[i+1][0])]
 159                 #last slice
 160                 elif slices[i][1] < len(sequence)-1:
 161                     s = s + sequence[(slices[i][1]+1):(len(sequence))]
 162     return s,coordstr
 163
 164
 165 def runGlobPlot():
 166     try:
 167         smoothFrame = 10
 168         DOM_joinFrame = 15
 169         DOM_peakFrame = 74
 170         DIS_joinFrame = 4
 171         DIS_peakFrame = 5
 172         file = str(sys.argv[1])
 173         db = open(file,'r')
 174     except:
 175         print 'Usage:'
 176         print '         ./GlobPipe.py FASTAfile'
 177         raise SystemExit
 178     parser = Fasta.RecordParser()
 179     iterator = Fasta.Iterator(db,parser)
 180     while 1:
 181         try:
 182             cur_record = iterator.next()
 183             #uppercase is searchspace
 184             seq = upper(cur_record.sequence)
 185             # sum function
 186             sum_vector = Sum(seq,RL)
 187             # Run Savitzky-Golay
 188             smooth = SavitzkyGolay(`smoothFrame`,0, sum_vector)
 189             dydx_vector = SavitzkyGolay(`smoothFrame`,1, sum_vector)
 190             #test
 191             sumHEAD = sum_vector[:smoothFrame]
 192             sumTAIL = sum_vector[len(sum_vector)-smoothFrame:]
 193             newHEAD = []
 194             newTAIL = []
 195             for i in range(len(sumHEAD)):
 196                 try:
 197                     dHEAD = (sumHEAD[i+1]-sumHEAD[i])/2
 198                 except:
 199                     dHEAD = (sumHEAD[i]-sumHEAD[i-1])/2
 200                 try:
 201                     dTAIL = (sumTAIL[i+1]-sumTAIL[i])/2
 202                 except:
 203                     dTAIL = (sumTAIL[i]-sumTAIL[i-1])/2
 204                 newHEAD.append(dHEAD)
 205                 newTAIL.append(dTAIL)
 206             dydx_vector[:smoothFrame] = newHEAD
 207             dydx_vector[len(dydx_vector)-smoothFrame:] = newTAIL
 208             globdoms, globdis = getSlices(dydx_vector, DOM_joinFrame, DOM_peakFrame, DIS_joinFrame, DIS_peakFrame)
 209             s_domMask, coordstrDOM = reportSlicesTXT(globdoms, seq, 'DOM')
 210             s_final, coordstrDIS = reportSlicesTXT(globdis, s_domMask, 'DIS')
 211             sys.stdout.write('>'+cur_record.title+'\n')
 212             sys.stdout.write('# '+coordstrDOM+'\n')
 213             sys.stdout.write('# '+coordstrDIS+'\n')
 214
 215 # UNCOMMENT THIS IF NEED TO PRODUCE PER RESEDUE VALUES
 216             sys.stdout.write('# RESIDUE' + '\t' + 'DYDX' + '\t' + 'RAW' + '\t' +'SMOOTHED\n')
 217             for i in range(len(dydx_vector)):
 218 # dydx (positive values seems to indicate disorder in rows more than ~6 chars)  raw    smoothed
 219                sys.stdout.write(seq[i]+'\t'+fpformat.fix(dydx_vector[i],4)+ '\t'+fpformat.fix(smooth[i],4)+'\t'+fpformat.fix(sum_vector[i],4)+ '\n')
 220
 221 #            print s_final
 222             print '\n'
 223         except AttributeError:
 224             break
 225     return
 226
 227 runGlobPlot()