added symm site

python/Symmetric_Site.py

@@ -0,0 +1,351 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Nov 14 16:53:18 2024.
+
+@author: js2746
+"""
+from Site import Site
+import numpy as np
+from utils import calculate_hist_mode, calculate_hist_mean, calculate_dG
+
+
+class Symmetric_Site:
+    """
+    An aggregation of multiple binding sites on/in a protein/inclusion. User \
+    defines the base_site Site object first (including setting the bin_coords!)\
+    and then provides it to the Symmetric_Site constructor.
+
+    Attributes
+    ----------
+    temperature : float
+        The temperature of your system in K.
+
+    Calculated Properties
+    ---------------------
+    symmetry : int
+        The N-fold symmetry desired. I.E. 5 would yield 5 Sites.
+    site_list : list
+        The list of constituent Site objects that make up this Symmetric_Site.
+    site_counts_histogram : numpy ndarray
+        One-dimensional ndarray where the histogrammed ligand bead counts are \
+        stored. e.g. [12, 5, 0, 0, 1, 0] would correspond to 12 frames having \
+        zero beads in the Site, 5 frames having one bead in the Site, 0 frames \
+        having 2, 3, or 5 beads in the site, and 1 frame having 4 beads in the \
+        Site.
+    bulk_counts_histogram : numpy ndarray
+        One-dimensional ndarray where the histogrammed ligand bead counts are \
+        stored. e.g. [12, 5, 0, 0, 1, 0] would correspond to 12 frames having \
+        zero beads in the bulk patch, 5 frames having one bead in the patch, 0 \
+        frames having 2, 3, or 5 beads in the patch, and 1 frame having 4 beads\
+        in the patch.
+    n_peak : int
+        The mode of the bulk histogram. Indicates the cut-off for P_unocc.
+    dG : float
+        The binding affinity of the lipid for the Site, in kcal/mol.
+    dG_std : float
+        The standard deviation of the mean binding affinity for the Sites that\
+        comprise this Symmetric_Site.
+    """
+
+    def __init__(self, symmetry, base_site, Ntheta):
+        """
+        Create a Symmetric_Site by creating clones of the base_site and \
+        rotating them symmetrically around the origin.
+
+        Parameters
+        ----------
+        symmetry : int
+            The N-fold symmetry desired. I.E. 5 would yield 5 Sites.
+        base_site : Site
+            The original Site object that should be cloned and rotated.
+
+        """
+        assert isinstance(symmetry, int), "symmetry must be an integer."
+        assert Ntheta % symmetry == 0, "This symmetry does not evenly divide \
+            across the number of theta bins."
+        assert isinstance(base_site, Site), "base_site must be a Site."
+        assert base_site.bin_coords is not None, "The base_site needs to be fully defined before creating a Symmetric_Site."
+        self._symmetry = symmetry
+        self._site_list = self._make_symmetric_sites(base_site, Ntheta)
+        assert len(self.site_list) == symmetry, "Number of Sites does not match symmetry."
+        self.temperature = base_site.temperature
+
+    def __iter__(self):
+        """Iterate through the site_list."""
+        for site in self.site_list:
+            yield site
+
+    @property
+    def symmetry(self):
+        """
+        Tell me the symmetry, but don't let me change the symmetry.
+
+        Returns
+        -------
+        int
+            The N-fold symmetry of the Symmetric_Site.
+
+        """
+        return self._symmetry
+
+    @property
+    def site_list(self):
+        """
+        Tell me the site_list, but don't let me change the site_list.
+
+        Returns
+        -------
+        list
+            List of constituent Site objects that comprise this Symmetric_Site.
+
+        """
+        return self._site_list
+
+    @property
+    def site_counts_histogram(self):
+        """
+        Tell me the current counts, in histogram form, for the Symmetric_Site.
+
+        Returns
+        -------
+        site_counts_histogram : numpy ndarray
+            One-dimensional ndarray where the histogrammed ligand bead counts \
+            are stored. e.g. [12, 5, 0, 0, 1, 0] would correspond to 12 frames \
+            having zero beads in the Site, 5 frames having one bead in the \
+            Site, 0 frames having 2, 3, or 5 beads in the site, and 1 frame \
+            having 4 beads in the Site.
+
+        """
+        return _aggregate_site_counts_histograms(self.site_list)
+
+    @property
+    def bulk_counts_histogram(self):
+        """
+        Tell me the current counts, in histogram form, for the Symmetric_Site. \
+        In practice, this is just the bulk_counts_histogram for the base_site.
+
+        Returns
+        -------
+        bulk_counts_histogram : numpy ndarray
+            One-dimensional ndarray where the histogrammed ligand bead counts \
+            are stored. e.g. [12, 5, 0, 0, 1, 0] would correspond to 12 frames \
+            having zero beads in the bulk patch, 5 frames having one bead in \
+            the patch, 0 frames having 2, 3, or 5 beads in the patch, and 1 \
+            frame having 4 beads in the patch.
+
+        """
+        return _check_bulk_counts_histogram(self.site_list)
+
+    @property
+    def n_peak(self):
+        """
+        Tell me what the n_peak is.
+
+        Returns
+        -------
+        int
+            The mode of the bulk distribution in a patch of membrane that has \
+            equal accessible area to the site.
+
+        """
+        return calculate_hist_mode(self.bulk_counts_histogram)
+
+    @property
+    def dG(self):
+        """
+        Calculate the binding affinity of the lipid for this Symmetric_Site, \
+        including the bulk correction factor dG_ref.
+
+        Returns
+        -------
+        float
+            The total binding affinity, in kcal/mol.
+
+        """
+        n_peak = self.n_peak
+        dG_site = calculate_dG(self.site_counts_histogram, n_peak, self.temperature)
+        dG_ref = calculate_dG(self.bulk_counts_histogram, n_peak, self.temperature)
+        return dG_site - dG_ref
+
+    @property
+    def dG_std(self):
+        """
+        Calculate the standard deviation of the delta G values across the \
+        constituent Sites that comprise this Symmetric_Site.
+
+        Returns
+        -------
+        float
+            The standard deviation
+
+        """
+        dGs = []
+        for site in self.site_list:
+            dGs.append(site.dG)
+        return np.std(dGs)
+
+    def update_counts_histogram(self, bulk, counts_data):
+        """
+        Update the counts histograms for all constituent Sites.
+
+        Parameters
+        ----------
+        bulk : boolean
+            If True, update the counts histogram for the bulk patch. If False,\
+            update the counts histogram for the site.
+        counts_data : ndarray
+            If bulk=True, provide 1D nddarray containing bulk counts. \
+            If bulk=False, provide the 3D ndarray containing binned counts.
+
+        Returns
+        -------
+        None.
+
+        """
+        for site in self.site_list:
+            site.update_counts_histogram(bulk, counts_data)
+
+    def predict_accessible_area(self, bulk_area, mode=True):
+        """
+        Predict the accessible area of the site. A reasonable method is to \
+        multiply the area of the bulk patch you just analyzed by the ratio of\
+        the means (or modes) for the site distribution and the bulk \
+        distribution. This will put you in the ballpark of the bulk patch area.
+
+        Parameters
+        ----------
+        bulk_area : float
+            The area of the bulk patch previously analyzed in square Angstroms.
+        mode : boolean
+            If True, use the site and bulk modes rather than the means. Default\
+            is True. If False use means instead of modes (untested feature).
+
+        Returns
+        -------
+        predicted_accessible_area : float
+            The area of the bulk patch you should analyze next to try and more\
+            closely match the site distribution. Units are in square Angstroms.
+
+        """
+        if mode:
+            site = calculate_hist_mode(self.site_counts_histogram)
+            bulk = calculate_hist_mode(self.bulk_counts_histogram)
+        else:
+            site = calculate_hist_mean(self.site_counts_histogram)
+            bulk = calculate_hist_mean(self.bulk_counts_histogram)
+        predicted_accessible_area = bulk_area * (site / bulk)
+        return predicted_accessible_area
+
+    def _make_symmetric_sites(self, base_site, Ntheta):
+        """
+        Create identical sites to the base_site, rotated symmetrically around \
+        the origin.
+
+        Parameters
+        ----------
+        base_site : Site
+            The Site object that you want to replicate symmetrically.
+        Ntheta : int
+            The total number of theta bins in the circle.
+
+        Returns
+        -------
+        site_list : list of Sites
+            The list of all Sites that comprise this Symmetric_Site.
+
+        """
+        base_site.name = base_site.name + '_1'
+        site_list = [base_site]
+        for site_number in range(1, self.symmetry):
+            site_name = base_site.name + '_' + str(site_number + 1)
+            new_site = Site(site_name, base_site.leaflet, base_site.temperature)
+            new_site.bin_coords = self._rotate_bin_coords(base_site.bin_coords, Ntheta, site_number)
+            site_list.append(new_site)
+        return site_list
+
+    def _rotate_bin_coords(self, bin_coords, Ntheta, site_number):
+        """
+        Rotate the provided bin_coords around the circle.
+
+        Parameters
+        ----------
+        bin_coords : list of tuples
+            The bins that belong to this site in (r, theta) format. e.g. \
+            [(2, 10), (2, 11), (2, 12)] would correspond to the 11th, 12th, and \
+            13th theta bins in the 3rd radial bin from the origin. Bin coordinates \
+            are zero-indexed by convention.
+        Ntheta : int
+            The number of theta bins in the circle.
+        site_number : int
+            Which constituent site this is.
+
+        Returns
+        -------
+        rotated_bin_coords : list of tuples
+            Should match input bin_coords in length and first tuple component; \
+            second tuple components should all be shifted (rotated).
+
+        """
+        rotated_bin_coords = []
+        for each_bin in bin_coords:
+            r_bin, theta_bin = each_bin
+            shift = Ntheta // self.symmetry
+            rotated_theta_bin = theta_bin + shift * site_number
+            if rotated_theta_bin >= Ntheta:
+                rotated_theta_bin -= Ntheta
+            rotated_bin_coords.append((r_bin, rotated_theta_bin))
+        return rotated_bin_coords
+
+
+def _aggregate_site_counts_histograms(site_list):
+    """
+    Cycle through all the sites and add their counts histograms together.
+
+    Parameters
+    ----------
+    site_list : list
+        List of Sites.
+
+    Returns
+    -------
+    counts : ndarray
+        1D numpy array of histogrammed bead counts.
+
+    """
+    hist_lengths = []
+    for site in site_list:
+        hist_length = site.site_counts_histogram.shape[0]
+        hist_lengths.append(hist_length)
+    max_len = max(hist_lengths)
+    counts = np.zeros(max_len)
+    for site in site_list:
+        counts_to_add = site.site_counts_histogram.copy()
+        if counts_to_add.shape[0] < max_len:
+            padding = max_len - counts_to_add.shape[0]
+            counts_to_add = np.pad(counts_to_add, (0, padding), mode='constant', constant_values=0)
+        counts += counts_to_add
+    return counts
+
+
+def _check_bulk_counts_histogram(site_list):
+    """
+    Cycle through each Site and make sure the bulk_counts_histograms all match.\
+    Return one of them.
+
+    Parameters
+    ----------
+    site_list : list
+        List of Sites.
+
+    Returns
+    -------
+    bulk : ndarray
+        1D numpy array of histogrammed bead counts from the bulk distribution.
+
+    """
+    first_site = site_list[0]
+    bulk = first_site.bulk_counts_histogram.copy()
+    for site in site_list[1:]:
+        assert bulk.all() == site.bulk_counts_histogram.all(), "One or more sites have different bulk histograms. This shouldn't be possible."
+    return bulk

python/plotting.py

@@ -12,7 +12,7 @@
 import matplotlib.transforms as mtransforms
 from matplotlib.colors import ListedColormap, Normalize
 from Site import Site
-# from Symmetric_Site import Symmetric_Site
+from Symmetric_Site import Symmetric_Site
 from utils import calculate_hist_mode
 
 
@@ -58,10 +58,10 @@ def outline_site_new(ax, site, grid_dims):
     if isinstance(site, Site):
         for each_bin in site.bin_coords:
             ax = outline_bin(ax, each_bin, grid_dims)
-    # elif isinstance(site, Symmetric_Site):
-    #     for each_site in site.site_list:
-    #         for each_bin in each_site.bin_coords:
-    #             ax = outline_bin(ax, each_bin, grid_dims)
+    elif isinstance(site, Symmetric_Site):
+        for each_site in site.site_list:
+            for each_bin in each_site.bin_coords:
+                ax = outline_bin(ax, each_bin, grid_dims)
     else:
         Exception("site must be a Site or Symmetric_Site.")
     return ax