RacemicMixtureBindingModel and Data

.gitignore

@@ -8,6 +8,7 @@
 *.egg-info
 dist
 build
+workdir
 eggs
 parts
 bin

README.md

@@ -61,13 +61,19 @@ Sample files have been included to test if the library is functional. You can fi
 
 An example for a two-component binding model 
 ```
-bitc_util.py mcmc bitc/testdata/sample.itc -w workdir -v -m TwoComponent --niters=20000 --nburn=1000 --nthin=10 --nfit=100
+python scripts/bitc_util.py mcmc bitc/testdata/sample.itc -w workdir -v -m TwoComponent --niters=20000 --nburn=1000 --nthin=10 --nfit=100
+```
+
+UNDER CONSTRUCTION: Calculate the mean and standard deviation of the log likelihood of your data given previously saved samples from a posterior.
+
+```
+python scripts/bitc_util.py bitc/testdata/sample.itc -w workdir -v -m TwoComponent --calc_logp workdir/sample.h5
 ```
 
 Analyse `.itc` files without mcmc
 
 ```
-bitc_util.py bitc/testdata/sample.itc bitc/testdata/sample2.itc -w workdir
+python scripts/bitc_util.py bitc/testdata/sample.itc bitc/testdata/sample2.itc -w workdir
 ```
 
 An example for a competitive binding model.

bitc/models.py

@@ -3,12 +3,14 @@
 PyMC models to describe ITC binding experiments
 """
 
+import numpy
+import pymc
+
 import copy
 import logging
 from math import exp, log
 
-import numpy
-import pymc
+from bitc.units import ureg, Quantity
 import scipy.integrate
 
 from bitc.units import ureg, Quantity
@@ -23,7 +25,6 @@ class RescalingStep(pymc.StepMethod):
 
     """
     Rescaling StepMethod for sampling correlated changes in ligand and receptor concentration
-
     """
 
     def __init__(self, dictionary, beta, max_scale=1.03, interval=100, verbose=0):
@@ -162,7 +163,6 @@ def __init__(self):
     def _add_unit_to_guesses(value, maximum, minimum, unit):
         """
         Add units to inital guesses for priors
-
         :param value: mean value of the guess
         :type value: float
         :param maximum: maximum for the guess
@@ -230,13 +230,16 @@ def _normal_observation_with_units(name, q_n_model, q_ns, tau, unit):
         """Define a set of normally distributed observations, while stripping units
         :rtype : pymc.Normal
         """
