Improve AddHs for molecules read from PDB #1647

mwojcikowski · 2017-11-06T14:16:28Z

What does this implement/fix? Explain your changes.

This PR fixes behavior of AddHs for molecules read from PDB file. Current behavior is to ignore residue information. What's more, the Hs written to PDB file are always HETATM and of UNL residue.

Proposed changes add AtomPDBResidueInfo based on the information collected from the root atom.
I decided against adding a flag, since it doesn't break default behavior and is kind of expected when working with proteins.

bp-kelley · 2017-11-13T18:33:59Z

Code/GraphMol/AddHs.cpp

@@ -386,6 +387,13 @@ void addHs(RWMol &mol, bool explicitOnly, bool addCoords,
      mol.addBond(aidx, newIdx, Bond::SINGLE);
      mol.getAtomWithIdx(newIdx)->updatePropertyCache();
      if (addCoords) setHydrogenCoords(&mol, newIdx, aidx);
+      AtomPDBResidueInfo *info = (AtomPDBResidueInfo *)(newAt->getMonomerInfo());
+      if (info && info->getMonomerType() == AtomMonomerInfo::PDBRESIDUE) {
+        AtomPDBResidueInfo *newInfo = new AtomPDBResidueInfo(" H  ", newIdx, "", info->getResidueName(),


Using the atom index for the serial number might not work in practice as they are not related. I would find the current maximum serial number and add one to it.

It might be nice is the hydrogens have unique names, especially if one is using/making a refinement dictionary. In general, this is done by using the number of hydrogen atoms on the residue, i.e. " H1", " H2" again you will need to count the existing hydrogens and names first.

Note that in practice, the hydrogen index "wraps around" due to a mis reading of the initial PDB spec way back in the day. Hydrogen 123 is written "3H12" for instance. There was no reason for the H to be in column 1, but someone codified that and now we're stuck :)

The names only have to be unique for the residue, not for the protein.

mwojcikowski · 2017-11-16T20:47:46Z

@bp-kelley All requested changes were implemented. I have also refactored the code to be in one place. The hydrogens are numbered one per-resiude basis + wrap-around. Effectively it supports 1099 hydrogens to be uniquely numbered (the 10XY are wrapped around to YH0X). Plus the serial number, as you suggested, is incremented maximum serial number seen among other atoms.

greglandrum · 2017-11-17T05:30:34Z

Code/GraphMol/AddHs.cpp

@@ -409,6 +422,47 @@ void addHs(RWMol &mol, bool explicitOnly, bool addCoords,
    }
    // update the atom's derived properties (valence count, etc.)
    newAt->updatePropertyCache();
+


Not critical, but it seems like this code should be in a separate function (in an anonymous namespace).

mwojcikowski · 2017-11-17T11:20:06Z

@greglandrum Done. I've fixed a bug: the H idx was reset for every atom, now it is for each residue.

mwojcikowski · 2018-01-09T14:35:04Z

Ping @bp-kelley @greglandrum

greglandrum

I like the change and think it's quite useful, but I don't like that it's called for every molecule, even ones that have not been generated from PDB files. That's an unneccessary performance hit. I think it would be better to have a new function in the public API: AddHsWithResidueInfo() that calls AddHs() and then this.

mwojcikowski · 2018-01-15T10:08:26Z

Why not expose AssignHsPDBInfo to the world, then? User could decide should it be called.

greglandrum · 2018-01-15T11:57:53Z

That's more or less what I was suggesting, but if you call it AddHsWithResidueInfo() then you could imagine the same call working for other types of polymers (not just things from PDB files).

greglandrum · 2018-01-15T11:59:09Z

Note that I'm not suggesting that you actually make it general right now, just that we have a stub in place so that we can do that in the future.

mwojcikowski · 2018-01-15T15:23:53Z

Yet another idea - what if we execute AssignHsResidueInfo inside the AddHs function, but optionally (new param AssignResidueInfo=False). This way we retain fast(er) AddHs for non-PDB molecules, give an opportunity to assign the info (user's choice), and finally we don't need to copy constructor/helper/wrapper code connected to AddHs (Java, python, duplicated API which further changes should be also modified).

mwojcikowski · 2018-01-16T11:56:42Z

@greglandrum I've modified the code to my previous post.

greglandrum

Looks good. The one change suggested is to change the argument name.

greglandrum · 2018-01-17T04:16:37Z

Code/GraphMol/Wrap/MolOps.cpp

@@ -931,8 +931,11 @@ struct molops_wrapper {
    - addCoords: (optional) if this toggle is set, The Hs will have 3D coordinates\n\
      set.  Default value is 0 (no 3D coords).\n\
 \n\
-    - onlyOnHs: (optional) if this sequence is provided, only these atoms will be\n\
+    - onlyOnAtoms: (optional) if this sequence is provided, only these atoms will be\n\


Nice catch. Thanks for that fix. :-)

greglandrum · 2018-01-17T04:17:09Z

Code/GraphMol/AddHs.cpp

 void addHs(RWMol &mol, bool explicitOnly, bool addCoords,
-           const UINT_VECT *onlyOnAtoms) {
+           const UINT_VECT *onlyOnAtoms, bool residueInfo) {


I'd prefer the argument to be called addResidueInfo as I think it's more explicit.

mwojcikowski · 2018-01-18T09:06:00Z

Done

greglandrum · 2018-01-18T16:11:11Z

Thanks @mwojcikowski !

@bp-kelley : I'll take care of getting this merged into modern_cxx over the next couple of days.

* Add AtomPDBResidueInfo during molops::addhs * Test addhs and AtomPDBResidueInfo * Little cleanup * Fix serial * Refactor code, add unique Hs labels * Move code to separate function + bugfix. * Make function optional (AddHs residueInfo param) * Rename argument

mwojcikowski added 2 commits November 6, 2017 13:51

Add AtomPDBResidueInfo during molops::addhs

6aa5469

Test addhs and AtomPDBResidueInfo

b93f5eb

mwojcikowski changed the title ~~[MRG] Improve AddHs for molecules read freom PDB~~ [MRG] Improve AddHs for molecules read from PDB Nov 6, 2017

bp-kelley requested changes Nov 13, 2017

View reviewed changes

mwojcikowski added 3 commits November 16, 2017 17:47

Little cleanup

3939b94

Fix serial

752bfa0

Refactor code, add unique Hs labels

df0662e

mwojcikowski changed the title ~~[MRG] Improve AddHs for molecules read from PDB~~ Improve AddHs for molecules read from PDB Nov 16, 2017

greglandrum reviewed Nov 17, 2017

View reviewed changes

Move code to separate function + bugfix.

9674816

mwojcikowski force-pushed the fix-pdb-addhs branch from cfff095 to 9674816 Compare November 17, 2017 15:43

greglandrum mentioned this pull request Jan 10, 2018

Treat bonds in PDB CONECT records explicitly, but make blacklisted ones zero-order. #1658

Merged

greglandrum requested changes Jan 15, 2018

View reviewed changes

Make function optional (AddHs residueInfo param)

ea711e6

mwojcikowski force-pushed the fix-pdb-addhs branch from 37559ed to ea711e6 Compare January 16, 2018 08:32

greglandrum requested changes Jan 17, 2018

View reviewed changes

Rename argument

5cd0055

greglandrum approved these changes Jan 17, 2018

View reviewed changes

greglandrum added the enhancement label Jan 18, 2018

greglandrum added this to the 2018_03_1 milestone Jan 18, 2018

greglandrum merged commit 50299cb into rdkit:master Jan 18, 2018

mwojcikowski deleted the fix-pdb-addhs branch February 6, 2018 13:36

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improve AddHs for molecules read from PDB #1647

Improve AddHs for molecules read from PDB #1647

mwojcikowski commented Nov 6, 2017

bp-kelley Nov 13, 2017

mwojcikowski commented Nov 16, 2017

greglandrum Nov 17, 2017

mwojcikowski commented Nov 17, 2017

mwojcikowski commented Jan 9, 2018

greglandrum left a comment

mwojcikowski commented Jan 15, 2018

greglandrum commented Jan 15, 2018

greglandrum commented Jan 15, 2018

mwojcikowski commented Jan 15, 2018

mwojcikowski commented Jan 16, 2018

greglandrum left a comment

greglandrum Jan 17, 2018

greglandrum Jan 17, 2018

mwojcikowski commented Jan 18, 2018

greglandrum commented Jan 18, 2018

Improve AddHs for molecules read from PDB #1647

Improve AddHs for molecules read from PDB #1647

Conversation

mwojcikowski commented Nov 6, 2017

What does this implement/fix? Explain your changes.

bp-kelley Nov 13, 2017

Choose a reason for hiding this comment

mwojcikowski commented Nov 16, 2017

greglandrum Nov 17, 2017

Choose a reason for hiding this comment

mwojcikowski commented Nov 17, 2017

mwojcikowski commented Jan 9, 2018

greglandrum left a comment

Choose a reason for hiding this comment

mwojcikowski commented Jan 15, 2018

greglandrum commented Jan 15, 2018

greglandrum commented Jan 15, 2018

mwojcikowski commented Jan 15, 2018

mwojcikowski commented Jan 16, 2018

greglandrum left a comment

Choose a reason for hiding this comment

greglandrum Jan 17, 2018

Choose a reason for hiding this comment

greglandrum Jan 17, 2018

Choose a reason for hiding this comment

mwojcikowski commented Jan 18, 2018

greglandrum commented Jan 18, 2018