Source code for lipyds.analysis.apl

"""
Lipid APL
=========

Classes
-------

.. autofunction:: lipid_area

.. autoclass:: AreaPerLipid
    :members:

"""

from collections import defaultdict
import itertools
from typing import Optional, Union

import numpy as np
from MDAnalysis.analysis.distances import capped_distance
from MDAnalysis.core.universe import Universe
from MDAnalysis.core.groups import AtomGroup
from numpy.typing import ArrayLike

from .base import LeafletAnalysisBase
from ..leafletfinder.utils import get_centers_by_residue
from ..lib.cutils import unwrap_around


[docs]def lipid_area(headgroup_coordinate: ArrayLike, neighbor_coordinates: ArrayLike, other_coordinates: Optional[ArrayLike]=None, box: Optional[ArrayLike]=None, plot: bool=False) -> float: """ Calculate the area of a lipid by projecting it onto a plane with neighboring coordinates and creating a Voronoi diagram. Parameters ---------- headgroup_coordinate: numpy.ndarray Coordinate array of shape (3,) or (n, 3) of the central lipid neighbor_coordinates: numpy.ndarray Coordinate array of shape (n, 3) of neighboring lipids to the central lipid. These coordinates are used to construct the plane of best fit. other_coordinates: numpy.ndarray (optional) Coordinate array of shape (n, 3) of neighboring atoms to the central lipid. These coordinates are *not* used to construct the plane of best fit, but are projected onto it. box: numpy.ndarray (optional) Box of minimum cell, used for unwrapping coordinates. plot: bool (optional) Whether to plot the resulting Voronoi diagram. Returns ------- area: float Area of the central lipid Raises ------ ValueError If a Voronoi cell cannot be constructed for the central lipid, usually because too few neighboring lipids have been given. """ from scipy.spatial import Voronoi # preprocess coordinates headgroup_coordinate = np.asarray(headgroup_coordinate) if len(headgroup_coordinate.shape) > 1: if box is not None: headgroup_coordinates = unwrap_around(headgroup_coordinate.copy(), headgroup_coordinate[0], box) headgroup_coordinate = headgroup_coordinates.mean(axis=0) if len(neighbor_coordinates) < 2: return np.nan if box is not None: neighbor_coordinates = unwrap_around(neighbor_coordinates.copy(), headgroup_coordinate, box) if other_coordinates is not None: other_coordinates = np.asarray(other_coordinates).copy() other_coordinates = unwrap_around(other_coordinates, headgroup_coordinate, box) points = np.concatenate([[headgroup_coordinate], neighbor_coordinates]) points -= headgroup_coordinate center = points.mean(axis=0) points -= center Mt_M = np.matmul(points.T, points) u, s, vh = np.linalg.linalg.svd(Mt_M) # project onto plane if other_coordinates is not None: points = np.r_[points, other_coordinates-center] xy = np.matmul(points, vh[:2].T) xy -= xy[0] # voronoi vor = Voronoi(xy) headgroup_cell_int = vor.point_region[0] headgroup_cell = vor.regions[headgroup_cell_int] if plot: from scipy.spatial import voronoi_plot_2d import matplotlib.pyplot as plt fig = voronoi_plot_2d(vor, show_vertices=False, line_alpha=0.6) plt.plot([vor.points[0][0]], [vor.points[0][1]], 'r+', markersize=12) plt.show() if not all(vertex != -1 for vertex in headgroup_cell): raise ValueError("headgroup not bounded by Voronoi cell points: " f"{headgroup_cell}. " "Try including more neighbor points") # x and y should be ordered clockwise x, y = np.array([vor.vertices[x] for x in headgroup_cell]).T area = np.dot(x[:-1], y[1:]) - np.dot(y[:-1], x[1:]) area += (x[-1] * y[0] - y[-1] * x[0]) lipid_area = 0.5 * np.abs(area) return lipid_area
[docs]class AreaPerLipid(LeafletAnalysisBase): """ Calculate the area of each lipid by projecting it onto a plane with neighboring coordinates and creating a Voronoi diagram. Parameters ---------- universe: AtomGroup or Universe :class:`~MDAnalysis.core.universe.Universe` or :class:`~MDAnalysis.core.groups.AtomGroup` to operate on. select: str (optional) A :meth:`Universe.select_atoms` selection string for atoms that define the lipid head groups, e.g. "name PO4" or "name P*" select_other: str (optional) A :meth:`Universe.select_atoms` selection string for atoms that should be incorporated in the area calculation but that you do not want to calculat areas for. cutoff: float (optional) Cutoff distance (ångström) to look for neighbors cutoff_other: float (optional) Cutoff distance (ångström) to look for neighbors in the ``other`` selection. This is generally shorter than ``cutoff`` -- e.g. you may look for only lipid headgroups in ``select``, but all protein atoms in ``select_other``. max_neighbors: int (optional) Maximum neighbors to consider for area calculation. Past 40 the area tends to converge, and the lower the number of neighbors, the more efficient the calculation. **kwargs Passed to :class:`~lipyds.analysis.base.LeafletAnalysisBase` """ def __init__(self, universe: Union[AtomGroup, Universe], cutoff: float=50, cutoff_other: Optional[float]=15, select_other: Optional[str]="protein", max_neighbors: int=100, **kwargs): super().__init__(universe, **kwargs) if select_other is None: self.other = self.universe.atoms[[]] else: self.other = universe.select_atoms(select_other) - self.residues.atoms if len(self.other): self._get_other_coordinates = self._get_other else: self._get_other_coordinates = lambda x, y, z: None self.cutoff = cutoff if cutoff_other is None: cutoff_other = cutoff self.cutoff_other = cutoff_other self.unique_ids = np.unique(self.ids) self.max_neighbors = max_neighbors def _prepare(self): self.areas = np.zeros((self.n_frames, self.n_residues)) self.areas_by_leaflet = np.zeros((self.n_frames, self.n_leaflets, self.n_residues)) self.areas_by_leaflet[:] = np.nan self.areas_by_attr = [] for i in range(self.n_leaflets): dct = {} for each in self.unique_ids: dct[each] = [] self.areas_by_attr.append(dct) def _get_other(self, central_xyz: ArrayLike, other: ArrayLike, box: Optional[ArrayLike]=None) -> Optional[ArrayLike]: pairs2 = capped_distance(central_xyz, other, self.cutoff_other, box=box, return_distances=False) if len(pairs2): other_xyz = other[np.unique(pairs2[:, 1])] else: other_xyz = None return other_xyz def _single_frame(self): other = self.other.positions box = self.get_box() for lf_i in range(self.n_leaflets): ag = self.leaflet_atomgroups[lf_i] coords = get_centers_by_residue(ag) pairs, dists = capped_distance(coords, coords, self.cutoff, box=box, return_distances=True) if not len(pairs): continue # think this is faster than constructing massive voronoi diagram splix = np.where(np.ediff1d(pairs[:, 0]))[0] + 1 plist = np.split(pairs, splix) dlist = np.split(dists, splix) d_order = [np.argsort(x) for x in dlist] plist = [p[x[:self.max_neighbors+1]] for p, x in zip(plist, d_order)] for pairs_ in plist: central_i = pairs_[0, 0] central_coord = coords[central_i] neighbor_coords = coords[pairs_[1:, 1]] other_coords = self._get_other_coordinates(central_coord, other, box) try: area = lipid_area(central_coord, neighbor_coords, other_coordinates=other_coords, box=box) except ValueError: area = np.nan resindex = ag.residues[central_i].resindex residue_index = self._resindex_to_analysis_order[resindex] residue_label = self.ids[residue_index] self.areas[self._frame_index][residue_index] = area self.areas_by_leaflet[self._frame_index][lf_i][residue_index] = area self.areas_by_attr[lf_i][residue_label].append(area)