fix inclusion of PRO in secondary structure (#5065)

* fix inclusion of PRO in secondary structure

- fix #4913
- Fixes incorrect assignment of secondary structure to proline residues in
   DSSP by porting upstream PyDSSP 0.9.1 fix
- ported fix from PyDSSP 0.9.1 by @ShintaroMinami to analysis.dssp.DSSP
  (see also ShintaroMinami/PyDSSP#2)
- updated dssp test files to match PyDSSP output starting 0.9.1
- refactored `get_hbond_map` to correctly reflect residues and not atoms
- updated docs
- minimal regression tests
- updated CHANGELOG

---------

Co-authored-by: Egor Marin <[email protected]>

Author: orbeckst

package/CHANGELOG

@@ -14,13 +14,13 @@ The rules for this file:
 
 
 -------------------------------------------------------------------------------
-??/??/?? IAlibay
+??/??/?? IAlibay, orbeckst, marinegor
 
  * 2.11.0
 
 Fixes
-
-Enhancements
+ * Fixes incorrect assignment of secondary structure to proline residues in
+   DSSP by porting upstream PyDSSP 0.9.1 fix (Issue #4913)
 
 Changes
 

package/MDAnalysis/analysis/dssp/dssp.py

@@ -112,27 +112,27 @@
    :inherited-members:
 
    .. attribute:: results.dssp
-   
-      Contains the time series of the DSSP assignment as a 
+
+      Contains the time series of the DSSP assignment as a
       :class:`numpy.ndarray` array of shape ``(n_frames, n_residues)`` where each row
-      contains the assigned secondary structure character for each residue (whose 
+      contains the assigned secondary structure character for each residue (whose
       corresponding resid is stored in :attr:`results.resids`). The three characters
       are ['H', 'E', '-'] and representi alpha-helix, sheet and loop, respectively.
 
    .. attribute:: results.dssp_ndarray
-   
+
       Contains the one-hot encoding of the time series of the DSSP assignment
-      as a :class:`numpy.ndarray` Boolean array of shape ``(n_frames, n_residues, 3)`` 
+      as a :class:`numpy.ndarray` Boolean array of shape ``(n_frames, n_residues, 3)``
       where for each residue the encoding is stored as ``(3,)`` shape
-      :class:`numpy.ndarray` of Booleans so that ``True`` at index 0 represents loop 
-      ('-'), ``True`` at index 1 represents helix ('H'), and ``True`` at index 2 
+      :class:`numpy.ndarray` of Booleans so that ``True`` at index 0 represents loop
+      ('-'), ``True`` at index 1 represents helix ('H'), and ``True`` at index 2
       represents sheet 'E'.
 
       .. SeeAlso:: :func:`translate`
-      
+
 
    .. attribute:: results.resids
-   
+
       A :class:`numpy.ndarray` of length ``n_residues`` that contains the residue IDs
       (resids) for the protein residues that were assigned a secondary structure.
 
@@ -144,12 +144,14 @@
 .. autofunction:: translate
 """
 
-from typing import Union
+from typing import Optional, Union
+
 import numpy as np
-from MDAnalysis import Universe, AtomGroup
 
+from MDAnalysis import AtomGroup, Universe
+
+from ...due import Doi, due
 from ..base import AnalysisBase, ResultsGroup
-from ...due import due, Doi
 
 due.cite(
     Doi("10.1002/bip.360221211"),
@@ -163,17 +165,17 @@
 
 try:  # pragma: no cover
     from pydssp.pydssp_numpy import (
-        assign,
         _get_hydrogen_atom_position,
+        assign,
     )
 
     HAS_PYDSSP = True
 
 except ModuleNotFoundError:
     HAS_PYDSSP = False
     from .pydssp_numpy import (
-        assign,
         _get_hydrogen_atom_position,
+        assign,
     )
 
 
@@ -280,6 +282,11 @@ class DSSP(AnalysisBase):
        Enabled **parallel execution** with the ``multiprocessing`` and ``dask``
        backends; use the new method :meth:`get_supported_backends` to see all
        supported backends.
+
+    .. versionchanged:: 2.10.0
+       Change treatment of proline and follow pydssp 0.9.1 (prolines are now explicitly
+       forbidden to participate in the hydrogen bond network). Previous version could yield
+       wrong assignment of prolines.
     """
 
     _analysis_algorithm_is_parallelizable = True
@@ -315,6 +322,7 @@ def __init__(
             ]
             for t in heavyatom_names
         }
+        self._donor_mask: Optional[np.ndarray] = ag.residues.resnames != "PRO"
         self._hydrogens: list["AtomGroup"] = [
             res.atoms.select_atoms(f"name {hydrogen_name}")
             for res in ag.residues
@@ -391,7 +399,7 @@ def _get_coords(self) -> np.ndarray:
 
     def _single_frame(self):
         coords = self._get_coords()
-        dssp = assign(coords)
+        dssp = assign(coords, donor_mask=self._donor_mask)
         self.results.dssp_ndarray.append(dssp)
 
     def _conclude(self):

package/MDAnalysis/analysis/dssp/pydssp_numpy.py

@@ -74,20 +74,20 @@ def _unfold(a: np.ndarray, window: int, axis: int):
 
 
 def _get_hydrogen_atom_position(coord: np.ndarray) -> np.ndarray:
-    """Fills in hydrogen atoms positions if they are abscent, under the
+    """Fills in hydrogen atoms positions if they are absent, under the
     assumption that C-N-H and H-N-CA angles are perfect 120 degrees,
     and N-H bond length is 1.01 A.
 
     Parameters
     ----------
     coord : np.ndarray
-        input coordinates in Angstrom, shape (n_atoms, 4, 3),
+        input coordinates in Angstrom, shape (n_residues, 4, 3),
         where second axes corresponds to (N, CA, C, O) atom coordinates
 
     Returns
     -------
     np.ndarray
-        coordinates of additional hydrogens, shape (n_atoms-1, 3)
+        coordinates of additional hydrogens, shape (n_residues-1, 3)
 
     .. versionadded:: 2.8.0
     """
@@ -118,6 +118,7 @@ def _get_hydrogen_atom_position(coord: np.ndarray) -> np.ndarray:
 
 def get_hbond_map(
     coord: np.ndarray,
+    donor_mask: np.ndarray = None,
     cutoff: float = DEFAULT_CUTOFF,
     margin: float = DEFAULT_MARGIN,
     return_e: bool = False,
@@ -128,8 +129,15 @@ def get_hbond_map(
     ----------
     coord : np.ndarray
         input coordinates in either (n, 4, 3) or (n, 5, 3) shape
-        (without or with hydrogens). If hydrogens are not present, then
-        ideal positions (see :func:_get_hydrogen_atom_positions) are used.
+        (without or with hydrogens respectively), where ``n`` is number of residues.
+        If hydrogens are not present, then ideal positions (see :func:_get_hydrogen_atom_positions)
+        are used.
+    donor_mask : np.array
+         Mask out any hydrogens that should not be considered (in particular HN
+         in PRO). If ``None`` then all H will be used (behavior up to 2.10.0).
+
+         .. versionadded:: 2.10.0
+
     cutoff : float, optional
         cutoff, by default DEFAULT_CUTOFF
     margin : float, optional
@@ -144,8 +152,12 @@ def get_hbond_map(
 
 
     .. versionadded:: 2.8.0
+
+    .. versionchanged:: 2.10.0
+       Support masking of hydrogen donors via `donor_mask` (especially needed
+       for ignoring HN on proline residues). Backport of PRO fix from pydssp 0.9.1.
     """
-    n_atoms, n_atom_types, _ = coord.shape
+    n_residues, n_atom_types, _xyz = coord.shape
     assert n_atom_types in (
         4,
         5,
@@ -161,13 +173,13 @@ def get_hbond_map(
             "Number of atoms should be 4 (N,CA,C,O) or 5 (N,CA,C,O,H)"
         )
     # after this:
-    # h.shape == (n_atoms, 3)
-    # coord.shape == (n_atoms, 4, 3)
+    # h.shape == (n_residues, 3)
+    # coord.shape == (n_residues, 4, 3)
 
     # distance matrix
     n_1, c_0, o_0 = coord[1:, 0], coord[0:-1, 2], coord[0:-1, 3]
 
-    n = n_atoms - 1
+    n = n_residues - 1
     cmap = np.tile(c_0, (n, 1, 1))
     omap = np.tile(o_0, (n, 1, 1))
     nmap = np.tile(n_1, (1, 1, n)).reshape(n, n, 3)
@@ -191,18 +203,32 @@ def get_hbond_map(
         return e
 
     # mask for local pairs (i,i), (i,i+1), (i,i+2)
-    local_mask = ~np.eye(n_atoms, dtype=bool)
-    local_mask *= ~np.diag(np.ones(n_atoms - 1, dtype=bool), k=-1)
-    local_mask *= ~np.diag(np.ones(n_atoms - 2, dtype=bool), k=-2)
+    local_mask = ~np.eye(n_residues, dtype=bool)
+    local_mask *= ~np.diag(np.ones(n_residues - 1, dtype=bool), k=-1)
+    local_mask *= ~np.diag(np.ones(n_residues - 2, dtype=bool), k=-2)
+    # mask for donor H absence (Proline)
+    donor_mask = (
+        np.array(donor_mask).astype(float)
+        if donor_mask is not None
+        else np.ones(n_residues, dtype=float)
+    )
+    donor_mask = np.tile(donor_mask[:, np.newaxis], (1, n_residues))
     # hydrogen bond map (continuous value extension of original definition)
     hbond_map = np.clip(cutoff - margin - e, a_min=-margin, a_max=margin)
     hbond_map = (np.sin(hbond_map / margin * np.pi / 2) + 1.0) / 2
-    hbond_map = hbond_map * local_mask
+
+    assert hbond_map.shape == local_mask.shape == donor_mask.shape
+
+    hbond_map *= local_mask
+    hbond_map *= donor_mask
 
     return hbond_map
 
 
-def assign(coord: np.ndarray) -> np.ndarray:
+def assign(
+    coord: np.ndarray,
+    donor_mask: np.ndarray = None,
+) -> np.ndarray:
     """Assigns secondary structure for a given coordinate array,
     either with or without assigned hydrogens
 
@@ -214,6 +240,12 @@ def assign(coord: np.ndarray) -> np.ndarray:
         (N, CA, C, O) atoms coordinates (if k=4), or (N, CA, C, O, H) coordinates
         (when k=5).
 
+    donor_mask : np.array
+         Mask out any hydrogens that should not be considered (in particular HN
+         in PRO). If ``None`` then all H will be used (behavior up to 2.9.0).
+
+         .. versionadded:: 2.10.0
+
     Returns
     -------
     np.ndarray
@@ -222,9 +254,13 @@ def assign(coord: np.ndarray) -> np.ndarray:
 
 
     .. versionadded:: 2.8.0
+
+    .. versionchanged:: 2.10.0
+       Support masking of donors.
+
     """
     # get hydrogen bond map
-    hbmap = get_hbond_map(coord)
+    hbmap = get_hbond_map(coord, donor_mask=donor_mask)
     hbmap = np.swapaxes(hbmap, -1, -2)  # convert into "i:C=O, j:N-H" form
 
     # identify turn 3, 4, 5

testsuite/MDAnalysisTests/analysis/test_dssp.py

@@ -1,8 +1,9 @@
-import pytest
 import glob
-import MDAnalysis as mda
 
+import MDAnalysis as mda
+import pytest
 from MDAnalysis.analysis.dssp import DSSP, translate
+
 from MDAnalysisTests.datafiles import DSSP as DSSP_FOLDER
 from MDAnalysisTests.datafiles import TPR, XTC
 
@@ -32,8 +33,6 @@ def test_trajectory(client_DSSP):
     assert (
         first_frame[:10] != last_frame[:10] == avg_frame[:10] == "-EEEEEE---"
     )
-    protein = mda.Universe(TPR, XTC).select_atoms("protein")
-    run = DSSP(protein).run(**client_DSSP, stop=10)
 
 
 def test_atomgroup(client_DSSP):

testsuite/MDAnalysisTests/data/dssp/1eteA.pdb.dssp

@@ -1 +1 @@
---------------HHHHHHHHH-------EEEEE--E------HHHHHHHHHHHHHHHHHH---HHHHHHHHHHHHHHHHHHH-----------EEEEEHHHHHHHHHHHHHHHH--------HHH----E-- 1eteA.pdb
+--------------HHHHHHHHH-------EEEEE--E------HHHHHHHHHHHHHHHHHH---HHHHHHHHHHHHHHHHHHH-----------EEEEEHHHHHHHHHHHHHHH---------HHH----E-- 1eteA.pdb

testsuite/MDAnalysisTests/data/dssp/2fvvA.pdb.dssp

@@ -1 +1 @@
-----E-----EEEEEEEEE------EEEEEE-------EEE-EEE------HHHHHHHHHHHHH-EEEEE--E-EEEEE----EEEEEEEEEE-EE----HHHHHH---EEEEEHHHHHHHH----HHHHH---- 2fvvA.pdb
+----E-----EEEEEEEEE------EEEEEE-------EE--EEE------HHHHHHHHHHHHH-EEEEE--E-EEEEE----EEEEEEEEEE-EE----HHHHHH---EEEEEHHHHHHHH----HHHHH---- 2fvvA.pdb

testsuite/MDAnalysisTests/data/dssp/2j49A.pdb.dssp

@@ -1 +1 @@
----HHHHHHHHHHHHH---HHHHHHHHHHHHHHHHHHHHHHHHH-HHHHHHHHHHHHHHHHHHHHHHHH-------HHHHHH-HHHHHHH---EEEEE-HHHHHHHHHHHH--HHHHHHHHHHHHHHHEEEEE- 2j49A.pdb
+---HHHHHHHHHHHHH------HHHHHHHHHHHHHHHHHHHHHH-HHHHHHHHHHHHHHHHHHHHHHHH-------HHHHHH-HHHHHHH---EEEEE-HHHHHHHHHHHH--HHHHHHHHHHHHHHHEEEEE- 2j49A.pdb

testsuite/MDAnalysisTests/data/dssp/3a4rA.pdb.dssp

@@ -1 +1 @@
--------EEEEEE------EEEEEE----EHHHHHHHHHHHH-------EEEE--EE-----EHHHH------EEEEE- 3a4rA.pdb
+-------EEEEEE------EEEEEE-----HHHHHHHHHHHH-------EEEE--EE------HHHH------EEEEE- 3a4rA.pdb

testsuite/MDAnalysisTests/data/dssp/3aqgA.pdb.dssp

@@ -1 +1 @@
--EE-----EEEE--------E-EEEEEE----E-EEEEEE---E---EE----EEEEEE------E-EEEEEE-----EEEEEE----EEEEE----EEEEEE-----EEE-EEEEEEE--EE-EEEEEEEE- 3aqgA.pdb
+-EE-----EEEE--------E-EEEEEE----E-EEEEEE---E---EE----EEEEEE------E-EEEEEE-----EEEEEE----EEEEE----EEEEE-------EE-EEEEEEE--EE-EEEEEEE-- 3aqgA.pdb

testsuite/MDAnalysisTests/data/dssp/3e8mA.pdb.dssp

@@ -1 +1 @@
------EEEE-----------EE-----EE-----HHHHHHHHHHH---EEEEE-----HHHHHHHHH----EEE-----HHHHHHHHHHHH---HHHEEEE---HHHHHHH----EEE------HHHH--------------HHHHHHHHH----HHHHHHH-- 3e8mA.pdb
+-----EEEE-----------EE-----EE-----HHHHHHHHHHH----EEEE-----HHHHHHHHH----EEE-----HHHHHHHHHHHH---HHHEEEE---HHHHHHH----EEE------HHHH--------------HHHHHHHHH----HHHHHHH-- 3e8mA.pdb