forked from Pymol-Scripts/Pymol-script-repo
-
Notifications
You must be signed in to change notification settings - Fork 0
/
forster_distance_calculator.py
319 lines (287 loc) · 17.4 KB
/
forster_distance_calculator.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
'''
Described at PyMOL wiki:
http://www.pymolwiki.org/index.php/forster_distance_calculator
#-------------------------------------------------------------------------------
# Name: Forster
# Purpose: Forster resonance energy transfer calculator.
# Input is manufactor provided spectres of Donor emission and
# acceptor excitation spectrum.
#
# Carl Boswell and co. have made a new homepage with a long list of dyes which can be downloaded.
# With a graphics program, they have traced several spectre of dyes from the literature and made this easily public at:
# http://www.spectra.arizona.edu/ I highly recommend this homepage.
# With these Spectra, the script can calculate the Forster Distance for different dyes from different companies.
# Download "one spectrum at the time" by "deselecting" one of the spectre in the right side of the graph window.
# Then get the datafile with the button in the lower right corner.
##
# Made from
# http://en.wikipedia.org/wiki/F%C3%B6rster_resonance_energy_transfer#Theoretical_basis
# {R_0}^6 = \frac{9\,Q_0 \,(\ln 10) \kappa^2 \, J}{128 \, \pi^5 \,n^4 \, N_A}
#
# Author: Troels Emtekaer Linnet
#
# Created: 29/03/2011
# Copyright: (c) tlinnet 2011
# Licence: Free for all
-------------------------------------------------------------------------------
#Ref(1)
# Biochemistry 1997, 36, 11261-11272
# M. Pilar Lillo, Joseph M. Beechem, Barbara K. Szpikowska, Mark A. Sherman, and Maria T. Mas
#Design and Characterization of a Multisite Fluorescence Energy-Transfer System for Protein Folding Studies: A Steady-State and Time-Resolved Study of Yeast Phosphoglycerate Kinase
NOTES:
Datafiles: Two column file. Space separated. Numbers are "." dot noted. First column is in nanometers nm. Second column is arbitrary units of fluorescence/emission.
If you have collected Donor Exitation and Acceptor Emission, they can be collected and plotted in gnuplot automatically. Set: Compare"yes"
xunit="nm": Enter x values in "nm" or "cm".
A_e_Max_Y : Acceptor maximum molar extinction coefficient. In units of M-1 cm-1. Approx 60000 - 1500000
A_e_Max_X : Enter at which wavelength (in nm) the maximum absorption occurs.
Qd=0.8 : Fluorescence quantum yield of the donor in the absence of the acceptor. Qd = neta_fl = n_fl / n_abs = n_emi / n_exi
k2 = 2.0/3.0 Dipole orientation factor.
n = 1.33 : Refractive index of the medium. water=1.33, protein=1.4, n2MGuHCl=1.375 Ref(1)
NA = 6.02214179e+023 # (units: Number*mol-1 )Avogadros number
'''
from __future__ import print_function
try:
from pymol import cmd
runningpymol = 'yes'
except:
runningpymol = 'no'
pass
import os
import platform
import math
def forster(D_Exi="ATTO488Exi.txt", D_Emi="ATTO488Emi.txt", A_Exi="ATTO590Exi.txt", A_Emi="ATTO590Emi.txt", A_e_Max_Y=120000, A_e_Max_X=594, Qd=0.8, k2=0.66667, n=1.33, Compare="yes", xunit="nm"):
A_e_Max_Y = float(A_e_Max_Y)
A_e_Max_X = float(A_e_Max_X)
Qd = float(Qd)
k2 = float(k2)
n = float(n)
NA = 6.02214179e+023
print(k2, Qd)
printAll = "ye" # To print out all info
fileDexi, extDexi = os.path.splitext(D_Exi)
fileDemi, extDemi = os.path.splitext(D_Emi)
fileAexi, extAexi = os.path.splitext(A_Exi)
fileAemi, extAemi = os.path.splitext(A_Emi)
overlapname = fileDemi + "-" + fileAexi + "-overlap.dat"
overlapfile = open(overlapname, "w")
overlapgnuplotname = fileDemi + "-" + fileAexi + "-overlap.plt"
overlapgnuplotfile = open(overlapgnuplotname, "w")
print("\nI have opened two files for you: \n%s and %s" % (overlapname, overlapgnuplotname))
print("The .plt should be opened with gnuplot to make the graphs.")
print("The created graphs are .eps files.")
print("They can be converted to pdf with the program: epstopdf or eps2pdf")
print('Part of LaTeX: C:\Program Files (x86)\MiKTeX 2.9\miktex' + "\\" + "bin")
print("Or download here: http://tinyurl.com/eps2pdf")
DonorEmi = open(D_Emi, "r")
AcceptorExi = open(A_Exi, "r")
lineDemi = DonorEmi.readlines()
lineAexi = AcceptorExi.readlines()
Demi = []
Aexi = []
for i in lineDemi:
if not i.strip(): # If line cannot get stripped(does not exist), then continue
continue
else: # If line can get stripped
if testfloat(str.split(i)[0]):
Demi.append([float(str.split(i)[0]), float(str.split(i)[1])])
AreaDemi = numintegrator(Demi)
print("Nummerical integration of Donor emission spectrum, used for normalization, gives: Area=", AreaDemi)
for i in lineAexi:
if not i.strip():
continue
else:
if testfloat(str.split(i)[0]):
Aexi.append([float(str.split(i)[0]), float(str.split(i)[1])])
if float(str.split(i)[0]) == float(A_e_Max_X):
Epsiloncorrection = [float(A_e_Max_X), float(str.split(i)[0]), float(str.split(i)[1])]
# Making the overlap
OverlapDataPoints = []
OverlapSum = 0.0
# For comparing two floating numbers, one have to be carefully. Setting error allowing difference
eallow = 0.00000001
for i in range(len(Demi)):
for j in range(len(Aexi)):
if Demi[i][0] - eallow < Aexi[j][0] and Demi[i][0] + eallow > Aexi[j][0]:
Overlap = (Demi[i][1] * Aexi[j][1] * float(A_e_Max_Y) * math.pow(Demi[i][0], 4)) / (AreaDemi * Epsiloncorrection[2])
OverlapSum = OverlapSum + Overlap
OverlapDataPoints.append([Demi[i][0], Demi[i][1], Aexi[j][0], Aexi[j][1], Overlap, OverlapSum])
AreaOverlap = numintegrator(OverlapDataPoints, 0, 4)
Prefactor = ForsterPrefactor6(Qd, k2, n, NA, printAll)
ForsterAng = ForsterCalc(Prefactor, AreaOverlap, xunit, printAll)
# Outputting data
overlapfile.write("Emi-wavelength Emi-value-norm1 Emi-value-normA Exi-wavelength Exi-value-norm1 Exti-coefficient Overlap Overlap-Sum\n")
for line in range(len(OverlapDataPoints)):
textline = "%4.1f %24.4f %15.4e %14.1f %15.4e %16.4e %12.4e %13.4e" % (OverlapDataPoints[line][0], OverlapDataPoints[line][1], float(OverlapDataPoints[line][1] / AreaDemi), OverlapDataPoints[line][2], OverlapDataPoints[line][3], float(A_e_Max_Y * OverlapDataPoints[line][3] / Epsiloncorrection[2]), float(OverlapDataPoints[line][4]), float(OverlapDataPoints[line][5]))
overlapfile.write(textline + "\n")
# Make gnuplot plot file
overlapgnuplotfile.write("reset" + "\n")
overlapgnuplotfile.write("cd " + "'" + os.getcwd() + "'" + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("set xrange [400:800]" + "\n")
overlapgnuplotfile.write("set ytics nomirror" + "\n")
overlapgnuplotfile.write("set y2tics" + "\n")
if xunit == "cm":
overlapgnuplotfile.write("set xlabel 'Wavelength (cm)'" + "\n")
else:
overlapgnuplotfile.write("set xlabel 'Wavelength (nm)'" + "\n")
overlapgnuplotfile.write("set size ratio 0.5" + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("A_e_Max_Y = " + str(A_e_Max_Y) + "\n")
overlapgnuplotfile.write("A_e_Max_X = " + str(A_e_Max_X) + "\n")
overlapgnuplotfile.write("AreaDemi = " + str(AreaDemi) + "\n")
overlapgnuplotfile.write("AreaOverlap = " + str(AreaOverlap) + "\n")
overlapgnuplotfile.write("ForsterAng= " + str(ForsterAng) + "\n")
overlapgnuplotfile.write("\n")
if Compare == "yes":
overlapgnuplotfile.write("#########################Graph 1#############################" + "\n")
overlapgnuplotfile.write('set title ' + '"' + fileDemi + "-" + fileAexi + ' FRET Donor/Acceptor spectre"' + "\n")
overlapgnuplotfile.write('set ylabel "Donor Fluorescence Intensity F_{D}({/Symbol l}) \\n Acceptor Extinction coefficient {/Symbol e}({/Symbol l})"' + "\n")
overlapgnuplotfile.write('set y2label "F_{D}({/Symbol l}){/Symbol e}({/Symbol l})^{norm1}{/Symbol l}^{4}"' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write('set label 1 "Donor Emission Area= %g", AreaDemi at graph 0.7, -0.15' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write('set term postscript eps enhanced color' + "\n")
overlapgnuplotfile.write('set output ' + '"1-' + fileDemi + "-" + fileAexi + '-overlap-all-spectre.eps"' + "\n")
overlapgnuplotfile.write('plot ' + '"' + fileDexi + extDexi + '" using 1:2 title ' + '"' + fileDexi + ' exitation" with lines,\\' + "\n")
# overlapgnuplotfile.write('"'+overlapname+'" using 1:2 title '+'"'+fileDemi+' emission" with lines,\\'+"\n")
# overlapgnuplotfile.write('"'+overlapname+'" using 4:5 title '+'"'+fileAexi+' exitation" with lines,\\'+"\n")
overlapgnuplotfile.write('"' + fileDemi + extDemi + '" using 1:2 title ' + '"' + fileDemi + ' emission" with lines,\\' + "\n")
overlapgnuplotfile.write('"' + fileAexi + extAexi + '" using 1:2 title ' + '"' + fileAexi + ' exitation" with lines,\\' + "\n")
overlapgnuplotfile.write('"' + fileAemi + extAemi + '" using 1:2 title ' + '"' + fileAemi + ' emission" with lines,\\' + "\n")
overlapgnuplotfile.write('"' + overlapname + '" using 1:($2*$5*$1**4) title "D/A Overlap :y2" with lines axis x1y2' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("## Show in window, x11 for Linux" + "\n")
overlapgnuplotfile.write("#set term x11" + "\n")
overlapgnuplotfile.write("#set term windows" + "\n")
overlapgnuplotfile.write("#replot" + "\n")
overlapgnuplotfile.write("#pause -1" + "\n")
overlapgnuplotfile.write("unset label" + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("#########################Graph 2#############################" + "\n")
overlapgnuplotfile.write('set title ' + '"' + fileDemi + "-" + fileAexi + ' FRET Donor/Acceptor overlap"' + "\n")
overlapgnuplotfile.write('set ylabel "Donor Fluorescence Intensity \\n Normalized by F_{D}({/Symbol l})^{normA}=F_{D}({/Symbol l})/F_{Area}"' + "\n")
overlapgnuplotfile.write('set y2label "Acceptor Extinction coefficient [M^{-1}cm^{-1}] \\n Normalized to {/Symbol e}({/Symbol l})"' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write('set label 1 "{/Symbol e}(max)= %g", A_e_Max_Y at graph 0.63, 0.65' + "\n")
overlapgnuplotfile.write('set label 2 "at %g ' + xunit + '", A_e_Max_X at graph 0.63, 0.60' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write('set term postscript eps enhanced color' + "\n")
overlapgnuplotfile.write('set output ' + '"2-' + fileDemi + "-" + fileAexi + '-overlap-normalized-spectre.eps"' + "\n")
overlapgnuplotfile.write('plot ' + '"' + overlapname + '" using 1:3 title ' + '"' + fileDemi + ', Normalized by area, emission" with lines,\\' + "\n")
overlapgnuplotfile.write('"' + overlapname + '" using 4:6 title ' + '"' + fileAexi + ' Extinction coefficient :y2" with lines axis x1y2' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("## Show in window, x11 for Linux" + "\n")
overlapgnuplotfile.write("#set term x11" + "\n")
overlapgnuplotfile.write("#set term windows" + "\n")
overlapgnuplotfile.write("#replot" + "\n")
overlapgnuplotfile.write("#pause -1" + "\n")
overlapgnuplotfile.write("unset label" + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("#########################Graph 3#############################" + "\n")
overlapgnuplotfile.write('set title ' + '"' + fileDemi + "-" + fileAexi + ' FRET Donor/Acceptor overlap integral"' + "\n")
if xunit == "cm":
overlapgnuplotfile.write('set ylabel "Donor/Acceptor overlap [cm^{3}M^{-1}] \\n F_{D}({/Symbol l})^{normA}{/Symbol e}({/Symbol l}){/Symbol l}^{4}"' + "\n")
overlapgnuplotfile.write('set y2label "Donor/Acceptor overlap integral [cm^{3}M^{-1}] \\n {/Symbol S} F_{D}({/Symbol l})^{normA}{/Symbol e}({/Symbol l}){/Symbol l}^{4}"' + "\n")
else:
overlapgnuplotfile.write('set ylabel "Donor/Acceptor overlap [cm^{-1}nm^{4}M^{-1}] \\n F_{D}({/Symbol l})^{normA}{/Symbol e}({/Symbol l}){/Symbol l}^{4}"' + "\n")
overlapgnuplotfile.write('set y2label "Donor/Acceptor overlap integral [cm^{-1}nm^{4}M^{-1}] \\n {/Symbol S} F_{D}({/Symbol l})^{normA}{/Symbol e}({/Symbol l}){/Symbol l}^{4}"' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write('set label 1 "Overlap integral:" at graph 0.55, 0.65' + "\n")
overlapgnuplotfile.write('set label 2 "{/Symbol S}= %g", AreaOverlap at graph 0.55, 0.60' + "\n")
overlapgnuplotfile.write('set label 3 "Forster Distance:" at graph 0.55, 0.50' + "\n")
overlapgnuplotfile.write('set label 5 "R_{0}= %g Angstrom", ForsterAng at graph 0.55, 0.45' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write('set term postscript eps enhanced color' + "\n")
overlapgnuplotfile.write('set output ' + '"3-' + fileDemi + "-" + fileAexi + '-overlap-integral.eps"' + "\n")
overlapgnuplotfile.write('plot ' + '"' + overlapname + '" using 1:7 title "Overlap: F_{D}({/Symbol l})^{normA}{/Symbol e}({/Symbol l}) {/Symbol l}^{4}" with lines,\\' + "\n")
overlapgnuplotfile.write('"' + overlapname + '" using 1:8 title "Integral: {/Symbol S} F_{D}({/Symbol l})^{normA} {/Symbol e}({/Symbol l}){/Symbol l}^{4} :y2" with lines axis x1y2' + "\n")
overlapgnuplotfile.write("\n")
overlapgnuplotfile.write("## Show in window, x11 for Linux" + "\n")
overlapgnuplotfile.write("#set term x11" + "\n")
overlapgnuplotfile.write("#set term windows" + "\n")
overlapgnuplotfile.write("#replot" + "\n")
overlapgnuplotfile.write("#pause -1" + "\n")
overlapgnuplotfile.write("unset label")
overlapgnuplotfile.close()
overlapfile.close()
return(ForsterAng)
if runningpymol != 'no':
cmd.extend("forster", forster)
def ForsterConstFactor6(NA, printAll):
vForsterConstFactor6 = (9 * math.log(10)) / (128 * math.pow(math.pi, 5) * NA)
if printAll == 'yes':
print("Forster constant pre-factor is:", str(vForsterConstFactor6), "(units: mol)")
return vForsterConstFactor6
def ForsterVariableFactor6(Qd, k2, n, printAll):
vForsterVariableFactor6 = (k2 * Qd) / (math.pow(n, 4))
if printAll == 'yes':
print("Forster variable pre-factor is:", str(vForsterVariableFactor6), "(units: NIL)")
return vForsterVariableFactor6
def ForsterPrefactor6(Qd, k2, n, NA, printAll):
vForsterPrefactor6 = ForsterConstFactor6(NA, printAll) * ForsterVariableFactor6(Qd, k2, n, printAll)
if printAll == 'yes':
print("Forster pre-factor is:", str(vForsterPrefactor6), "(units: mol)")
return vForsterPrefactor6
def ForsterCalcnm(fFPreFactor6, fAreaOverlap, printAll):
if printAll == 'yes':
print("Overlap sum is: ", str(fAreaOverlap), "(units: cm-1 nm^4 L mol-1)")
Forster6 = fFPreFactor6 * fAreaOverlap
if printAll == 'yes':
print("Forster distance 6th power:", str(Forster6), "(units: cm-1 nm^4 L), obs(1L=1e-3m^3)")
Forster6m = Forster6 * 100 * math.pow(1e-9, 4) * 1e-3 # 1e-3 is conversion from 1L = 1e-3 m3
if printAll == 'yes':
print("Forster distance 6th power:", str(Forster6m), "(units: meter m^6)")
Forster6Ang = Forster6m * math.pow(1e10, 6.0)
if printAll == 'yes':
print("Forster distance Angstrom 6th power:", "%e" % (Forster6Ang), "(units: Angstrom^6)")
ForsterAng = math.pow(Forster6Ang, 1.0 / 6.0)
print("Forster distance:", str(ForsterAng), "(units: Angstrom)")
return ForsterAng
def ForsterCalccm(fFPreFactor6, fAreaOverlap, printAll):
if printAll == 'yes':
print("Overlap sum is: ", str(fAreaOverlap), "(units: cm^3 L mol-1)")
Forster6 = fFPreFactor6 * fAreaOverlap
if printAll == 'yes':
print("Forster distance 6th power:", str(Forster6), "(units: cm^3 L), obs(1L=1e-3m^3)")
Forster6m = Forster6 * math.pow(1e-2, 3) * 1e-3 # 1e-3 is conversion from 1L = 1e-3 m3
if printAll == 'yes':
print("Forster distance 6th power:", str(Forster6m), "(units: meter m^6)")
Forster6cm = Forster6m * math.pow(1e2, 6.0)
if printAll == 'yes':
print("Forster distance cm 6th power:", "%e" % (Forster6cm), "(units: cm^6)")
Forster6Ang = Forster6m * math.pow(1e10, 6.0)
if printAll == 'yes':
print("Forster distance Angstrom 6th power:", "%e" % (Forster6Ang), "(units: Angstrom^6)")
ForsterAng = math.pow(Forster6m, 1.0 / 6.0)
print("Forster distance:", str(ForsterAng), "(units: Angstrom)")
return ForsterAng
def ForsterCalc(fFPreFactor6, fAreaOverlap, xunit, printAll):
if xunit == "nm":
Value = ForsterCalcnm(fFPreFactor6, fAreaOverlap, printAll)
if xunit == "cm":
Value = ForsterCalccm(fFPreFactor6, fAreaOverlap, printAll)
return Value
def testfloat(x):
try:
v = float(x)
return x
except:
return False
def numintegrator(fluarray, col1=0, col2=1):
xprev = 0
xpres = 0
yprev = 0
ypres = 0
summing = 0
for i in range(len(fluarray)):
# Have to skip first datapoint
if i > 0:
xprev = xpres
yprev = ypres
xpres = fluarray[i][col1]
ypres = fluarray[i][col2]
summing = yprev * (xpres - xprev) + (ypres - yprev) * (xpres - xprev) / 2.0 + summing
else:
xpres = fluarray[i][col1]
ypres = fluarray[i][col2]
return summing