Skip to content

Commit

Permalink
docs: small fixes from test procedure
Browse files Browse the repository at this point in the history
  • Loading branch information
JoepVanlier committed Jan 8, 2025
1 parent 6d06752 commit 78a58f4
Show file tree
Hide file tree
Showing 17 changed files with 66 additions and 60 deletions.
35 changes: 14 additions & 21 deletions docs/examples/bead_coupling/coupling.rst
Original file line number Diff line number Diff line change
Expand Up @@ -293,17 +293,10 @@ Next, we define a function that calibrates all axes with both passive and active

`calculate_calibrations` now returns the distances between the beads and all the calibration factors obtained with passive and active calibration.

>>> calculate_calibrations(lk.File(calibration_data[0][0]), bead_diameter=2.1, temperature=25)
{'dx': 22.107661709308474,
'dy': -0.07111526715106109,
'passive': {'1x': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d024aaa0>,
'2x': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d06df130>,
'1y': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d0297010>,
'2y': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d04992d0>},
'active': {'1x': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d056fe80>,
'2x': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d053b370>,
'1y': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d04992a0>,
'2y': <lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults at 0x2d049cfd0>}}
>>> single_distance_result = calculate_calibrations(lk.File(calibration_data[0][0]), bead_diameter=2.1, temperature=25)
... single_distance_result["passive"]["1x"]

.. image:: single_result.png

Let's calculate calibration factors for all the data in this dataset.
Note that this cell may take a while to execute as it is performing `240` calibrations::
Expand Down Expand Up @@ -347,7 +340,7 @@ Now that we have those results, let's define some functions to conveniently extr
We can now show the effect of coupling on active calibration in practice using the analyzed data::

parameters = {
"Rd": "Displacement sensitivity [$\mu$m/V]",
"Rd": r"Displacement sensitivity [$\mu$m/V]",
"Rf": "Force sensitivity [pN/V]",
"kappa": "Stiffness [pN/nm]",
}
Expand All @@ -358,7 +351,7 @@ We can now show the effect of coupling on active calibration in practice using t
dx, dy = extract_distances(experiment[0])
plt.plot(dx, extract_parameter(experiment[0], "active", "2x", param), ".", label="active")
plt.plot(dx, extract_parameter(experiment[0], "passive", "2x", param), "x", label="passive")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel(param_description)

plt.tight_layout()
Expand Down Expand Up @@ -405,7 +398,7 @@ To show how well this model fits the data, we can plot it alongside the ratio of
plt.plot(dx, ac / pc, f"C{trap}.")
plt.plot(dx_c, 1 / c, "k--")
plt.axvline(bead_diameter, color="lightgray", linestyle="--")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel("$R_{d, ac} / R_{d, passive}$")
plt.title(f"Displacement sensitivity ratio {axis} AC/PC")

Expand All @@ -415,7 +408,7 @@ To show how well this model fits the data, we can plot it alongside the ratio of
plt.plot(dx, ac / pc, f"C{trap}.")
plt.plot(dx_c, c, "k--")
plt.axvline(bead_diameter, color="lightgray", linestyle="--")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel("$R_{f, ac} / R_{f, passive}$")
plt.title(f"Force sensitivity ratio {axis} AC/PC")

Expand All @@ -425,8 +418,8 @@ To show how well this model fits the data, we can plot it alongside the ratio of
plt.plot(dx, ac / pc, f"C{trap}.", label=f"{trap}{axis}")
plt.plot(dx_c, c**2, "k--")
plt.axvline(bead_diameter, color="lightgray", linestyle="--")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.ylabel("$\kappa_{ac} / \kappa_{passive}$")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel(r"$\kappa_{ac} / \kappa_{passive}$")
plt.title(f"Stiffness sensitivity ratio {axis} AC/PC")

