fix inclusion of PRO in secondary structure #5065
Conversation
- fix #4913 - ported fix from PyDSSP 0.9.1 by @ShintaroMinami to analysis.dssp.DSSP (see also ShintaroMinami/PyDSSP#2) - new kwarg ignore_proline_donor=True for DSSP (the new default changes the behavior and implements the fix, False recovers old behavior); the kwarg also exists in PyDSSP - updated docs - minimal regression tests - updated CHANGELOG
orbeckst
force-pushed
the
update-dssp-proline-fix
branch
from
b38e9db to
005100f
Compare
| self._donor_mask: Optional[np.ndarray] = ( | ||
| ag.residues.resnames != "PRO" if ignore_proline_donor else None | ||
| ) |
There was a problem hiding this comment.
This may not be correct. The code runs ... but I am not sure if I should be masking corresponding atoms.
| 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). |
There was a problem hiding this comment.
These docs should be more specific and state the shape. I just quickly guessed the shape from https://github.com/ShintaroMinami/PyDSSP/blob/e251a43ff8622fe0a555313b1567edce45e789e8/scripts/pydssp#L30
donor_mask = sequence != 'PRO' if args.ignore_proline_donor else None| if donor_mask is not None | ||
| else np.ones(n_atoms, dtype=float) | ||
| ) | ||
| donor_mask = np.tile(donor_mask[:, np.newaxis], (1, n_atoms)) |
There was a problem hiding this comment.
Is the donor_mask (one element for each residue) really correct for this tiling????
| # 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 | ||
| hbond_map *= local_mask | ||
| hbond_map *= donor_mask |
There was a problem hiding this comment.
Is this correct? The original code uses https://github.com/ShintaroMinami/PyDSSP/blob/e251a43ff8622fe0a555313b1567edce45e789e8/pydssp/pydssp_numpy.py#L72
hbond_map = hbond_map * repeat(donor_mask, 'l1 l2 -> b l1 l2', b=b)(with einops.repeat()). Note that we create our donor_mask with tile so it may already be the right size and shape.
| @@ -13,15 +13,52 @@ | |||
| "pdb_filename", glob.glob(f"{DSSP_FOLDER}/?????.pdb.gz") | |||
| ) | |||
| def test_file_guess_hydrogens(pdb_filename, client_DSSP): | |||
| # run 2.9.0 tests (which include PRO) | |||
| # ignore_proline_donor=False | |||
| # TODO: update reference data for ignore_proline_donor=True | |||
There was a problem hiding this comment.
We should really have correct reference data. About half of the files do not show a difference between ignore_proline_donor=True and ignore_proline_donor=False.
| protein = mda.Universe(TPR, XTC).select_atoms("protein") | ||
| run = DSSP(protein).run(**client_DSSP, stop=10) |
There was a problem hiding this comment.
These lines seemed superfluous as the atomgroup approach is tested separately.
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## develop #5065 +/- ##
========================================
Coverage 93.63% 93.63%
========================================
Files 177 177
Lines 22033 22037 +4
Branches 3115 3115
========================================
+ Hits 20631 20635 +4
Misses 948 948
Partials 454 454 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
Fixes #4913
Changes made in this Pull Request:
PR Checklist
package/CHANGELOGfile updated?package/AUTHORS? (If it is not, add it!)Developers Certificate of Origin
I certify that I can submit this code contribution as described in the Developer Certificate of Origin, under the MDAnalysis LICENSE.
📚 Documentation preview 📚: https://mdanalysis--5065.org.readthedocs.build/en/5065/