1
$\begingroup$

I know the algorithm for creating a Hydrogen atom and adding to a residue:

Point3d  create_hydrogen(Point3d C, Point3d N, Point3d CA, Point3d H) 
{
  H.set(N);
  H -= C;
  H.norm();

  Point3d tmp2(N);
  tmp2 -= CA;
  tmp2.norm();

  H += tmp2;
  H.norm(N_H_BOND_LENGTH);
  H += N;

  return H;
}

Point3d  create_hydrogen(Residue prev_residue, Residue residue, Point3d H) 
{
  return create_hydrogen(prev_residue.find_atom(" C  "), 
                              residue.find_atom(" N  "),
                              residue.find_atom(" CA "), 
                                   H);
}

Atom create_hydrogen(Residue prev_residue, Residue residue) 
{
  Atom H = make_shared_atom(0, " H  ", 1);
  create_hydrogen(prev_residue, residue, H);
  residue.add(H);
  return H;
}

I didn't find this implementation in BioPython.

Is this already implemented in BioPython? If YES, in which module?

If NO, what would be the implementation like in BioPython?

$\endgroup$
19
  • $\begingroup$ Search for Bio.PDB.Structure.Structure object in biopython code $\endgroup$
    – pippo1980
    Nov 1, 2021 at 6:57
  • $\begingroup$ I would suggest again using already existing tools to add hydrogens. You said your PI had concerns with this, but I must say that it is not a random process like letting kids adding baubles to a Christmas tree: for sp3 heavy atoms bound by two heavy atoms and for terminal heavy atom in sp2 the hydrogens have a defined geometric position. Terminal amines, hydroxyls and methyls will require more advanced calculations, but bar for ST turns these will not contribute to SS. $\endgroup$ Nov 1, 2021 at 9:38
  • $\begingroup$ @MatteoFerla, I understand. I talked to him, and he doesn't agree. This is his project and his decision is final. :) What can I do? $\endgroup$
    – user366312
    Nov 1, 2021 at 10:08
  • $\begingroup$ @pippo1980, Search for Bio.PDB.Structure.Structure object in biopython code --- there is nothing useful in that module. $\endgroup$
    – user366312
    Nov 1, 2021 at 10:09
  • $\begingroup$ Switch to openbabel ? Look how the Arpeggio guys handled the visualization of ligand protein interaction of pdb : github.com/harryjubb/arpeggio/blob/master/arpeggio.py $\endgroup$
    – pippo1980
    Nov 1, 2021 at 10:59

2 Answers 2

1
$\begingroup$

The other answer is correct. But I thought I'd give a few pointers in how one can understand how to do something with a given Python module —teach a man to fish kind of thing...

In a Jupyter Notebook or IPython shell (don't), you can do:

  • help(function) print out the docstring
  • dir(object) shows all the public and magic attributes and methods —note a class will not show the instance methods
  • type(object) or object.__class__ gives you what the class is, but object.__class__.__mro__ tells you what it inherits.
  • in the standard library module inspect is the function inspect.getsource(function), which for a notebook is messy due to newlines so needs to be print_code = lambda fun: print(inspect.getsource(fun)) and print_code(fun).

In your case you read a file...

from Bio import PDB
parser = PDB.PDBParser()
structure : PDB.Structure.Structure = parser.get_structure("tryp-cage", "1L2Y.first.pdb")

So parser.get_structure "magically" makes a structure, eh? This means it can add stuff! As mentioned print_code(parser.get_structure) will reveal its secrets, or within it call some self._private_fun(), which can be inspected (self = parser) all the way down the rabbit hole. Until you find a nice constructor that declares nicely an atom instance or whatever you ever fancy searching for —that is the great thing of open source code.

In the specific case of PDB module, everything inherits Entity, which has the attribute child_list, which is what the iterators (__iter__ is the topmost) read without making a copy. The end of the above rabbit hole may just add a child Entity-inheriting class instance to this list, but it also do other things that is worth keeping an eye out.

In the case of residue.add(atom), creating the Bio.PDB.Atom.Atom instance is slightly problematic as there's no typehinting or :type xxx: declarations in the docstring. But doing type(atom.bfactor) etc. on the attributes of an existing atom will reveal how to do it. One the serial number is a PDB thing, which is sequential for the model.

last_serial = list(structure.get_residues())[-1].child_list[-1].serial_number + 1
residue.add(PDB.Atom.Atom(name=' H  ', coord=np.array([1,1,1]), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number=last_serial+1,element='H'))

The coordinates will not be [1,1,1], but the point that you have calculate. Say for H, this would be coplanar with C, CA and N, that is 1 and smidge Å away from N at an angle not quite in the middle as found in amide bonds.

$\endgroup$
12
  • 1
    $\begingroup$ If the other answer is correct, what do you think about the algorithm I posted? This algorithm is taken from the c++ source code I have in my hand. $\endgroup$
    – user366312
    Nov 1, 2021 at 16:10
  • 1
    $\begingroup$ As you would realize, I can write Python and C++ program, but I am just confused about hydrogens. :) $\endgroup$
    – user366312
    Nov 1, 2021 at 16:11
  • $\begingroup$ Biopython documentation is not that rich especially for ones like me that use google as primary/unique docs search $\endgroup$
    – pippo1980
    Nov 1, 2021 at 16:45
  • $\begingroup$ I have expanded the answer with an example, showing that what I wrote is how one finds the correct solution without documentation. In terms of how to do the maths, you have to do that in numpy as its simple maths —I do not know if the C++ code is 100% correct as it is projecting on the line between C-CA midpoint and N, which I am not sure is the correct angle. Totally easy to verify though $\endgroup$ Nov 1, 2021 at 17:04
  • 1
    $\begingroup$ The CB atom is not chiral in valine, but the atom names make it "name chiral" (not a real word, but just like isotopic chirality): looking at valine with the HB atom forewards, you have CG2 on the side of the N (left) and CG1 on the C (right) in a predictable counterclockwise fashion. Isoleucine is chiral on the CB atom and the CG1 is on the left as it is more complex due to CD1. As a result they do not superpose well by atom names, but do if CG2 is mapped to CG1 and CG1: CG2. $\endgroup$ Nov 2, 2021 at 16:49
