osprai/trunk/ba_class.py

"""
ba_class
--------

Biosensor Array class for storing SPRI data.
This contains an array of SPR sensorgrams with their
associated models and parameters.
This can also contain the microarray geometery.

.. moduleauthor:: Christopher Lausted, 
                  Institute for Systems Biology, 
                  Yuhang Wan,
                  OSPRAI developers.   
                  
Examples::

  >>> ba1 = BiosensorArray(800,1500)     ## Allocate object with 800 spots of 1500 time points.
  >>> ba1.roi[799].time[1499] = 7502.1   ## Set time in seconds.
  >>> ba1.roi[799].value[1499] = 0.50    ## Set SPR reflectance signal. 
  >>> ba1.roi[799].name = "anti-IgG"     ## Set name of microarray feature.
  >>> print ba1.roi[799].name
  anti-IgG
  >>> ba1 = BiosensorArray(1, 300)
  >>> ba1.roi[0].simulate(bulkRI=1000, kon=1e4, koff=1e-2, conc=1e-6)
  >>> print round(ba1.roi[0].value[299], 1)
  160.6
"""
__version__ = "111012"


## Import libraries.
import numpy as np
import cPickle as pickle
from datetime import date


class BiosensorArray(object):
    """
    The Biosensor Array class holds all ROIs (spots)
    from an SPRI microarray run.
    
    ================ =============== =======================================
    Property         Type            Description
    ---------------- --------------- ---------------------------------------
    roi              list of object  List of RegOfInterest objects.
    rois             integer         Size of roi (Number of ROIs).
    dpoints          integer         Number of datapoints in each ROI.
    primarydatafiles list of string  List of file names.
    comments         list of string  List of comments regarding experiments.
    ================ =============== =======================================
    """
    
    def __init__(self, roi_size, datapoint_size):
        """This object is dimensioned upon initialization.""" ## Might improve speed.
        ## Use a list comprehension to dimension list of classes.
        self.roi = [RegOfInterest(i,datapoint_size) for i in range(roi_size)] 
        ## Remember the dimensions for quick bounds-checking later.
        self.rois = roi_size
        self.dpoints = datapoint_size
        ## Use these to keep track of laboratory notes.
        self.primarydatafiles = [""]  ## Like "spritdata.txt"
        self.comments = [""]          ## Like "Antibody data from 1 Jan 2010"
        return
        
    def xy2uid(self, gx, gy, x, y):
        """Given roi coordinates, find corresponding uid (index)."""
        coor1 = (gx, gy, x, y)
        for roi in self.roi:
            coor2 = (roi.gridx, roi.gridy, roi.spotx, roi.spoty)
            if (coor1==coor2): return roi.uid  ## Success.
        return -1  ## Failure.
        
    def set_plot_all(self):
        """Choose to plot every sensorgram."""
        for roi in self.roi: roi.plottable = True
        return
        
    def set_plot_list(self, ilist):
        """Choose a list of sensorgrams to plot"""
        for roi in self.roi: roi.plottable = False
        for i in ilist: self.roi[i].plottable = True
        return
       
    def save(self, filename="default.sba"):
        """
        Save (serialize) this ba object as a file.
        Use the extension '.sba' and designate it file format 1.
        Use cPickle protocol 2 (binary).
        """
        if (filename[-4:] != '.sba'): filename = filename + ".sba"
        fp = open(filename, 'wb')
        protocol = 2  ## Use cPickle protocol 2 (binary).
        formatver = 1  ## My format version 1.
        pickle.dump(formatver, fp, protocol)
        pickle.dump(self.comments, fp, protocol)
        pickle.dump(self.primarydatafiles, fp, protocol)
        pickle.dump(self.dpoints, fp, protocol)
        pickle.dump(self.rois, fp, protocol)
        pickle.dump(self.roi, fp, protocol)
        
        fp.close()
        return
        """ End of save() method"""

    def load(self, filename="default.sba"):
        """
        Load data from sba file into this ba object.
        """
        fp = open(filename, 'rb')
        formatver = pickle.load(fp)
        if (formatver == 1):
            self.comments = pickle.load(fp)
            self.primarydatafiles = pickle.load(fp)
            self.dpoints = pickle.load(fp)
            self.rois = pickle.load(fp)
            self.roi = pickle.load(fp)
        else:
            print "Error: file format version is not 1 but " + formatver
        fp.close()
        return
        """ End of load() method"""
        
    """End of BiosensorArray class definition."""