+        # TODO: Ingore N_ignore first injection heats in the likelihood calculation,
+        #       where N_ignore is a variable passed here somehow.
+        # TODO: Make sure that q_n_model and tau are arrays?
+        # TODO: Edit this to ignore N first injection heats in q_n_model and tau?
         return pymc.Normal(name, mu=q_n_model, tau=tau, observed=True, value=q_ns / unit)
 
     @staticmethod
     def _uniform_prior(name, value, maximum, minimum):
         """Define a uniform prior without units
            Added for consistency with other Bindingmodel
-
         :rtype : pymc.Uniform
         """
         return pymc.Uniform(name,
@@ -248,7 +251,6 @@ def _uniform_prior(name, value, maximum, minimum):
     @staticmethod
     def _uniform_prior_with_units(name, value, maximum, minimum, unit):
         """Define a uniform prior, while stripping units
-
         :rtype : pymc.Uniform
         """
         return pymc.Uniform(name,
@@ -262,7 +264,6 @@ def _uniform_prior_with_guesses_and_units(name, value, maximum, minimum, prior_u
         """
         Take initial values, add units or convert units to the right type,
         returns a pymc uniform prior
-
         :rtype : pymc.Uniform
         """
         # Guess has units
@@ -363,7 +364,6 @@ def __init__(self, experiment):
     def expected_injection_heats(V0, DeltaVn, P0, Ls, DeltaG, DeltaH, DeltaH_0, beta, N):
         """
         Expected heats of injection for two-component binding model.
-
         ARGUMENTS
         V0 - cell volume (liter)
         DeltaVn - injection volumes (liter)
@@ -374,12 +374,9 @@ def expected_injection_heats(V0, DeltaVn, P0, Ls, DeltaG, DeltaH, DeltaH_0, beta
         DeltaH_0 - heat of injection (cal)
         beta - inverse temperature * gas constant (mole / kcal)
         N - number of injections
-
         Returns
         -------
         expected injection heats -
-
-
         """
         # TODO Units that go into this need to be verified
         # TODO update docstring with new input
@@ -455,7 +452,12 @@ def _create_rescaling_sampler(self, Ls_stated, P0_stated, experiment):
     def _create_metropolis_sampler(self, Ls_stated, P0_stated, experiment):
         """Create an MCMC sampler for the two component model.
         """
-        mcmc = pymc.MCMC(self, db='ram')
+        dbname = experiment.name+'.h5'
+        import os
+        if os.path.isfile(dbname):
+          mcmc = pymc.MCMC(self, db='ram')
+        else:
+          mcmc = pymc.MCMC(self, db='hdf5', dbname=dbname, dbmode='a')
         mcmc.use_step_method(pymc.Metropolis, self.DeltaG)
         mcmc.use_step_method(pymc.Metropolis, self.DeltaH)
         mcmc.use_step_method(pymc.Metropolis, self.DeltaH_0)
@@ -517,7 +519,6 @@ class CompetitiveBindingModel(BindingModel):
     def __init__(self, experiments, receptor, concentration_uncertainty=0.10):
         """
         ARGUMENTS
-
         experiments (list of Experiment) -
         instrument Instrument that experiment was carried out in (has to be one)
         receptor (string) - name of receptor species
@@ -634,54 +635,35 @@ def __init__(self, experiments, receptor, concentration_uncertainty=0.10):
     def equilibrium_concentrations(Ka_n, C0_R, C0_Ln, V, c0=None):
         """
         Compute the equilibrium concentrations of each complex species for N ligands competitively binding to a receptor.
-
         ARGUMENTS
-
         Ka_n (numpy N-array of float) - Ka_n[n] is the association constant for receptor and ligand species n (1/M)
         x_R (float) - the total number of moles of receptor in the sample volume
         x_n (numpy N-array of float) - x_n[n] is the total number of moles of ligand species n in the sample volume
         V (float) - the total sample volume (L)
-
         RETURNS
-
         C_n (numpy N-array of float) - C_n[n] is the concentration of complex of receptor with ligand species n
-
         EXAMPLES
-
         >>> V = 1.4303e-3 # volume (L)
         >>> x_R = V * 510.e-3 # receptor
         >>> x_Ln = numpy.array([V * 8.6e-6, 200.e-6 * 55.e-6]) # ligands
         >>> Ka_n = numpy.array([1./(400.e-9), 1./(2.e-11)]) # association constants
         >>> from bitc.models import CompetitiveBindingModel
         >>> C_PLn = CompetitiveBindingModel.equilibrium_concentrations(Ka_n, x_R, x_Ln, V)
-
+        
         NOTES
-
         Each complex concentration C_n must obey the relation
-
         Ka_n[n] = C_RLn[n] / (C_R * C_Ln[n])           for n = 1..N
-
         with conservation of mass constraints
-
         V * (C_Ln[n] + C_RLn[n]) = x_Ln[n]             for n = 1..N
-
         and
-
         V * (C_R + C_RLn[:].sum()) = x_R
-
         along with the constraints
-
         0 <= V * C_RLn[n] <= min(x_Ln[n], x_R)         for n = 1..N
         V * C_RLn[:].sum() <= x_R
-
         We can rearrange these expressions to give
-
         V * C_R * C_Ln[n] * Ka_n[n] - V * C_RLn[n] = 0
-
         and eliminate C_Ln[n] and C_R to give
-
         V * (x_R/V - C_RLn[:].sum()) * (x_Ln[n]/V - C_RLn[n]) * Ka_n[n] - V * C_RLn[n] = 0    for n = 1..N
-
         """
 
         x_R = C0_R * V
@@ -921,8 +903,27 @@ def _zero_for_missing__concentrations(self, experiment):
             if species not in experiment.true_syringe_concentration:
                 experiment.true_syringe_concentration[species] = 0.0
 
+class RacemicMixtureBindingModel(CompetitiveBindingModel):
+    """
+    Racemic Mixture Binding Model
+
+    """
+    def equilibrium_concentrations(self, Ka_n, C0_R, C0_Ln, V, c0=None):
+      a = 1./Ka_n[0] + 1./Ka_n[1] + C0_Ln[0] + C0_Ln[1] - C0_R
+      b = 1./Ka_n[1]*(C0_Ln[0]-C0_R) + 1./Ka_n[0]*(C0_Ln[1]-C0_R) + 1./Ka_n[0]*1./Ka_n[1]
+      c = (-1./Ka_n[0])*1./Ka_n[1]*C0_R
+      d = numpy.sqrt(a*a-3*b)
+      theta = numpy.arccos(((-2.*a**3)+9.*a*b-27.*c)/(2.*(d**3)))
+      PA = C0_Ln[0]*(2.*d*numpy.cos(theta/3.)-a)/(3.*1./Ka_n[0]+(2.*d*numpy.cos(theta/3.)-a))
+      PB = C0_Ln[1]*(2.*d*numpy.cos(theta/3.)-a)/(3.*1./Ka_n[1]+(2.*d*numpy.cos(theta/3.)-a))
+      return numpy.array([PA, PB])
+
+#      P = (-a/3)+2./3.*(a**2.-3.*b)**.5*numpy.cos(theta/3.)
+#
+#      A = C0_Ln[0] - PA
+#      B = C0_Ln[1] - PB
 
 # Container of all models that this module provides for use
 known_models = {'TwoComponent': TwoComponentBindingModel,
                 'Competitive': CompetitiveBindingModel,
-                }
+                'RacemicMixture': RacemicMixtureBindingModel}

bitc/parser.py

@@ -19,6 +19,7 @@ def optparser(argv=sys.argv[1:]):
 Usage:
   ITC.py <datafiles>... [-w <workdir> | --workdir=<workdir>] [-n <name> | --name=<name>] [-q <file> | --heats=<file>] [-i <ins> | --instrument=<ins> ] [-v | -vv | -vvv] [-r <file> | --report=<file>] [ -l <logfile> | --log=<logfile>]
   ITC.py mcmc <datafiles>...  (-m <model> | --model=<model>) [-w <workdir> | --workdir=<workdir>] [ -r <receptor> | --receptor=<receptor>] [-n <name> | --name=<name>] [-q <file> | --heats=<file>] [-i <ins> | --instrument=<ins> ] [ -l <logfile> | --log=<logfile>] [-v | -vv | -vvv] [--report=<file>] [options]
+  ITC.py <datafiles>... --calc_logp=<file> (-m <model> | --model=<model>) [-w <workdir> | --workdir=<workdir>] [ -r <receptor> | --receptor=<receptor>] [-n <name> | --name=<name>] [-q <file> | --heats=<file>] [-i <ins> | --instrument=<ins> ] [ -l <logfile> | --log=<logfile>] [-v | -vv | -vvv] [--report=<file>] [options]
   ITC.py (-h | --help)
   ITC.py --license
   ITC.py --version
@@ -41,6 +42,7 @@ def optparser(argv=sys.argv[1:]):
   --nburn=<n>                            No. of Burn-in iterations for mcmc sampling     [default: 500000]
   --nthin=<n>                            Thinning period for mcmc sampling               [default: 250]
   --report=<file>                        Output file with summary in markdown
+  --calc_logp=<file>                     Calculates the log posterior based on samples in the file
 """
     arguments = docopt(__usage__, argv=argv, version='ITC.py, pre-alpha')
     schema = Schema({'--heats': Or(None, And(str, os.path.isfile, Use(os.path.abspath))),  # str, verify that it exists
@@ -67,6 +69,7 @@ def optparser(argv=sys.argv[1:]):
                      # list and ensure it contains existing files
                      '<datafiles>': And(list, lambda inpfiles: [os.path.isfile(inpfile) for inpfile in inpfiles], Use(lambda inpfiles: [os.path.abspath(inpfile) for inpfile in inpfiles])),
                      'mcmc': bool,  # True or False are accepted
+                     '--calc_logp': Or(None, And(str, len)),
                      '--report': Or(None, Use(lambda f: open(f, 'w'))),
                      # Don't use, or open file with writing permissions
                      '--log': Or(None, str),  # Don't use, or str

data/ADP-AdK/ADP25AdK5.DAT

@@ -0,0 +1,32 @@
+DH	INJV	Xt	Mt	XMt	NDH
+-796.59036	8	0	0.212	0.65971	-3982.95179
+-437.85523	8	0.13908	0.21082	1.3231	-2189.27615
+-242.85234	8	0.27738	0.20965	1.99015	-1214.2617
+-196.75705	8	0.41491	0.20848	2.66087	-983.78527
+-183.99409	8	0.55166	0.20732	3.33526	-919.97043
+-183.15421	8	0.68762	0.20617	4.01333	-915.77105
+-185.4336	8	0.82281	0.20502	4.69506	-927.16801
+-189.43672	8	0.95723	0.20388	5.38046	-947.18359
+-192.87694	8	1.09086	0.20274	6.06954	-964.38468
+-197.66906	8	1.22372	0.20162	6.76228	-988.34532
+-200.44818	8	1.3558	0.20049	7.4587	-1002.2409
+-205.05612	8	1.4871	0.19938	8.15878	-1025.28061
+-207.76387	8	1.61762	0.19827	8.86254	-1038.81934
+-212.38287	8	1.74736	0.19716	9.56996	-1061.91436
+-215.14892	8	1.87633	0.19606	10.28106	-1075.74462
+-218.15501	8	2.00452	0.19497	10.99582	-1090.77504
+-221.60966	8	2.13193	0.19389	11.71426	-1108.04828
+-224.89212	8	2.25856	0.1928	12.43636	-1124.46061
+-227.32612	8	2.38441	0.19173	13.16214	-1136.63058
+-230.5682	8	2.50949	0.19066	13.89158	-1152.841
+-232.13219	8	2.63379	0.1896	14.6247	-1160.66094
+-235.35577	8	2.7573	0.18854	15.36148	-1176.77886
+-237.71969	8	2.88005	0.18748	16.10194	-1188.59846
+-240.20598	8	3.00201	0.18644	16.84607	-1201.02991
+-242.2048	8	3.12319	0.1854	17.59386	-1211.02398
+-245.04354	8	3.2436	0.18436	18.34533	-1225.2177
+-246.40309	8	3.36323	0.18333	19.10047	-1232.01545
+-249.28787	8	3.48208	0.1823	19.85927	-1246.43933
+-250.41005	8	3.60015	0.18128	20.62175	-1252.05023
+-252.48003	8	3.71745	0.18027	21.3879	-1262.40017
+	--	3.83397	0.17926	--