1
$\begingroup$

Not sure but here https://biopython.org/docs/1.75/api/Bio.PDB.Residue.html Bio.PDB.Residue module Your can find :

add(self, atom) :

      Add an Atom object.

Checks for adding duplicate atoms, and raises a PDBConstructionException if so.

Can’t try it just now but who knows it may work

Googling around I found out that if you plan to add Hydrogens to PDB files [I believe nowadays PDBx is the repository grade standard] you could do thar following PDB format rules. I am copying them from UCSF Chimera website (https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens):


Hydrogen Atoms

In brief, conventions for hydrogen atoms in version 3.0 PDB format are as follows:

Hydrogen atom records follow the records of all other atoms of a particular residue.
A hydrogen atom name starts with H. The next part of the name is based on the name of the connected nonhydrogen atom. For example, in amino acid residues, H is followed by the remoteness indicator (if any) of the connected atom, followed by the branch indicator (if any) of the connected atom; if more than one hydrogen is connected to the same atom, an additional digit is appended so that each hydrogen atom will have a unique name. Hydrogen atoms in standard nucleotides and amino acids (other than the rarely seen HXT) are named according to the IUPAC recommendations (Markley et al., Pure Appl Chem 70:117 (1998)). Names of hydrogen atoms in HETATM residues are determined in a similar fashion.
If the name of a hydrogen has four characters, it is left-justified starting in column 13; if it has fewer than four characters, it is left-justified starting in column 14. 

In the following excerpt from entry 1vm3, atom H is attached to atom N. Atom HA is attached to atom CA; the remoteness indicator A is the same for these atoms. Two hydrogen atoms are connected to CB, one is connected to CG, three are connected to CD1, and three are connected to CD2.

ATOM     10  N   LEU A   2       4.595   6.365   3.756  1.00  0.00           N
ATOM     11  CA  LEU A   2       4.471   5.443   2.633  1.00  0.00           C
ATOM     12  C   LEU A   2       5.841   5.176   2.015  1.00  0.00           C
ATOM     13  O   LEU A   2       6.205   4.029   1.755  1.00  0.00           O
ATOM     14  CB  LEU A   2       3.526   6.037   1.578  1.00  0.00           C
ATOM     15  CG  LEU A   2       2.790   4.919   0.823  1.00  0.00           C
ATOM     16  CD1 LEU A   2       3.803   3.916   0.262  1.00  0.00           C
ATOM     17  CD2 LEU A   2       1.817   4.196   1.769  1.00  0.00           C
ATOM     18  H   LEU A   2       4.169   7.246   3.704  1.00  0.00           H
ATOM     19  HA  LEU A   2       4.063   4.514   2.992  1.00  0.00           H
ATOM     20  HB2 LEU A   2       2.804   6.675   2.065  1.00  0.00           H
ATOM     21  HB3 LEU A   2       4.099   6.623   0.873  1.00  0.00           H
ATOM     22  HG  LEU A   2       2.234   5.353   0.004  1.00  0.00           H
ATOM     23 HD11 LEU A   2       4.648   4.447  -0.148  1.00  0.00           H
ATOM     24 HD12 LEU A   2       3.334   3.331  -0.516  1.00  0.00           H
ATOM     25 HD13 LEU A   2       4.137   3.260   1.052  1.00  0.00           H
ATOM     26 HD21 LEU A   2       0.941   3.892   1.216  1.00  0.00           H
ATOM     27 HD22 LEU A   2       1.522   4.860   2.568  1.00  0.00           H
ATOM     28 HD23 LEU A   2       2.296   3.323   2.188  1.00  0.00           H


```
$\endgroup$
1
  • $\begingroup$ As you can see from the algorithm, adding Hydrogen is not that straightforward. :D If it were that simple, I wouldn't have posted this question!! $\endgroup$
    – user366312
    Nov 1, 2021 at 15:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.