ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/ba_class.py
Revision: 65
Committed: Fri May 20 16:33:02 2011 UTC (7 years, 6 months ago) by clausted
File size: 10126 byte(s)
Log Message:
Added much testing to all modules using docstrings and doctest.  To help this, I also added two optional parameters to mclma(), verbose and maxiter.  For example, use "mclma(rois, verbose=False, maxiter=99)" to run the regression up to 99 iterations without echoing progress to the standard output.

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