/
Fit_EC50.py
executable file
·174 lines (101 loc) · 3.85 KB
/
Fit_EC50.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
#!/usr/bin/python
import sys
import os
import os.path
import commands
import time
import math
import copy
import pprint
import collections
import numpy as np
from collections import OrderedDict
from scipy.optimize import curve_fit
import scipy
from scipy import stats
from scipy.stats import norm
from numpy import random
#from numpy.random import Generator
try:
infilename1 = sys.argv[1]
# outfilename = sys.argv[3]
inital_guess = [float(sys.argv[2]), float(sys.argv[3]) ]
# ns_per_step = float (sys.argv[6])
# try:
# ratio = float (sys.argv[7])
# except:
# ratio = 1
except:
print "Usage:",sys.argv[0], "infile ns per step outfile"; sys.exit(1)
#########################################
print(' ')
print(' ')
print(' ')
print(' ')
print(' ')
ifile1 = open(infilename1,'r') # open index pdb, we need resid
ofile = open(str(infilename1+'.fit'),'w') # open file for writing
def func_hill(L, EC50 , Gmax):
return Gmax*L / (L + EC50 )
counter = 0
data = OrderedDict()
L_index = []
GV_shift_index_current = []
SEM_index = []
for line in ifile1:
columns = line.split()
if len(columns) == 3 : # and line[0:4] == 'ATOM':
L_current = float(columns[0])
GV_shift_current = float(columns[1])
SEM_current = float(columns[2])
L_index.append(L_current)
GV_shift_index_current.append(GV_shift_current)
SEM_index.append(SEM_current)
data[L_current] = OrderedDict()
data[L_current]['GV_shift'] = GV_shift_current
data[L_current]['SEM'] = SEM_current
print 'if not add SEM'
popt, pcov = curve_fit(func_hill, L_index, GV_shift_index_current)
#popt array
# Optimal values for the parameters so that the sum of the squared residuals of f
#pcov a2-D array
# The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters use perr = np.sqrt(np.diag(pcov)).
#
print 'EC50, Gmax, fit with raw mean ', popt
print 'if add SEM'
popt, pcov = curve_fit(func_hill, L_index, GV_shift_index_current, inital_guess, SEM_index, True)
print 'EC50, Gmax, fit with raw mean ', popt
#print 'estimated covariance of popt', pcov
perr = np.sqrt(np.diag(pcov))
print perr
##########
print ' now generate a data set with same mean on each L'
data_reconstruction = OrderedDict()
for L_current in L_index:
GV_shift_current = data[L_current]['GV_shift']
SEM_current = data[L_current]['SEM']
num_samples = 1000
data_reconstruction[L_current] = OrderedDict()
data_reconstruction[L_current]['GV_shift'] = np.random.normal(GV_shift_current, SEM_current, num_samples)
# data_reconstruction[L_current]['GV_shift'] = np.random.Generator.normal(GV_shift_current, SEM_current, num_samples) # we do not have thie model in current version
print L_current, np.mean(data_reconstruction[L_current]['GV_shift']), np.std(data_reconstruction[L_current]['GV_shift'])
########### now
print ' not fit from random choice '
# each time GV_shift_index_current is regenerate
num_samples = 10000
current_samples = 0
Gmax_list = []
EC50_list = []
while current_samples < num_samples:
GV_shift_index_current = []
for L_current in L_index:
GV_shift_current = random.choice( data_reconstruction[L_current]['GV_shift'] )
GV_shift_index_current.append(GV_shift_current)
popt, pcov = curve_fit(func_hill, L_index, GV_shift_index_current)
Gmax_list.append(popt[1])
EC50_list.append(popt[0])
if float(current_samples)%1000 == 0 :
print ' done frame ', current_samples
current_samples += 1
print 'EC50 sampled', np.mean(EC50_list), np.std(EC50_list)
print 'Gmax sampled', np.mean(Gmax_list), np.std(Gmax_list)