plt.tight_layout()
Expand Down Expand Up @@ -461,7 +454,7 @@ Some remaining variability is expected as the bead radius (which is subject to v
pc = extract_parameter(exp, "passive", f"{trap}{axis}", "Rd")
plt.plot(dx, ac * c / pc, f"C{trap}.")
plt.axvline(bead_diameter, color="lightgray", linestyle="--")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel("$R_{d, ac} / R_{d, passive}$")
plt.title(f"Displacement sensitivity ratio {axis} AC/PC")

Expand All @@ -470,7 +463,7 @@ Some remaining variability is expected as the bead radius (which is subject to v
pc = extract_parameter(exp, "passive", f"{trap}{axis}", "Rf")
plt.plot(dx, ac / c / pc, f"C{trap}.")
plt.axvline(bead_diameter, color="lightgray", linestyle="--")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel("$R_{f, ac} / R_{f, passive}$")
plt.title(f"Force sensitivity ratio {axis} AC/PC")

Expand All @@ -479,8 +472,8 @@ Some remaining variability is expected as the bead radius (which is subject to v
pc = extract_parameter(exp, "passive", f"{trap}{axis}", "kappa")
plt.plot(dx, ac / c**2 / pc, f"C{trap}.", label=f"{trap}{axis}")
plt.axvline(bead_diameter, color="lightgray", linestyle="--")
plt.xlabel("Bead-Bead Distance [$\mu$m]")
plt.ylabel("$\kappa_{ac} / \kappa_{passive}$")
plt.xlabel(r"Bead-Bead Distance [$\mu$m]")
plt.ylabel(r"$\kappa_{ac} / \kappa_{passive}$")
plt.title(f"Stiffness sensitivity ratio {axis} AC/PC")

plt.tight_layout()
Expand Down
3 changes: 3 additions & 0 deletions docs/examples/bead_coupling/single_result.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
40 changes: 20 additions & 20 deletions docs/examples/binding_lifetime/binding_lifetime.rst
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ Binding lifetime analysis
Determine the binding lifetime from a kymograph
-----------------------------------------------

In this Notebook, we will determine the binding lifetime of a fluorescently labeled protein binding to DNA. The binding lifetime is also referred to as the binding duration
In this Notebook, we will determine the binding lifetime of a fluorescently labeled protein binding to DNA. The binding lifetime is also referred to as the binding duration
and relates to the off rate as :math:`k_{off}` = 1/binding lifetime.

First, we will track the binding events. Then we collect the binding durations and use maximum likelihood fitting to fit an exponential function to the data.
Expand All @@ -33,7 +33,7 @@ Load the first file and plot the kymograph::

file1 = lk.File("test_data/kymo1.h5")
_, kymo1 = file1.kymos.popitem()

plt.figure()
kymo1.plot("g", aspect = 5, adjustment=lk.ColorAdjustment([0], [5]))

Expand All @@ -43,7 +43,7 @@ Load and plot the second kymograph::

file2 = lk.File("test_data/kymo2.h5")
_, kymo2 = file2.kymos.popitem()

plt.figure()
kymo2.plot("g", aspect = 5, adjustment=lk.ColorAdjustment([0], [5]))

Expand Down Expand Up @@ -100,7 +100,7 @@ The (loaded) tracks can be plotted on top of the original kymograph to visualize
.. image:: tracks1.png

Note that two tracks at t=0 were manually removed, because the starting points of these tracks are not visible. This means that we cannot determine the duration of these tracks.
The length of each track corresponds to the duration of a binding event. As can be seen from the above image, there is a large variation in track lengths.
The length of each track corresponds to the duration of a binding event. As can be seen from the above image, there is a large variation in track lengths.
By collecting all these track durations into a 'binding lifetime distribution', we can analyze the binding lifetime in more detail.

Combine tracks
Expand Down Expand Up @@ -139,14 +139,14 @@ Fit a single exponential to the dwell times and plot the result::
The fitted lifetime :math:`\tau = 4` seconds.

The parameter `n_components` indicates the number of exponential time scales in the fit, as further explained below.
The parameters `observed_minimum` and `discrete_model` are further explained in :ref:`dwelltime_analysis`.
The parameters `observed_minimum` and `discrete_model` are further explained in :ref:`dwelltime_analysis`.

Double exponential fit
^^^^^^^^^^^^^^^^^^^^^^

Sometimes, the distribution can best be fit by multiple exponential time scales.
These exponential time scales reveal something about the underlying mechanism of binding.
For example, the protein of interest binds with higher affinity to the target site, while it binds more transiently to off-target sites.
Sometimes, the distribution can best be fit by multiple exponential time scales.
These exponential time scales reveal something about the underlying mechanism of binding.
For example, the protein of interest binds with higher affinity to the target site, while it binds more transiently to off-target sites.
Such behavior has been observed for various proteins such as Cas9 [1]_.

In binding lifetime analysis, it is therefore important to test which number of exponentials optimally fits the data.
Expand All @@ -167,10 +167,10 @@ Fit a double exponential distribution to the binding lifetimes by setting `n_com
.. _double_exponential_fit:
.. image:: double_exponential_fit.png

The component :math:`a_1=0.94` with lifetime :math:`\tau_1 = 3.2` seconds, while component :math:`a_2=0.06` with lifetime :math:`\tau_2 = 24` seconds.
The component :math:`a_1=0.94` with lifetime :math:`\tau_1 = 3` seconds, while component :math:`a_2=0.059` with lifetime :math:`\tau_2 = 18` seconds.

Next we have to select which is the optimal model: 1 or 2 exponential time scales.
There are various methods for model selection. We will discuss 3 of them below.
There are various methods for model selection. We will discuss 3 of them below.

Confidence intervals and model comparison
-----------------------------------------
Expand All @@ -190,7 +190,7 @@ The :ref:`pop_confidence_intervals`, can be used to judge how precisely we can e
.. image:: profile1.png

The parameter to be fitted is given on the x-axis of the plots and the optimal value is where the curve is at its minimum.
The lower the :math:`\chi^2` value at the minimum, the better the fit.
The lower the :math:`\chi^2` value at the minimum, the better the fit.
The point where the profile crosses the dashed horizontal line is an indication for the 95% confidence interval.

The profile likelihood for the single exponent looks parabolic and is almost symmetric, which indicates that the estimate of the lifetime is precise.
Expand All @@ -206,7 +206,7 @@ The likelihood profile for the double exponential fit::

.. image:: profile2.png

For the double exopnential fit, the profiles look more skewed. The values of :math:`\chi^2` at the minimum are lower, which indicates a better fit.
For the double exopnential fit, the profiles look more skewed. The values of :math:`\chi^2` at the minimum are lower, which indicates a better fit.
However, the binding lifetime labeled 'lifetime 1' never crosses the horizontal line, which indicates that it does not really have an upper bound; this parameter can not be optimized for this data set.
When looking at the likelihood profiles, the single exponential fit is optimal.

Expand All @@ -216,7 +216,7 @@ Bootstrapping
Bootstrapping can be used to select the most suitable model and is a good method for determining the confidence intervals for the fitted parameters.
During bootstrapping, a random sample is taken from the original dataset and fitted. The fitted parameters are gathered in the bootstrapping distribution.
In the example below, we perform 10000 iterations, which means that 10000 times we take a sample from the data and fit the sample with a single exponential distribution.
The resulting 10000 binding lifetimes are plotted in the histogram.
The resulting 10000 binding lifetimes are plotted in the histogram.

Compute and plot the bootstrapping the distribution for the single exponential fit. This will take a while...::

Expand All @@ -239,15 +239,15 @@ Compute and plot the bootstrapping the distribution for the double exponential f

.. image:: bootstrap2.png

The bootstrapping distribution for the double exponential fit is sometimes bimodal and the component :math:`a_2` has a peak close to zero.
The bootstrapping distribution for the double exponential fit is sometimes bimodal and the component :math:`a_2` has a peak close to zero.
This indicates that for many of the bootstrap samples, the fraction associated with the second lifetime was really small and that the parameters of this second lifetime cannot be estimated reliably from the data.

According to the bootstrapping distributions, the single exponential fit is better suitable for the data.

Bayesian Information Criterion
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Typically adding more parameters (components) to a model, will make the fit better. However, having too many parameters can lead to overfitting.
Typically adding more parameters (components) to a model, will make the fit better. However, having too many parameters can lead to overfitting.
The Bayesian information Criterion (BIC) quantifies the quality of the fit by looking at the value of the likelihood function and penalizes the addition of parameters.

The BIC for the single and double exponential fit are respectively given by::
Expand All @@ -267,20 +267,20 @@ We fitted a single exponential and double exponential to the distribution of bin
The likelihood profile and bootstrapping indicated that when using a two-component model, we cannot reliably estimate the second lifetime nor the fraction of events that have this lifetime associated with them.
The BIC indicated that a double exponential is more suitable, but the difference between the small and large model is not very large.

Looking at Figure with the :ref:`double exponential fit <double_exponential_fit>`, there are only a few data points larger than 20 seconds that support the second exponential time scale.
Looking at Figure with the :ref:`double exponential fit <double_exponential_fit>`, there are only a few data points larger than 20 seconds that support the second exponential time scale.
Therefore, the data set is likely too small to support a second exponential time scale. (Fitting two exponentials without overfitting, typically requires a few hundred data points.)

With the current dataset, we conclude that the most suitable model is a single exponential as it gives us the most precise estimates. The fitted lifetime is :math:`\tau = 4` seconds with a 95% confidence interval of (3,5.2) seconds as determined by bootstrapping.
With the current dataset, we conclude that the most suitable model is a single exponential as it gives us the most precise estimates. The fitted lifetime is :math:`\tau = 4` seconds with a 95% confidence interval of (3,5.2) seconds as determined by bootstrapping.
However, given that we do see a hint that there may be a second lifetime involved, it would be worthwhile to gather more data in this case.

Splitting tracks by position
^^^^^^^^^^^^^^^^^^^^^^^^^^^^

When the target sites of the protein are known, the binding lifetimes can also be split by position and analyzed separately [1]_.
For example, to select all tracks from `kymo1_selection` that have an average position larger than 8 micron, type::
For example, to select all tracks from `kymo1_selection` that have an average position larger than 8 micron, type::

track_selection = tracks1[[np.mean(track.position) > 8 for track in tracks1]]

Similarly, we can have a two-sided interval. For example, tracks with a position between 5.5 and 6.2 micron can be obtained by::

track_selection = tracks1[[5.5 < np.mean(track.position) < 6.2 for track in tracks1]]
Expand Down
4 changes: 2 additions & 2 deletions docs/examples/binding_lifetime/bootstrap1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/examples/binding_lifetime/bootstrap2.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/examples/binding_lifetime/double_exponential_fit.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/examples/binding_lifetime/exponential_fit.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
1 change: 1 addition & 0 deletions docs/run_notebooks.py
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,7 @@ def run_notebooks(include_list, reset_cache, only_copy):
exclude_list = [
"nbwidgets", # Exclude the notebook widgets since those require interaction
"cas9_kymotracking",
"binding_lifetime",
]
base_dir = pathlib.Path(__file__).parent.parent.resolve()
nb_test_dir = base_dir / "nb_test"
Expand Down
4 changes: 4 additions & 0 deletions docs/theory/force_calibration/active.rst
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,10 @@
Active Calibration
------------------

.. only:: html

:nbexport:`Download this page as a Jupyter notebook <self>`

For certain applications, passive force calibration, as described above, is not sufficiently accurate.
Using active calibration, the accuracy of the calibration can be improved.
The reason for this is that active calibration uses fewer assumptions than passive calibration.
Expand Down
4 changes: 2 additions & 2 deletions docs/theory/force_calibration/figures/diffusion_spectra.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/theory/force_calibration/figures/hydro_fast.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/theory/force_calibration/figures/lorentzians.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 4 additions & 0 deletions docs/theory/force_calibration/passive.rst
Original file line number Diff line number Diff line change
@@ -1,6 +1,10 @@
Passive calibration
-------------------

.. only:: html

:nbexport:`Download this page as a Jupyter notebook <self>`

Passive calibration is also often referred to as thermal calibration and involves calibration without
moving the trap or stage. In passive calibration, the Brownian motion of the bead in the trap is
analyzed in order to find calibration factors for both the positional detection as well as the force.
Expand Down
4 changes: 2 additions & 2 deletions docs/tutorial/figures/population_dynamics/discrete_pop.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/tutorial/force_calibration/active_calibration.rst
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ Let's compare the active calibration result to passive calibration::
... num_points_per_block=200
... )
>>> print(passive_fit.stiffness)
0.11763849764570819
0.11751264110743381

This value is quite close to that obtained with active calibration above.

Expand All @@ -120,7 +120,7 @@ calculated from the physical input parameters)::
... num_points_per_block=200
... )
>>> print(passive_fit.stiffness)
0.08616565751377737
0.0860734724588009

.. _bead_bead_tutorial:

Expand Down
Loading

0 comments on commit 78a58f4

Please sign in to comment.