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 (8 years, 4 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 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 64 __version__ = "110518"
29 clausted 7
30 clausted 47
31 clausted 7 ## Import libraries.
32     import numpy as np
33     from datetime import date
34    
35    
36     class BiosensorArray(object):
37 clausted 27 """
38     The Biosensor Array class holds all ROIs (spots)
39     from an SPRI microarray run.
40    
41 clausted 47 ================ =============== =======================================
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 clausted 27 """
51    
52 clausted 7 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 clausted 47 return
63 clausted 7
64     def xy2uid(self, gx, gy, x, y):
65     """Given roi coordinates, find corresponding uid (index)."""
66     coor1 = (gx, gy, x, y)
67 clausted 50 for roi in self.roi:
68     coor2 = (roi.gridx, roi.gridy, roi.spotx, roi.spoty)
69     if (coor1==coor2): return roi.uid ## Success.
70 clausted 7 return -1 ## Failure.
71    
72     def set_plot_all(self):
73     """Choose to plot every sensorgram."""
74 clausted 50 for roi in self.roi: roi.plottable = True
75 clausted 47 return
76 clausted 7
77     def set_plot_list(self, ilist):
78     """Choose a list of sensorgrams to plot"""
79 clausted 50 for roi in self.roi: roi.plottable = False
80 clausted 7 for i in ilist: self.roi[i].plottable = True
81 clausted 47 return
82    
83     """End of BiosensorArray class definition."""
84 clausted 7
85    
86 clausted 64
87    
88 clausted 7 class RegOfInterest(object):
89 clausted 27 """
90     This Region Of Interest class can hold one SPR sensorgram.
91    
92 clausted 47 =========== =============== =============================================
93 clausted 27 Property Type Description
94 clausted 47 ----------- --------------- ---------------------------------------------
95 clausted 27 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 clausted 47 =========== =============== =============================================
121 clausted 27 """
122    
123 clausted 7 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 clausted 17 ## Curve fitting information. Model and model parameters.
151     ## Example model is reference to function like data=simple1to1(times,params)
152 clausted 28 ## Example params = {'Rmax': {'value':1, 'min':0, 'max':10, 'fixed':'float'} }
153 clausted 17 self.model = None ## Model describing this roi. Reference to a function.
154     self.params = [] ## Parameters for this model. A dictionary of dictionaries.
155 clausted 64 return
156 clausted 7
157 clausted 64
158 clausted 7 def time2dp(self, t):
159 clausted 64 """
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 clausted 13 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 clausted 7 t2 = abs(self.time[pos2] - t)
171     t1 = abs(self.time[pos1] - t)
172 clausted 13 ## Decide if time point just before or just after is closer.
173 clausted 7 if (t2<t1):
174     return pos2
175     else:
176     return pos1
177 clausted 64
178    
179 clausted 12 def time2val(self, t1, t2):
180 clausted 64 """
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 clausted 12 dp1, dp2 = self.time2dp(t1), self.time2dp(t2)
189 clausted 13 return [self.value[i] for i in range(dp1, dp2)]
190 clausted 64
191    
192 clausted 65 def simulate(self, bulkRI=1e3, kon=1e4, koff=1e-2, conc=1e-6):
193 clausted 58 """
194 clausted 64 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 clausted 61
203 clausted 58 alpha = c*kon+koff
204     Req = Rmax*c*kon / alpha
205     R = Req * [1 - exp(-t*alpha)]
206 clausted 61
207 clausted 58 """
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 clausted 12
230 clausted 47 """End of RegOfInterest definition."""
231 clausted 7
232 clausted 50
233 clausted 58 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 clausted 64 return
242 clausted 58
243    
244     if __name__ == "__main__":
245 clausted 65 """Use 'python ba_class.py' for code testing."""
246 clausted 58 _test()
247    
248    
249     ########################### End of module ############################