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 < Yuhang Wan, OSPRAI developers
5 < Last modified on 100409 (yymmdd)
6 <
7 < Examples:
8 < #ba1 = BiosensorArray(800,1500)     ## Allocate object with 800 spots of 1500 time points.
9 < #ba1.roi[799].time[1499] = 7502.1   ## Set time in seconds.
10 < #ba1.roi[799].value[1499] = 0.50    ## Set SPR reflectance signal.
11 < #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__ = "100409"
28 > __version__ = "110518"
29 >
30  
31   ## Import libraries.
32   import numpy as np
17 import matplotlib.pyplot as plt
33   from datetime import date
34  
35  
36   class BiosensorArray(object):
37 <    """The Biosensor Array class holds all ROIs/spots from an SPRI microarray run."""
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.
# Line 30 | Line 59
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 ob in self.roi:
68 <            coor2 = (ob.gridx, ob.gridy, ob.spotx, ob.spoty)
69 <            if (coor1==coor2): return ob.uid  ## Success.
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 ob in self.roi: ob.plottable = True
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 ob in self.roi: ob.plottable = False
79 >        for roi in self.roi: roi.plottable = False
80          for i in ilist: self.roi[i].plottable = True
81 <    
82 <    ## A pyplot feature for testing purposes.
83 <    def plot(self):
84 <        """Show plot of the selected sensorgrams."""
85 <        plt.clf()
55 <        plt.title('SPR Data plotted %s' % date.today().isoformat())
56 <        plt.xlabel('Time (s)')
57 <        plt.ylabel('SPR Response')
58 <        plt.grid(True)
59 <        ## Plot traces.
60 <        for ob in self.roi:
61 <            if (ob.plottable==True):
62 <                mylabel = "%i:%s" % (ob.index,ob.name)
63 <                plt.plot(ob.time, ob.value, label=mylabel)
64 <        plt.legend(loc='best')
65 <        plt.show()
66 <    ### End of BiosensorArray class definition.
81 >        return
82 >        
83 >    """End of BiosensorArray class definition."""
84 >
85 >
86  
87  
88   class RegOfInterest(object):
89 <    """This Region Of Interest class can hold one SPR sensorgram."""
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.
# Line 95 | Line 147
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 <        """Find datapoint closest to given time."""
160 <        pos2 = self.time.searchsorted(t)  ## Time point just after t.
161 <        pos1 = max(pos2-1,0)
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 <    ## End of RegOfInterest definition.
184 <
185 <
186 < ## Here are a few lines to test this class.
187 < if __name__ == '__main__':
188 <    print "Starting demo..."
189 <    print "Test plotting..."
117 <    x = BiosensorArray(10,10)
118 <    x.roi[0].time = np.arange(10)
119 <    x.roi[0].value = np.arange(10)**2
120 <    x.roi[1].time = np.arange(10)
121 <    x.roi[1].value = np.arange(10)**1.5
122 <    x.roi[9].time = np.arange(10)
123 <    x.roi[9].value = np.arange(10)**1
124 <    #x.set_plot_all()
125 <    x.set_plot_list([0,1,9])
126 <    x.plot()
127 <    print "Test xy2uid... 9 =",
128 <    x.roi[9].gridx = 2
129 <    print x.xy2uid(2,0,0,0)
130 <    print "Test time2dp... 5.6 =>",
131 <    print x.roi[9].time2dp(5.6)
132 <    print "Demo finished."
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      
135 ################################# End of module #################################
192 +    def simulate(self, bulkRI=1000., kon=10000., koff=.01, conc=.000001):
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 +    """
246 +    Code testing.  
247 +    Simply execute 'python ba_class.py'.
248 +    """
249 +    _test()
250 +
251 +
252 + ########################### End of module ############################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines