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 User Rev File contents
1 clausted 7 """
2 clausted 47 ba_class
3     --------
4 clausted 7
5 clausted 47 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 clausted 48
17 clausted 47 >>> 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 clausted 58 >>> 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 clausted 64 >>> print round(ba1.roi[0].value[299], 1)
26     160.6
27 clausted 7 """
28 clausted 66 __version__ = "111012"
29 clausted 7
30 clausted 47
31 clausted 7 ## Import libraries.
32     import numpy as np
33 clausted 66 import cPickle as pickle
34 clausted 7 from datetime import date
35    
36    
37     class BiosensorArray(object):
38 clausted 27 """
39     The Biosensor Array class holds all ROIs (spots)
40     from an SPRI microarray run.
41    
42 clausted 47 ================ =============== =======================================
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 clausted 27 """
52    
53 clausted 7 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 clausted 47 return
64 clausted 7
65     def xy2uid(self, gx, gy, x, y):
66     """Given roi coordinates, find corresponding uid (index)."""
67     coor1 = (gx, gy, x, y)
68 clausted 50 for roi in self.roi:
69     coor2 = (roi.gridx, roi.gridy, roi.spotx, roi.spoty)
70     if (coor1==coor2): return roi.uid ## Success.
71 clausted 7 return -1 ## Failure.
72    
73     def set_plot_all(self):
74     """Choose to plot every sensorgram."""
75 clausted 50 for roi in self.roi: roi.plottable = True
76 clausted 47 return
77 clausted 7
78     def set_plot_list(self, ilist):
79     """Choose a list of sensorgrams to plot"""
80 clausted 50 for roi in self.roi: roi.plottable = False
81 clausted 7 for i in ilist: self.roi[i].plottable = True
82 clausted 47 return
83 clausted 66
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 clausted 47
101 clausted 66 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 clausted 47 """End of BiosensorArray class definition."""
124 clausted 7
125    
126     class RegOfInterest(object):
127 clausted 27 """
128     This Region Of Interest class can hold one SPR sensorgram.
129    
130 clausted 47 =========== =============== =============================================
131 clausted 27 Property Type Description
132 clausted 47 ----------- --------------- ---------------------------------------------
133 clausted 27 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 clausted 47 =========== =============== =============================================
159 clausted 27 """
160    
161 clausted 7 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 clausted 17 ## Curve fitting information. Model and model parameters.
189     ## Example model is reference to function like data=simple1to1(times,params)
190 clausted 28 ## Example params = {'Rmax': {'value':1, 'min':0, 'max':10, 'fixed':'float'} }
191 clausted 17 self.model = None ## Model describing this roi. Reference to a function.
192     self.params = [] ## Parameters for this model. A dictionary of dictionaries.
193 clausted 64 return
194 clausted 7
195 clausted 64
196 clausted 7 def time2dp(self, t):
197 clausted 64 """
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 clausted 13 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 clausted 7 t2 = abs(self.time[pos2] - t)
209     t1 = abs(self.time[pos1] - t)
210 clausted 13 ## Decide if time point just before or just after is closer.
211 clausted 7 if (t2<t1):
212     return pos2
213     else:
214     return pos1
215 clausted 64
216    
217 clausted 12 def time2val(self, t1, t2):
218 clausted 64 """
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 clausted 12 dp1, dp2 = self.time2dp(t1), self.time2dp(t2)
227 clausted 13 return [self.value[i] for i in range(dp1, dp2)]
228 clausted 64
229    
230 clausted 65 def simulate(self, bulkRI=1e3, kon=1e4, koff=1e-2, conc=1e-6):
231 clausted 58 """
232 clausted 64 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 clausted 61
241 clausted 58 alpha = c*kon+koff
242     Req = Rmax*c*kon / alpha
243     R = Req * [1 - exp(-t*alpha)]
244 clausted 61
245 clausted 58 """
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 clausted 12
268 clausted 47 """End of RegOfInterest definition."""
269 clausted 7
270 clausted 66
271 clausted 58 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 clausted 64 return
280 clausted 58
281    
282     if __name__ == "__main__":
283 clausted 65 """Use 'python ba_class.py' for code testing."""
284 clausted 58 _test()
285    
286    
287     ########################### End of module ############################