class RegOfInterest(object):
    """
    This Region Of Interest class can hold one SPR sensorgram.
    
    =========== =============== =============================================
    Property    Type            Description
    ----------- --------------- ---------------------------------------------
    uid         integer         Unique ID integer. Never changes.
    index       integer         Changes if ba size changes (add/remove rois).
    time        nparray         Usually seconds.
    value       nparray         Arbitrary units.
    dpoints     integer         Can use for bounds checking.
    name        string          Name of ROI.
    desc        string          Description of ROI
    
    gridx       integer         Can be Block number from GAL file.  
    gridy       integer         Unused?
    spotx       integer         ROI column number.
    spoty       integer         ROI row number.
    bgroi       list of integer One or more background ROIs.
    calibM      float           Slope used for calibration
    calibB      float           Intercept used for calibration.
    plottable   boolean         Whether to plot
    
    injconc     list of float   One or more analyte concentrations.
    injstart    list of float   One or more injection start times.
    injstop     list of float   One or more injection start times.
    injrind     list of float   One or more expected refractive index jumps.
    
    flow        float           Flowrate of injection.
    model       ref to funct    Model describing this ROI. 
    params      dict of dict    Parameters for this model. 
    =========== =============== =============================================
    """
    
    def __init__(self, i, datapoint_size):
        """
        Each ROI in a list needs a unique ID number which will be i+1.
        The time and value arrays must be allocated using datapoint_size.
        """
        self.uid = i+1            ## Unique ID integer. Never changes.
        self.index = i            ## Changes if ba size changes (add/remove rois).
        self.time = np.zeros(datapoint_size, dtype=float)  ## Usually seconds.
        self.value = np.zeros(datapoint_size, dtype=float) ## Arbitrary units.
        self.dpoints = datapoint_size  ## Can use for bounds checking.
        ## Useful optional information. Usually from GAL or Key files.
        self.name = "Spot%i" % (i+1)
        self.desc = "Description"
        self.gridx = 0            ## Can be Block number from GAL file.  
        self.gridy = 0            ## Unused?
        self.spotx = 0            ## ROI column number.
        self.spoty = 0            ## ROI row number.
        self.bgroi = [0]          ## One or more background ROIs.
        self.calibM = 1           ## Slope used for calibration
        self.calibB = 0           ## Intercept used for calibration.
        self.plottable = True     ## Whether to plot
        ## Optional injection information. Usually from Clamp or ICM Method files.
        self.injconc = [0.0]     ## One or more analyte concentrations.
        self.injstart = [0.0]    ## One or more injection start times.
        self.injstop = [0.0]     ## One or more injection start times.
        self.injrind = [0.0]     ## One or more expected refractive index jumps.
        self.flow = 0            ## Flowrate of injection.
        ## Curve fitting information.  Model and model parameters.
        ## Example model is reference to function like data=simple1to1(times,params)
        ## Example params = {'Rmax': {'value':1, 'min':0, 'max':10, 'fixed':'float'} } 
        self.model = None    ## Model describing this roi. Reference to a function.
        self.params = []     ## Parameters for this model. A dictionary of dictionaries.
        return
        
      
    def time2dp(self, t):
        """
        Find datapoint closest to given time. For example::
        
          >>> ba1 = BiosensorArray(1, 300)
          >>> ba1.roi[0].simulate()
          >>> print ba1.roi[0].time2dp(40.333)
          40
        """
        pos2 = self.time.searchsorted(t)    ## Time point just after t.
        pos2 = min(pos2, len(self.time)-1)  ## Avoid indexing error.
        pos1 = max(pos2-1,0)                ## Time point just before t.
        t2 = abs(self.time[pos2] - t)
        t1 = abs(self.time[pos1] - t)
        ## Decide if time point just before or just after is closer.
        if (t2<t1): 
            return pos2
        else: 
            return pos1
    
    
    def time2val(self, t1, t2):
        """
        Return list of SPR values between t1 and t2. For example::
        
          >>> ba1 = BiosensorArray(1, 300)
          >>> ba1.roi[0].simulate()
          >>> print min(ba1.roi[0].time2val(35, 65))
          1000.0
        """
        dp1, dp2 = self.time2dp(t1), self.time2dp(t2)
        return [self.value[i] for i in range(dp1, dp2)]
    
    
    def simulate(self, bulkRI=1e3, kon=1e4, koff=1e-2, conc=1e-6):
        """
        Fill ROI with 300 s of simulated data from simple 1:1 binding.
        First, there is a brief calibration injection from 35 to 65 s
        that appears as a square wave of height ``bulkRI``.
        At time 100 s, there is a jump of 10% ``bulkRI`` followed by
        binding at a rate dependent on ``kon`` and ``conc``.
        At time 200 s, dissociation occurs at rate ``koff``.
        The maximum response (Rmax) is fixed at 1000 units.
        The model is based on these equations::
        
          alpha = c*kon+koff
          Req = Rmax*c*kon / alpha
          R = Req * [1 - exp(-t*alpha)]
          
        """
        ## Error checking.
        if (koff <= 0): raise Exception, "Simulation koff must be >0"
        if (kon < 0): raise Exception, "Simulation kon must be >=0"
        if (conc < 0): raise Exception, "Simulation conc must be >=0"
        ## Start with flat line for 300 sec.
        self.time = np.arange(300)
        self.value = np.zeros(300)
        ## Add bulk refractive index jumps.
        self.value[35:65] += bulkRI
        self.value[100:199] += (0.1 * bulkRI)
        ## Add binding.
        rmax = 1000
        alpha = conc * kon + koff
        req = rmax * conc * kon / alpha
        for t in range(100):
             self.value[t+100] += req * (1 - np.exp(-t * alpha))
        ## Add dissociation
        response = req * (1 - np.exp(-100 * alpha))
        for t in range(100):
            self.value[t+200] += response * np.exp(-t * koff)
        return
    
    """End of RegOfInterest definition."""

        
def _test():
    """
    Automatic Code testing with doctest.
    Doctest runs the example code in the docstrings.
    Enable use of ellipses as wildcard for returned text.
    """
    import doctest
    doctest.testmod(optionflags=doctest.ELLIPSIS)
    return


if __name__ == "__main__":
    """Use 'python ba_class.py' for code testing."""
    _test()


########################### End of module ############################
Revision:	66
Committed:	Thu Oct 13 19:49:53 2011 UTC (8 years, 1 month 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 = Rmaxckon / 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 ############################