--- /dev/null
+#!/usr/bin/env python
+# Copyright (C) 2003 Rune Linding - EMBL
+# GlobPlot TM
+# GlobPlot is licensed under the Academic Free license
+
+from string import *
+from sys import argv
+import sys,re
+import os
+from os import system,popen3
+
+relpath = re.sub("/GlobPlot.py$","",argv[0])
+# cwd = re.sub("/$","", os.getcwd())
+newpath = relpath + "/biopython-1.50"
+sys.path.append(newpath)
+
+import Bio
+from Bio import File
+from Bio import Fasta
+import fpformat
+import tempfile
+import math
+
+# Russell/Linding
+RL = {'N':0.229885057471264,'P':0.552316012226663,'Q':-0.187676577424997,'A':-0.261538461538462,'R':-0.176592654077609, \
+ 'S':0.142883029808825,'C':-0.0151515151515152,'T':0.00887797506611258,'D':0.227629796839729,'E':-0.204684629516228, \
+ 'V':-0.386174834235195,'F':-0.225572305974316,'W':-0.243375458622095,'G':0.433225711769886,'H':-0.00121743364986608, \
+ 'Y':-0.20750516775322,'I':-0.422234699606962,'K':-0.100092289621613,'L':-0.337933495925287,'M':-0.225903614457831}
+
+def Sum(seq,par_dict):
+ sum = 0
+ results = []
+ raws = []
+ sums = []
+ p = 1
+ for residue in seq:
+ try:
+ parameter = par_dict[residue]
+ except:
+ parameter = 0
+ if p == 1:
+ sum = parameter
+ else:
+ sum = sum + parameter#*math.log10(p)
+ ssum = float(fpformat.fix(sum,10))
+ sums.append(ssum)
+ p +=1
+ return sums
+
+def getSlices(dydx_data, DOM_join_frame, DOM_peak_frame, DIS_join_frame, DIS_peak_frame):
+ DOMslices = []
+ DISslices = []
+ in_DOMslice = 0
+ in_DISslice = 0
+ beginDOMslice = 0
+ endDOMslice = 0
+ beginDISslice = 0
+ endDISslice = 0
+ for i in range( len(dydx_data) ):
+ #close dom slice
+ if in_DOMslice and dydx_data[i] > 0:
+ DOMslices.append([beginDOMslice, endDOMslice])
+ in_DOMslice = 0
+ #close dis slice
+ elif in_DISslice and dydx_data[i] < 0:
+ DISslices.append([beginDISslice, endDISslice])
+ in_DISslice = 0
+ # elseif inSlice expandslice
+ elif in_DOMslice:
+ endDOMslice += 1
+ elif in_DISslice:
+ endDISslice += 1
+ # if not in slice and dydx !== 0 start slice
+ if dydx_data[i] > 0 and not in_DISslice:
+ beginDISslice = i
+ endDISslice = i
+ in_DISslice = 1
+ elif dydx_data[i] < 0 and not in_DOMslice:
+ beginDOMslice = i
+ endDOMslice = i
+ in_DOMslice = 1
+ #last slice
+ if in_DOMslice:
+ DOMslices.append([beginDOMslice, endDOMslice])
+ if in_DISslice:
+ DISslices.append([beginDISslice,endDISslice])
+ k = 0
+ l = 0
+ while k < len(DOMslices):
+ if k+1 < len(DOMslices) and DOMslices[k+1][0]-DOMslices[k][1] < DOM_join_frame:
+ DOMslices[k] = [ DOMslices[k][0], DOMslices[k+1][1] ]
+ del DOMslices[k+1]
+ elif DOMslices[k][1]-DOMslices[k][0]+1 < DOM_peak_frame:
+ del DOMslices[k]
+ else:
+ k += 1
+ while l < len(DISslices):
+ if l+1 < len(DISslices) and DISslices[l+1][0]-DISslices[l][1] < DIS_join_frame:
+ DISslices[l] = [ DISslices[l][0], DISslices[l+1][1] ]
+ del DISslices[l+1]
+ elif DISslices[l][1]-DISslices[l][0]+1 < DIS_peak_frame:
+ del DISslices[l]
+ else:
+ l += 1
+ return DOMslices, DISslices
+
+
+def SavitzkyGolay(window,derivative,datalist):
+ SG_bin = relpath +'/'+ 'sav_gol -V0 '
+ stdin, stdout, stderr = popen3(SG_bin + '-D' + str(derivative) + ' -n' + str(window)+','+str(window))
+ for data in datalist:
+ stdin.write(`data`+'\n')
+ try:
+ stdin.close()
+ except:
+ print stderr.readlines()
+ results = stdout.readlines()
+ stdout.close()
+ SG_results = []
+ for result in results:
+ SG_results.append(float(fpformat.fix(result,6)))
+ return SG_results
+
+
+def reportSlicesTXT(slices, sequence, maskFlag):
+ if maskFlag == 'DOM':
+ coordstr = 'GlobDoms '
+ elif maskFlag == 'DIS':
+ coordstr = 'Disorder '
+ else:
+ raise SystemExit
+ if slices == []:
+ #by default the sequence is in uppercase which is our search space
+ s = sequence
+ else:
+ # insert seq before first slide
+ if slices[0][0] > 0:
+ s = sequence[0:slices[0][0]]
+ else:
+ s = ''
+ for i in range(len(slices)):
+ #skip first slice
+ if i > 0:
+ coordstr = coordstr + ', '
+ coordstr = coordstr + str(slices[i][0]+1) + '-' + str(slices[i][1]+1)
+ #insert the actual slice
+ if maskFlag == 'DOM':
+ s = s + lower(sequence[slices[i][0]:(slices[i][1]+1)])
+ if i < len(slices)-1:
+ s = s + upper(sequence[(slices[i][1]+1):(slices[i+1][0])])
+ #last slice
+ elif slices[i][1] < len(sequence)-1:
+ s = s + lower(sequence[(slices[i][1]+1):(len(sequence))])
+ elif maskFlag == 'DIS':
+ s = s + upper(sequence[slices[i][0]:(slices[i][1]+1)])
+ #insert untouched seq between disorder segments, 2-run labelling
+ if i < len(slices)-1:
+ s = s + sequence[(slices[i][1]+1):(slices[i+1][0])]
+ #last slice
+ elif slices[i][1] < len(sequence)-1:
+ s = s + sequence[(slices[i][1]+1):(len(sequence))]
+ return s,coordstr
+
+
+def runGlobPlot():
+ try:
+ smoothFrame = 10
+ DOM_joinFrame = 15
+ DOM_peakFrame = 74
+ DIS_joinFrame = 4
+ DIS_peakFrame = 5
+ file = str(sys.argv[1])
+ db = open(file,'r')
+ except:
+ print 'Usage:'
+ print ' ./GlobPipe.py FASTAfile'
+ raise SystemExit
+ parser = Fasta.RecordParser()
+ iterator = Fasta.Iterator(db,parser)
+ while 1:
+ try:
+ cur_record = iterator.next()
+ #uppercase is searchspace
+ seq = upper(cur_record.sequence)
+ # sum function
+ sum_vector = Sum(seq,RL)
+ # Run Savitzky-Golay
+ smooth = SavitzkyGolay(`smoothFrame`,0, sum_vector)
+ dydx_vector = SavitzkyGolay(`smoothFrame`,1, sum_vector)
+ #test
+ sumHEAD = sum_vector[:smoothFrame]
+ sumTAIL = sum_vector[len(sum_vector)-smoothFrame:]
+ newHEAD = []
+ newTAIL = []
+ for i in range(len(sumHEAD)):
+ try:
+ dHEAD = (sumHEAD[i+1]-sumHEAD[i])/2
+ except:
+ dHEAD = (sumHEAD[i]-sumHEAD[i-1])/2
+ try:
+ dTAIL = (sumTAIL[i+1]-sumTAIL[i])/2
+ except:
+ dTAIL = (sumTAIL[i]-sumTAIL[i-1])/2
+ newHEAD.append(dHEAD)
+ newTAIL.append(dTAIL)
+ dydx_vector[:smoothFrame] = newHEAD
+ dydx_vector[len(dydx_vector)-smoothFrame:] = newTAIL
+ globdoms, globdis = getSlices(dydx_vector, DOM_joinFrame, DOM_peakFrame, DIS_joinFrame, DIS_peakFrame)
+ s_domMask, coordstrDOM = reportSlicesTXT(globdoms, seq, 'DOM')
+ s_final, coordstrDIS = reportSlicesTXT(globdis, s_domMask, 'DIS')
+ sys.stdout.write('>'+cur_record.title+'\n')
+ sys.stdout.write('# '+coordstrDOM+'\n')
+ sys.stdout.write('# '+coordstrDIS+'\n')
+
+# UNCOMMENT THIS IF NEED TO PRODUCE PER RESEDUE VALUES
+ sys.stdout.write('# RESIDUE' + '\t' + 'DYDX' + '\t' + 'RAW' + '\t' +'SMOOTHED\n')
+ for i in range(len(dydx_vector)):
+# dydx (positive values seems to indicate disorder in rows more than ~6 chars) raw smoothed
+ 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')
+
+# print s_final
+ print '\n'
+ except AttributeError:
+ break
+ return
+
+runGlobPlot()
git://source.jalview.org / jabaws.git / blobdiff