1 |
""" |
2 |
ba_class |
3 |
-------- |
4 |
|
5 |
Biosensor Array class for storing SPRI data. |
6 |
This contains an array of SPR sensorgrams with their |
7 |
associated models and parameters. |
8 |
This can also contain the microarray geometery. |
9 |
|
10 |
.. moduleauthor:: Christopher Lausted, |
11 |
Institute for Systems Biology, |
12 |
Yuhang Wan, |
13 |
OSPRAI developers. |
14 |
|
15 |
Examples:: |
16 |
|
17 |
>>> ba1 = BiosensorArray(800,1500) ## Allocate object with 800 spots of 1500 time points. |
18 |
>>> ba1.roi[799].time[1499] = 7502.1 ## Set time in seconds. |
19 |
>>> ba1.roi[799].value[1499] = 0.50 ## Set SPR reflectance signal. |
20 |
>>> ba1.roi[799].name = "anti-IgG" ## Set name of microarray feature. |
21 |
>>> print ba1.roi[799].name |
22 |
anti-IgG |
23 |
>>> ba1 = BiosensorArray(1, 300) |
24 |
>>> ba1.roi[0].simulate(bulkRI=1000, kon=1e4, koff=1e-2, conc=1e-6) |
25 |
>>> print round(ba1.roi[0].value[299], 1) |
26 |
160.6 |
27 |
""" |
28 |
__version__ = "111012" |
29 |
|
30 |
|
31 |
## Import libraries. |
32 |
import numpy as np |
33 |
import cPickle as pickle |
34 |
from datetime import date |
35 |
|
36 |
|
37 |
class BiosensorArray(object): |
38 |
""" |
39 |
The Biosensor Array class holds all ROIs (spots) |
40 |
from an SPRI microarray run. |
41 |
|
42 |
================ =============== ======================================= |
43 |
Property Type Description |
44 |
---------------- --------------- --------------------------------------- |
45 |
roi list of object List of RegOfInterest objects. |
46 |
rois integer Size of roi (Number of ROIs). |
47 |
dpoints integer Number of datapoints in each ROI. |
48 |
primarydatafiles list of string List of file names. |
49 |
comments list of string List of comments regarding experiments. |
50 |
================ =============== ======================================= |
51 |
""" |
52 |
|
53 |
def __init__(self, roi_size, datapoint_size): |
54 |
"""This object is dimensioned upon initialization.""" ## Might improve speed. |
55 |
## Use a list comprehension to dimension list of classes. |
56 |
self.roi = [RegOfInterest(i,datapoint_size) for i in range(roi_size)] |
57 |
## Remember the dimensions for quick bounds-checking later. |
58 |
self.rois = roi_size |
59 |
self.dpoints = datapoint_size |
60 |
## Use these to keep track of laboratory notes. |
61 |
self.primarydatafiles = [""] ## Like "spritdata.txt" |
62 |
self.comments = [""] ## Like "Antibody data from 1 Jan 2010" |
63 |
return |
64 |
|
65 |
def xy2uid(self, gx, gy, x, y): |
66 |
"""Given roi coordinates, find corresponding uid (index).""" |
67 |
coor1 = (gx, gy, x, y) |
68 |
for roi in self.roi: |
69 |
coor2 = (roi.gridx, roi.gridy, roi.spotx, roi.spoty) |
70 |
if (coor1==coor2): return roi.uid ## Success. |
71 |
return -1 ## Failure. |
72 |
|
73 |
def set_plot_all(self): |
74 |
"""Choose to plot every sensorgram.""" |
75 |
for roi in self.roi: roi.plottable = True |
76 |
return |
77 |
|
78 |
def set_plot_list(self, ilist): |
79 |
"""Choose a list of sensorgrams to plot""" |
80 |
for roi in self.roi: roi.plottable = False |
81 |
for i in ilist: self.roi[i].plottable = True |
82 |
return |
83 |
|
84 |
def save(self, filename="default.sba"): |
85 |
""" |
86 |
Save (serialize) this ba object as a file. |
87 |
Use the extension '.sba' and designate it file format 1. |
88 |
Use cPickle protocol 2 (binary). |
89 |
""" |
90 |
if (filename[-4:] != '.sba'): filename = filename + ".sba" |
91 |
fp = open(filename, 'wb') |
92 |
protocol = 2 ## Use cPickle protocol 2 (binary). |
93 |
formatver = 1 ## My format version 1. |
94 |
pickle.dump(formatver, fp, protocol) |
95 |
pickle.dump(self.comments, fp, protocol) |
96 |
pickle.dump(self.primarydatafiles, fp, protocol) |
97 |
pickle.dump(self.dpoints, fp, protocol) |
98 |
pickle.dump(self.rois, fp, protocol) |
99 |
pickle.dump(self.roi, fp, protocol) |
100 |
|
101 |
fp.close() |
102 |
return |
103 |
""" End of save() method""" |
104 |
|
105 |
def load(self, filename="default.sba"): |
106 |
""" |
107 |
Load data from sba file into this ba object. |
108 |
""" |
109 |
fp = open(filename, 'rb') |
110 |
formatver = pickle.load(fp) |
111 |
if (formatver == 1): |
112 |
self.comments = pickle.load(fp) |
113 |
self.primarydatafiles = pickle.load(fp) |
114 |
self.dpoints = pickle.load(fp) |
115 |
self.rois = pickle.load(fp) |
116 |
self.roi = pickle.load(fp) |
117 |
else: |
118 |
print "Error: file format version is not 1 but " + formatver |
119 |
fp.close() |
120 |
return |
121 |
""" End of load() method""" |
122 |
|
123 |
"""End of BiosensorArray class definition.""" |
124 |
|
125 |
|
126 |
class RegOfInterest(object): |
127 |
""" |
128 |
This Region Of Interest class can hold one SPR sensorgram. |
129 |
|
130 |
=========== =============== ============================================= |
131 |
Property Type Description |
132 |
----------- --------------- --------------------------------------------- |
133 |
uid integer Unique ID integer. Never changes. |
134 |
index integer Changes if ba size changes (add/remove rois). |
135 |
time nparray Usually seconds. |
136 |
value nparray Arbitrary units. |
137 |
dpoints integer Can use for bounds checking. |
138 |
name string Name of ROI. |
139 |
desc string Description of ROI |
140 |
|
141 |
gridx integer Can be Block number from GAL file. |
142 |
gridy integer Unused? |
143 |
spotx integer ROI column number. |
144 |
spoty integer ROI row number. |
145 |
bgroi list of integer One or more background ROIs. |
146 |
calibM float Slope used for calibration |
147 |
calibB float Intercept used for calibration. |
148 |
plottable boolean Whether to plot |
149 |
|
150 |
injconc list of float One or more analyte concentrations. |
151 |
injstart list of float One or more injection start times. |
152 |
injstop list of float One or more injection start times. |
153 |
injrind list of float One or more expected refractive index jumps. |
154 |
|
155 |
flow float Flowrate of injection. |
156 |
model ref to funct Model describing this ROI. |
157 |
params dict of dict Parameters for this model. |
158 |
=========== =============== ============================================= |
159 |
""" |
160 |
|
161 |
def __init__(self, i, datapoint_size): |
162 |
""" |
163 |
Each ROI in a list needs a unique ID number which will be i+1. |
164 |
The time and value arrays must be allocated using datapoint_size. |
165 |
""" |
166 |
self.uid = i+1 ## Unique ID integer. Never changes. |
167 |
self.index = i ## Changes if ba size changes (add/remove rois). |
168 |
self.time = np.zeros(datapoint_size, dtype=float) ## Usually seconds. |
169 |
self.value = np.zeros(datapoint_size, dtype=float) ## Arbitrary units. |
170 |
self.dpoints = datapoint_size ## Can use for bounds checking. |
171 |
## Useful optional information. Usually from GAL or Key files. |
172 |
self.name = "Spot%i" % (i+1) |
173 |
self.desc = "Description" |
174 |
self.gridx = 0 ## Can be Block number from GAL file. |
175 |
self.gridy = 0 ## Unused? |
176 |
self.spotx = 0 ## ROI column number. |
177 |
self.spoty = 0 ## ROI row number. |
178 |
self.bgroi = [0] ## One or more background ROIs. |
179 |
self.calibM = 1 ## Slope used for calibration |
180 |
self.calibB = 0 ## Intercept used for calibration. |
181 |
self.plottable = True ## Whether to plot |
182 |
## Optional injection information. Usually from Clamp or ICM Method files. |
183 |
self.injconc = [0.0] ## One or more analyte concentrations. |
184 |
self.injstart = [0.0] ## One or more injection start times. |
185 |
self.injstop = [0.0] ## One or more injection start times. |
186 |
self.injrind = [0.0] ## One or more expected refractive index jumps. |
187 |
self.flow = 0 ## Flowrate of injection. |
188 |
## Curve fitting information. Model and model parameters. |
189 |
## Example model is reference to function like data=simple1to1(times,params) |
190 |
## Example params = {'Rmax': {'value':1, 'min':0, 'max':10, 'fixed':'float'} } |
191 |
self.model = None ## Model describing this roi. Reference to a function. |
192 |
self.params = [] ## Parameters for this model. A dictionary of dictionaries. |
193 |
return |
194 |
|
195 |
|
196 |
def time2dp(self, t): |
197 |
""" |
198 |
Find datapoint closest to given time. For example:: |
199 |
|
200 |
>>> ba1 = BiosensorArray(1, 300) |
201 |
>>> ba1.roi[0].simulate() |
202 |
>>> print ba1.roi[0].time2dp(40.333) |
203 |
40 |
204 |
""" |
205 |
pos2 = self.time.searchsorted(t) ## Time point just after t. |
206 |
pos2 = min(pos2, len(self.time)-1) ## Avoid indexing error. |
207 |
pos1 = max(pos2-1,0) ## Time point just before t. |
208 |
t2 = abs(self.time[pos2] - t) |
209 |
t1 = abs(self.time[pos1] - t) |
210 |
## Decide if time point just before or just after is closer. |
211 |
if (t2<t1): |
212 |
return pos2 |
213 |
else: |
214 |
return pos1 |
215 |
|
216 |
|
217 |
def time2val(self, t1, t2): |
218 |
""" |
219 |
Return list of SPR values between t1 and t2. For example:: |
220 |
|
221 |
>>> ba1 = BiosensorArray(1, 300) |
222 |
>>> ba1.roi[0].simulate() |
223 |
>>> print min(ba1.roi[0].time2val(35, 65)) |
224 |
1000.0 |
225 |
""" |
226 |
dp1, dp2 = self.time2dp(t1), self.time2dp(t2) |
227 |
return [self.value[i] for i in range(dp1, dp2)] |
228 |
|
229 |
|
230 |
def simulate(self, bulkRI=1e3, kon=1e4, koff=1e-2, conc=1e-6): |
231 |
""" |
232 |
Fill ROI with 300 s of simulated data from simple 1:1 binding. |
233 |
First, there is a brief calibration injection from 35 to 65 s |
234 |
that appears as a square wave of height ``bulkRI``. |
235 |
At time 100 s, there is a jump of 10% ``bulkRI`` followed by |
236 |
binding at a rate dependent on ``kon`` and ``conc``. |
237 |
At time 200 s, dissociation occurs at rate ``koff``. |
238 |
The maximum response (Rmax) is fixed at 1000 units. |
239 |
The model is based on these equations:: |
240 |
|
241 |
alpha = c*kon+koff |
242 |
Req = Rmax*c*kon / alpha |
243 |
R = Req * [1 - exp(-t*alpha)] |
244 |
|
245 |
""" |
246 |
## Error checking. |
247 |
if (koff <= 0): raise Exception, "Simulation koff must be >0" |
248 |
if (kon < 0): raise Exception, "Simulation kon must be >=0" |
249 |
if (conc < 0): raise Exception, "Simulation conc must be >=0" |
250 |
## Start with flat line for 300 sec. |
251 |
self.time = np.arange(300) |
252 |
self.value = np.zeros(300) |
253 |
## Add bulk refractive index jumps. |
254 |
self.value[35:65] += bulkRI |
255 |
self.value[100:199] += (0.1 * bulkRI) |
256 |
## Add binding. |
257 |
rmax = 1000 |
258 |
alpha = conc * kon + koff |
259 |
req = rmax * conc * kon / alpha |
260 |
for t in range(100): |
261 |
self.value[t+100] += req * (1 - np.exp(-t * alpha)) |
262 |
## Add dissociation |
263 |
response = req * (1 - np.exp(-100 * alpha)) |
264 |
for t in range(100): |
265 |
self.value[t+200] += response * np.exp(-t * koff) |
266 |
return |
267 |
|
268 |
"""End of RegOfInterest definition.""" |
269 |
|
270 |
|
271 |
def _test(): |
272 |
""" |
273 |
Automatic Code testing with doctest. |
274 |
Doctest runs the example code in the docstrings. |
275 |
Enable use of ellipses as wildcard for returned text. |
276 |
""" |
277 |
import doctest |
278 |
doctest.testmod(optionflags=doctest.ELLIPSIS) |
279 |
return |
280 |
|
281 |
|
282 |
if __name__ == "__main__": |
283 |
"""Use 'python ba_class.py' for code testing.""" |
284 |
_test() |
285 |
|
286 |
|
287 |
########################### End of module ############################ |