diff --git a/cinnabar/femap.py b/cinnabar/femap.py index e7b13783..77dffe6f 100644 --- a/cinnabar/femap.py +++ b/cinnabar/femap.py @@ -550,6 +550,16 @@ def to_legacy_graph(self) -> nx.DiGraph: return g def draw_graph(self, title: str = "", filename: Union[str, None] = None): + """ + Draw the graph using matplotlib. + + Parameters + ---------- + title : str, optional + Title for the graph. + filename : str or None, optional + If provided, the graph will be saved to this file. If None, the graph will be displayed. + """ plt.figure(figsize=(10, 10)) graph = self.to_legacy_graph() diff --git a/cinnabar/stats.py b/cinnabar/stats.py index 8c414fdf..292f7a08 100644 --- a/cinnabar/stats.py +++ b/cinnabar/stats.py @@ -148,8 +148,6 @@ def mle(graph: nx.DiGraph, factor: str = "f_ij", node_factor: Union[str, None] = The number 'factor' is the node attribute on which the MLE will be calculated, where d'factor' will be used as the standard error of the factor - We assume the free energy of node 0 is zero. - Reference : https://pubs.acs.org/doi/abs/10.1021/acs.jcim.9b00528 Xu, Huafeng. "Optimal measurement network of pairwise differences." Journal of Chemical Information and Modeling 59.11 (2019): 4720-4728. @@ -172,8 +170,7 @@ def mle(graph: nx.DiGraph, factor: str = "f_ij", node_factor: Union[str, None] = Returns ------- f_i : np.array with shape (n_ligands,) - f_i[i] is the absolute free energy of ligand i in kT - f_i[0] = 0 + f_i[i] is the absolute free energy of ligand i in kcal/mol C : np.array with shape (n_ligands, n_ligands) C[i,j] is the covariance of the free energy estimates of i and j diff --git a/docs/concepts/estimators.rst b/docs/concepts/estimators.rst new file mode 100644 index 00000000..a42efaac --- /dev/null +++ b/docs/concepts/estimators.rst @@ -0,0 +1,72 @@ +.. _estimators: + +=============================== +Absolute Free Energy Estimators +=============================== + +Relative free energy calculations produce :math:`\Delta\Delta G` (differences between two ligands). To compare these with experiment, +or to rank ligands by affinity, we need absolute free energies (:math:`\Delta G`) for each ligand. + +This requires an **estimator**: a method that takes the network of relative free energies and produces absolute values. + + +Maximum Likelihood Estimation (MLE) +----------------------------------- + +The Maximum Likelihood Estimation (MLE) [1]_ method is the **default** estimator used in cinnabar to obtain absolute +free energies (:math:`\Delta G`) from a network of relative free energies (:math:`\Delta\Delta G`). + + +The Core Idea +~~~~~~~~~~~~~~ +To place every ligand on a common absolute scale, we need to find a set of :math:`\Delta G` values that best explain all +relative differences (:math:`\Delta\Delta G`) simultaneously. The MLE method does this by asking: + + What set of :math:`\Delta G` values makes the observed data most likely, given the reported uncertainties? + +This framing naturally integrates all edges and cycles in the graph simultaneously. + +The Likelihood Function +~~~~~~~~~~~~~~~~~~~~~~~ + +Suppose we have a network with two ligands ``i`` and ``j``, with observed relative free energy :math:`\Delta\Delta G_{ij}` and +uncertainty :math:`\sigma_{ij}`. The model assumes each measurement is normally distributed: + +.. math:: + + \Delta\Delta G_{ij} \approx \mathcal{N}(\Delta G_j - \Delta G_i, \sigma_{ij}^2) + +The **likelihood** is the product of probabilities for all edges in the graph. The MLE procedure finds the set of :math:`\Delta G` +values that maximises this likelihood (or equivalently, minimises the negative log-likelihood). + +Uncertainty Propagation +~~~~~~~~~~~~~~~~~~~~~~~ + +Input uncertainties (:math:`\sigma_{ij}`) are explicitly included in the likelihood function. This means more precise edges +(smaller uncertainty) have greater weight in determining the solution. However, this does mean that high confidence but low +accuracy edges can impact the entire network and so robust uncertainty estimates on input data are crucial. + + +Centering of Results +~~~~~~~~~~~~~~~~~~~~ + +The absolute :math:`\Delta G` scale is arbitrary: adding a constant to all :math:`\Delta G` values does not change any relative differences :math:`\Delta\Delta G`. +As a result, the MLE solution is typically centred around zero (or another chosen reference). To compare with experimental +values, an experimental shift must be applied. By default ``cinnabar`` will align the mean of predicted and +experimental :math:`\Delta G` in the plotting functions. + + +Limitations +~~~~~~~~~~~ + +- The MLE method **can not** use multiple independent measurements of the same edge to improve precision automatically. + Each edge must be represented by a single :math:`\Delta\Delta G` and uncertainty. If multiple measurements are available, + they should be combined externally (e.g. via weighted averaging) before input to the estimator. + + + +References +~~~~~~~~~~~ + +.. [1] Xu, H., 2019. Optimal measurement network of pairwise differences. Journal of Chemical Information and Modeling, 59(11), pp.4720-4728. + diff --git a/docs/concepts/femap.rst b/docs/concepts/femap.rst new file mode 100644 index 00000000..ac2f5b78 --- /dev/null +++ b/docs/concepts/femap.rst @@ -0,0 +1,43 @@ +============================ +The ``FEMap`` Data Structure +============================ + +The :class:`.FEMap` is the **core abstraction in cinnabar**. +It provides a unified representation of free energy data, connecting the inputs from relative binding free energy +calculations with the **analyses** and **visualizations** that cinnabar enables. + +Why a Graph? +------------ + +Relative free energy calculations are inherently **relational**: they compare two ligands at a time, producing a free +energy difference (:math:`\Delta\Delta G`) with an associated uncertainty. If we want to reason about a whole series of ligands, which +may be connected by multiple pairwise comparisons, we need to connect these edges in a logical way. + +A **graph** is a natural way to represent this data, where: + +- **Nodes** represent individual ligands. +- **Absolute values** (:math:`\Delta G` experimental or calculated) can be atached to nodes as attributes. +- **Edges** represent pairwise free energy differences (:math:`\Delta\Delta G`) between ligands. + +This graph representation is powerful: it allows integration of relative and absolute data, and it provides the +foundation for robust statistical analysis and visualization. + + +From Relative to Absolute +~~~~~~~~~~~~~~~~~~~~~~~~~ + +Although we calculate relative free energies directly, many applications require absolute binding free energies for +each ligand in order to rank the ligands for prioritization. + +The :class:`.FEMap` provides convenience methods to support this transformation via a :ref:`maximum likelihood estimation (MLE) ` [1]_ +method by default. These methods take advantage of the entire network of relative free energies to infer absolute values. +Thus they require a graph which is at least weakly connected, with at least one path between any two ligands in the graph. +The ``FEMap`` provides utilities to check (:meth:`.FEMap.check_weakly_connected`) and visualise (:meth:`.FEMap.draw_graph`) the connectivity of the network. + + +References +~~~~~~~~~~~ + +.. [1] Xu, H., 2019. Optimal measurement network of pairwise differences. Journal of Chemical Information and Modeling, 59(11), pp.4720-4728. + + diff --git a/docs/concepts/index.rst b/docs/concepts/index.rst new file mode 100644 index 00000000..2cc5feb3 --- /dev/null +++ b/docs/concepts/index.rst @@ -0,0 +1,13 @@ +Core Concepts +============= + +This section covers the core concepts of the software, offering a deep dive into the key abstractions, data structures and theory that underpin its functionality. +For detailed step-by-step tutorials on using the software, please refer to the tutorials section. + + +.. toctree:: + :maxdepth: 1 + + femap + estimators + plotting diff --git a/docs/concepts/plotting.rst b/docs/concepts/plotting.rst new file mode 100644 index 00000000..f9141df2 --- /dev/null +++ b/docs/concepts/plotting.rst @@ -0,0 +1,68 @@ +======================= +Plotting Best Practices +======================= + +Visualisation plays a crucial role in assessing the accuracy of binding free energy calculations. The same +raw data can be represented in many different ways, and the choice of plot or analysis metric determines which aspects of performance are +highlighted. This page explains the plot types available in ``cinnabar``, and the best practices they follow. For more +detailed information we recommend reading the companion article [1]_. + +Edgewise :math:`\Delta\Delta G` Plots +------------------------------------- + +Edgewise relative free energy plots allow direct comparison of calculated and experimental :math:`\Delta\Delta G` values. +These plots are best for method developers, as they show how well individual transformations +are predicted and highlight outliers. These plots can be generated from an :class:`.FEMap` using the :meth:`plot_DDGs ` method. + +Best practices: + +- **Only error statistics (RMSE, MUE)** are shown by default. Correlation measures (:math:`R^{2}`, :math:`\rho`) are not meaningful, as the direction of a relative transformation is arbitrary. +- **Error bars** are shown on to represent uncertainty in both calculated and experimental values. +- **Statistics uncertainty estimates** are bootstrapped (1000 samples with replacement) to provide 95% confidence intervals. + + +Absolute :math:`\Delta G` Plots +------------------------------- + +Absolute free energy plots compare calculated and experimental :math:`\Delta G` values for each ligand. Absolute +:math:`\Delta G` predictions can be obtained from a connected network of relative free energies via an :ref:`estimator `. +These plots are useful for both method developers and users, as they show how well the overall ranking of ligands is predicted. +As they depend on the entire network they accumulate errors from multipule aspects of the calculation, including +edge accuracy, uncertainty quantification, perturbation network design and estimator, giving a complete picture of protocol performance. +These plots can be generated from an :class:`.FEMap` using the :meth:`plot_DGs ` method. + +Best practices: + +- **Both error (RMSE, MUE) and correlation** (:math:`R^{2}`, :math:`\rho`) statistics are shown by default, as absolute data are directional and correlation is meaningful. +- **Error bars** are shown on to represent uncertainty in both calculated and experimental values. +- **Statistics uncertainty estimates** are bootstrapped (1000 samples with replacement) to provide 95% confidence intervals. +- **Mean centering** is applied by default to align the mean of calculated and experimental values. This is needed when analyzing the outcomes of relaltive free energy simulations as the absolute scale is arbitrary. + +Pairwise (all-to-all) :math:`\Delta\Delta G` Plots +-------------------------------------------------- + +Pairwise or all-to-all relative free energy plots compare all possible pairwise :math:`\Delta\Delta G` values between ligands in a dataset. +These are generated by calculating pairwise differences from the estimated absolute :math:`\Delta G` values. +These plots are especially useful for method developers because they remove analysis biases introduced by the perturbation network design: +every method is compared on the same set of pairwise values, even if the underlying networks differ. These plots can be generated from an :class:`.FEMap` using the :meth:`plot_all_DDGs ` method. + +Best practices: + +- Follow the same guidelines as for edgewise plots (error statistics only, with uncertainties represented). +- Use these plots for fair, network-independent comparisons between methods. + + +Summary of Best Practices +------------------------- + +- Use :math:`\Delta\Delta G` plots for edge-level diagnostics. +- Use absolute :math:`\Delta G` plots for global performance and ranking. +- Always represent uncertainties, both at the individual estimate level and on reported statistics (e.g. RMSE). +- When comparing across methods, use all-to-all pairwise :math:`\Delta\Delta G` to enable fair comparisons. + + +References +~~~~~~~~~~~ + +.. [1] Hahn, D.F., Bayly, C.I., Boby, M.L., Macdonald, H.E.B., Chodera, J.D., Gapsys, V., Mey, A.S., Mobley, D.L., Benito, L.P., Schindler, C.E. and Tresadern, G., 2022. Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [article v1. 0]. Living journal of computational molecular science, 4(1), p.1497. + diff --git a/docs/getting_started.rst b/docs/getting_started.rst deleted file mode 100644 index adff8ba5..00000000 --- a/docs/getting_started.rst +++ /dev/null @@ -1,4 +0,0 @@ -Getting Started -=============== - -``cinnabar`` is a package to aid with the consistent and reliable plotting of results from free energy calculations. Fair assessment and comparison of computational methods relies on robust and reproduceable metrics, which `freeenergy framework` aims to provide. diff --git a/docs/index.rst b/docs/index.rst index 81a3369a..d304234d 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -30,7 +30,7 @@ Contents :maxdepth: 1 installation - getting_started + concepts/index Reference