ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/ba_class.py
(Generate patch)
# Line 1 | Line 1
1   """
2 < ba: Biosensor Array class for storing SPRI data.
3 < Christopher Lausted, Institute for Systems Biology
4 < Last modified on 100408 (yymmdd)
5 <
6 < Examples:
7 < #ba1 = BiosensorArray(800,1500)     ## Allocate object with 800 spots of 1500 time points.
8 < #ba1.roi[799].time[1499] = 7502.1   ## Set time in seconds.
9 < #ba1.roi[799].value[1499] = 0.50    ## Set SPR reflectance signal.
10 < #ba1.roi[799].name = "anti-IgG"     ## Set name of microarray feature.
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__ = "100408"
28 > __version__ = "111012"
29 >
30  
31   ## Import libraries.
32   import numpy as np
33 < import matplotlib.pyplot as plt
33 > import cPickle as pickle
34   from datetime import date
35  
36  
37   class BiosensorArray(object):
38 <    """The Biosensor Array class holds all ROIs/spots from an SPRI microarray run."""
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.
# Line 29 | Line 60
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 ob in self.roi:
69 <            coor2 = (ob.gridx, ob.gridy, ob.spotx, ob.spoty)
70 <            if (coor1==coor2): return ob.uid  ## Success.
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 ob in self.roi: ob.plottable = True
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 ob in self.roi: ob.plottable = False
80 >        for roi in self.roi: roi.plottable = False
81          for i in ilist: self.roi[i].plottable = True
82 <    
83 <    ## A pyplot feature for testing purposes.
84 <    def plot(self):
85 <        """Show plot of the selected sensorgrams."""
86 <        plt.clf()
87 <        plt.title('SPR Data plotted %s' % date.today().isoformat())
88 <        plt.xlabel('Time (s)')
89 <        plt.ylabel('SPR Response')
90 <        plt.grid(True)
91 <        ## Plot traces.
92 <        for ob in self.roi:
93 <            if (ob.plottable==True):
94 <                mylabel = "%i:%s" % (ob.index,ob.name)
95 <                plt.plot(ob.time, ob.value, label=mylabel)
96 <        plt.legend(loc='best')
97 <        plt.show()
98 <    ### End of BiosensorArray class definition.
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 <    """This Region Of Interest class can hold one SPR sensorgram."""
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.
# Line 94 | Line 185
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 <        """Find datapoint closest to given time."""
198 <        pos2 = self.time.searchsorted(t)  ## Time point just after t.
199 <        pos1 = max(pos2-1,0)
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 <    ## End of RegOfInterest definition.
222 <
223 <
224 < ## Here are a few lines to test this class.
225 < if __name__ == '__main__':
226 <    print "Starting demo..."
227 <    print "Test plotting..."
116 <    x = BiosensorArray(10,10)
117 <    x.roi[0].time = np.arange(10)
118 <    x.roi[0].value = np.arange(10)**2
119 <    x.roi[1].time = np.arange(10)
120 <    x.roi[1].value = np.arange(10)**1.5
121 <    x.roi[9].time = np.arange(10)
122 <    x.roi[9].value = np.arange(10)**1
123 <    #x.set_plot_all()
124 <    x.set_plot_list([0,1,9])
125 <    x.plot()
126 <    print "Test xy2uid... 9 =",
127 <    x.roi[9].gridx = 2
128 <    print x.xy2uid(2,0,0,0)
129 <    print "Test time2dp... 5.6 =>",
130 <    print x.roi[9].time2dp(5.6)
131 <    print "Demo finished."
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 < ################################# End of module #################################
231 < """
232 < class classSprArray:
233 <    time = np.zeros(0, dtype=float)  # 0,5,10...seconds.
234 <    sproi = []  # List of classSprRoi objects.
235 <    _calibrated = False  # Private variable.
236 <    _bgsubtracted = False  # Private variable.
237 <    
238 <    ## Apply calibration factors to all data.
239 <    def calibrate(self):
240 <        self._calibrated = True
241 <        for i in self.sproi:
242 <            i.oldvalue = i.value.copy()
243 <            i.value = i.value * i.calibM + i.calibB
244 <    def calibrated(self):
245 <        return self._calibrated
246 <        
247 <    ## Apply background subtraction
248 <    def bgsubtract(self):
249 <        self._bgsubtracted = True
250 <        for i in self.sproi:
251 <            i.oldvalue = i.value.copy()
252 <        for i in self.sproi:
253 <            if (i.bgroi<=len(self.sproi)):
254 <                i.value = i.value - self.sproi[i.bgroi].oldvalue
255 <            else:
256 <                print "Index out of range in classSprArray"
257 <    def bgsubtracted(self):
258 <        return self._bgsubtracted
259 <        
260 <    def addsproi(self):
261 <        newroi = classSprRoi()
262 <        self.sproi.append(newroi)
263 <
264 <    ## End classSprArray definition.
265 <
266 < Unused examples:
267 <
268 < class Person(object):
269 <    def get_full_name(self):
270 <        return "%s %s" % (self.first_name, self.last_name)
271 <    def set_full_name(self, full_name):
272 <        self.first_name, self.last_name = full_name.split()
273 <    full_name = property(get_full_name, set_full_name)
274 <
275 < class Example(object):
276 <    def __init__(self):
277 <        self.data = []
278 <    
279 <    @property
280 <    def something(self):
281 <        return whatever
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  
187 """
287 + ########################### End of module ############################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines