From ca6122cd8c22470203420903b1c784bebbdd703d Mon Sep 17 00:00:00 2001 From: AnthonyBogetti Date: Mon, 1 Dec 2025 16:01:12 -0500 Subject: [PATCH 1/2] fixes to new-hinge branch --- prody/dynamics/analysis.py | 142 +++++++++++++++++++++++- prody/tests/dynamics/test_get_hinges.py | 97 ++++++++++++++++ 2 files changed, 237 insertions(+), 2 deletions(-) create mode 100644 prody/tests/dynamics/test_get_hinges.py diff --git a/prody/dynamics/analysis.py b/prody/dynamics/analysis.py index 987359a76..af26a5841 100644 --- a/prody/dynamics/analysis.py +++ b/prody/dynamics/analysis.py @@ -21,9 +21,9 @@ __all__ = ['calcCollectivity', 'calcCovariance', 'calcCrossCorr', 'calcFractVariance', 'calcSqFlucts', 'calcRMSFlucts', 'calcMostMobileNodes', 'calcTempFactors', - 'calcProjection', 'calcCrossProjection', + 'calcProjection', 'calcCrossProjection', 'calcSpecDimension', 'calcPairDeformationDist', - 'calcDistFlucts', 'calcHinges', 'calcHitTime', + 'calcDistFlucts', 'calcHinges', 'getHinges', 'calcHitTime', 'calcAnisousFromModel', 'calcScipionScore', 'calcHemnmaScore'] #'calcEntropyTransfer', 'calcOverallNetEntropyTransfer'] @@ -646,6 +646,144 @@ def identify(v): return sorted(set(hinge_list)) return hinges + +def getHinges(gnm, n_modes=None, threshold=15, space=None, trim=False): + """ + [HZ384] H Zhang, M Gur, I Bahar (2024) Global hinge sites of proteins as target sites for drug binding Proc Natl Acad Sci USA 121 (49), e2414333121 + + Updated hinge identification based on [HZ384]. This will: + 1) If GNM model is provided without specification of n_modes, number of modes will be determined to achieve 33% cumulative variance. + If selected GNM modes are provided directly, all modes will be used for hinge detection. + 2) Identify hinge regions based on crossovers and residues within band, whose magnitude is specified by the threshold. + 3) Merge overlapping or adjacent hinge regions. + 4) Reduce transient hinge regions. + 5) Trim hinges at N- or C-terminal ends. + 6) Return list of hinges by mode. + + :param modes: GNM model or ModeSet + :type modes: GNM or Mode instance + + :param n_modes: Number of modes to consider. Defaults to all. + :type n_modes: int or None + + :param threshold: Threshold controlling band width for hinge region identification. + :type threshold: float + + :param space: spacing between hinges. Defaults to None. Higher space value will minimize the number of local hinges while retaining global hinges. + :type space: int or None + + :param trim: Whether to hinges at the N- or C-terminal ends by Protein length // 20. Defaults to False. + :type trim: bool + + :return: List of sorted hinge indices + :rtype: [list[int], list[int], ...] + + """ + + if isinstance(gnm, ModeSet) is False: + + if isinstance(gnm, GNMBase) is False: + raise TypeError("Input must be a GNM model or ModeSet instance.") + else: + if gnm._kirchhoff is None: + raise ValueError("Kirchhoff matrix not built. Build Kirchhoff matrix before calculating modes.") + + if n_modes is None: + ### Default to auto-select based on 33% cumulative variance ### + cumin =0.33 + n_modes = 1 + gnm.calcModes(n_modes='all') + fv = calcFractVariance(gnm) + cumulative_variances = np.cumsum(fv) + n_modes = np.argmax(cumulative_variances >= cumin) + 1 + cuvar = cumulative_variances[n_modes - 1] + LOGGER.info(f"Auto-selected {n_modes} modes for {cumin} cumulative variance.") + + else: + gnm.calcModes(n_modes=n_modes) + LOGGER.info(f"Using {n_modes} modes.") + + vecs = gnm[:n_modes].getEigvecs() + else: + ### Use all provided modes for Hinge detection ### + vecs = gnm.getArray() + + def detectHinges(v, threshold): + """Detect hinge-like regions within single eigenvector.""" + + crx = np.nonzero(np.diff(np.sign(v)))[0] + band = np.sqrt(1 / len(v)) / threshold + regs = [] + + for i in crx: + r = [i, i + 1] + + if abs(v[i]) > band and abs(v[i + 1]) > band: + regs.append(r) + else: + j = i - 1 + while j >= 0 and abs(v[j]) < band: + r.insert(0, j) + j -= 1 + j = i + 2 + while j < len(v) and abs(v[j]) < band: + r.append(j) + j += 1 + regs.append(r) + + ### Merge overlapping or adjacent regions (separated by space value) ### + + regs = sorted(regs, key=lambda x: x[0]) + merged = [regs[0]] + s = 1 + space if space is not None else 0 + + for curr in regs[1:]: + prev = merged[-1] + if curr[0] <= prev[-1] + s: + merged[-1] = list(range(prev[0], max(prev[-1], curr[-1]) + 1)) + else: + merged.append(curr) + + ### Select hinges from regions ### + fil=[] + p, _ = vecs.shape + + for reg in merged: + if len(reg) <= p // 10: + #if v[reg[0] - 1] * v[reg[-1] + 1] > 0: + if v[reg[0] - 1] * v[reg[-1]] > 0: + continue + + if len(reg) >= s + 5: + fil.append(reg[0]) + fil.append(reg[-1]) + + if len(reg) <= 5: + av = [abs(v[i]) for i in reg] + smv = np.argsort(av)[:1] + fil.extend(reg[i] for i in smv) + + return [int(x) for x in fil] + + ### Detect hinges for each mode ### + n_hinges = [] + p, m = vecs.shape + + + for i in range(m): + hinges = detectHinges(vecs[:, i], threshold) + n_hinges.append(hinges) + + if trim == False: + return n_hinges + else: + to_trim = p // 20 + n_hinges = [[h for h in hinges if h >= to_trim and h < (p - to_trim)] for hinges in n_hinges] + + return n_hinges + + + def calcHitTime(model, method='standard'): """Returns the hit and commute times between pairs of nodes calculated based on a :class:`.NMA` object. diff --git a/prody/tests/dynamics/test_get_hinges.py b/prody/tests/dynamics/test_get_hinges.py new file mode 100644 index 000000000..fb1c26bb2 --- /dev/null +++ b/prody/tests/dynamics/test_get_hinges.py @@ -0,0 +1,97 @@ +import unittest +import numpy as np +from prody import * + +# TODO: Import your function here +# from your_module_name import getHinges + +# --- THE UNIT TEST SUITE --- +class TestHingesWithRealGNM(unittest.TestCase): + + @classmethod + def setUpClass(cls): + """ + Runs once before all tests. + Downloads 1AKE and calculates modes to save time. + """ + print("\n[Setup] Downloading 1AKE and building GNM...") + + # 1. Parse PDB (Downloads if not found) + cls.pdb = parsePDB("1AKE", subset='ca') + + # 2. Build GNM + cls.gnm = GNM('1AKE') + cls.gnm.buildKirchhoff(cls.pdb) + cls.gnm.calcModes(n_modes='all') + + print("[Setup] GNM built successfully.") + + def test_output_structure(self): + """Verify the output is a list of lists of integers.""" + # Use explicit n_modes to keep it simple + hinges = getHinges(self.gnm, n_modes=2) + + self.assertIsInstance(hinges, list, "Output must be a list") + self.assertEqual(len(hinges), 2, "Must return results for exactly 2 modes") + self.assertIsInstance(hinges[0], list, "Each mode result must be a list") + + # Check if contents are integers (if any hinges found) + if len(hinges[0]) > 0: + self.assertIsInstance(hinges[0][0], (int, np.integer), "Hinge indices must be integers") + + def test_auto_mode_selection(self): + """Verify that passing n_modes=None triggers auto-selection.""" + # The function defaults to 33% cumulative variance + hinges = getHinges(self.gnm, n_modes=None) + + num_modes_selected = len(hinges) + + # Calculate expected number manually + fv = calcFractVariance(self.gnm) + cum_var = np.cumsum(fv) + expected_modes = np.argmax(cum_var >= 0.33) + 1 + + print(f"\n[Test Auto] Auto-selected {num_modes_selected} modes.") + self.assertEqual(num_modes_selected, expected_modes, + f"Should select {expected_modes} modes for 33% variance") + + def test_trim_functionality(self): + """Verify that 'trim=True' removes hinges at the very ends.""" + # 1AKE has ~214 residues. + # trim = length // 20 = 214 // 20 = 10 residues cut from each end. + + # First, run WITHOUT trim + hinges_untrimmed = getHinges(self.gnm, n_modes=1, trim=False)[0] + + # Then, run WITH trim + hinges_trimmed = getHinges(self.gnm, n_modes=1, trim=True)[0] + + # Calculate the forbidden zones + n_atoms = self.gnm.numAtoms() + trim_zone = n_atoms // 20 + + # Ensure no hinges in trimmed zones + for h in hinges_trimmed: + self.assertTrue(trim_zone <= h < (n_atoms - trim_zone), + f"Hinge at {h} should have been trimmed (Limit: {trim_zone} to {n_atoms-trim_zone})") + + # Verify we didn't lose Valid hinges from the middle + # (This assumes the un-trimmed version had middle hinges, which 1AKE Mode 1 usually does) + self.assertGreater(len(hinges_trimmed), 0, "1AKE Mode 1 should have valid central hinges") + + def test_threshold_sensitivity(self): + """Verify that a lower threshold finds MORE (or wider) hinges.""" + # Low threshold = wider band = stricter condition to be considered a hinge + # Actually, in your logic: + # band = sqrt(1/N) / threshold + # Higher threshold -> smaller band -> Easier to be OUTSIDE the band -> Hinge region is narrower/more specific? + # Let's test equality. Changing threshold MUST change output. + + h_std = getHinges(self.gnm, n_modes=1, threshold=15)[0] + h_extreme = getHinges(self.gnm, n_modes=1, threshold=1.0)[0] + + self.assertNotEqual(h_std, h_extreme, "Changing threshold should affect hinge detection results") + +if __name__ == '__main__': + # Verbosity=2 shows individual test status (OK/FAIL) + unittest.main(verbosity=2) From 3deb753ac1247a625f439eca3228b0e20e40d046 Mon Sep 17 00:00:00 2001 From: AnthonyBogetti Date: Mon, 1 Dec 2025 16:06:16 -0500 Subject: [PATCH 2/2] clean up code --- prody/dynamics/analysis.py | 1 - 1 file changed, 1 deletion(-) diff --git a/prody/dynamics/analysis.py b/prody/dynamics/analysis.py index af26a5841..b46aceb34 100644 --- a/prody/dynamics/analysis.py +++ b/prody/dynamics/analysis.py @@ -750,7 +750,6 @@ def detectHinges(v, threshold): for reg in merged: if len(reg) <= p // 10: - #if v[reg[0] - 1] * v[reg[-1] + 1] > 0: if v[reg[0] - 1] * v[reg[-1]] > 0: continue