ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/ba_class.py
Revision: 66
Committed: Thu Oct 13 19:49:53 2011 UTC (7 years, 11 months ago) by clausted
File size: 11559 byte(s)
Log Message:
Add methods BiosensorArray.save(filename) and .load(filename) which save serialized data files with the sba extension.
Line File contents
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 ############################