# Read PDB file, extract dihedral angles, modify dihedral angles, reconstruct Cartesian coordinates, and write PDB file

As the title summarizes, I am trying to:

1. Read a PDB file (for example, 1enh.pdb).
2. Extract the backbone dihedral angles (phi, psi).
3. Modify the dihedral angles (phi, psi) arbitrarily (for example, replace them with np.random.uniform(-np.pi, np.pi, size = (56, 2)) in 1enh, since there are 56 pairs of dihedrals).
4. Reconstruct the Cartesian coordinates of the modified backbone.
5. Write a PDB file based on the modified backbone.

From this question it seems the tools for what I want are available in BioPython. There is also a set_angle() method. But I only have a very limited knowledge of Python and BioPython. I am stuck at Steps 3 and 4, and dubious on Step 5.

This is what I have tried so far:

#!/usr/bin/python

import numpy as np
from Bio.PDB import PDBParser, PICIO, PDBIO

## Step 1

parser = PDBParser()
structure = parser.get_structure("", "1enh.pdb")

## Step 2

structure.atom_to_internal_coordinates()
chain = list(structure.get_chains())[0]
ic_chain = chain.internal_coord
dihedrals = ic_chain.dihedra

## Step 3 - ??

# Modifiy (phi,psi). For example, replace (phi,psi) with
newdihedrals = np.random.uniform(-np.pi, np.pi, size = (56, 2))

## Step 4 - ??

modified.internal_to_atom_coordinates()

## Step 5

io = PDBIO()
io.set_structure(modified)
io.save('1enh_modified.pdb', preserve_atom_numbering = True)


EDIT: I do not care about possible clashes if the side-chain is included in the modified PDB. In any case, writing a modified PDB file that only includes the backbone is also fine.

• POINT 5 READS : Write a PDB file based on the modified backbone. --> Just C-alpha or side-chains too ?? How to avoid clashes in case of the latter ? Aug 20, 2022 at 9:59
• Thanks for the comment. C-alpha would suffice. If side-chains are included, I do not really care about clashes, as I could determine them later and reject the modified structure. Aug 20, 2022 at 10:57
• Any Python package that checks PDB geometry just for clashes ?? Aug 20, 2022 at 11:12
• In PyMOL I remember it is possible to detect clashes, I do not know other functionality beyond it. In any case, my main interest is on tweaking the (phi,psi) angles and reconstruct the protein backbone in a PDB file. Side-chains could be present or not. Thanks. Aug 20, 2022 at 11:50

ok, here my input pdb, example_short.pdb :

ATOM      1  N   MET A   1      14.067  31.664  -1.639  1.00 80.87           N
ATOM      2  HN1 MET A   1      14.250  31.259  -0.699  1.00  0.00           H
ATOM      3  HN2 MET A   1      14.891  32.220  -1.945  1.00  0.00           H
ATOM      4  HN3 MET A   1      13.229  32.279  -1.596  1.00  0.00           H
ATOM      5  CA  MET A   1      13.824  30.564  -2.612  1.00 79.17           C
ATOM      6  HA  MET A   1      14.702  29.921  -2.671  1.00  0.00           H
ATOM      7  C   MET A   1      13.556  31.158  -3.989  1.00 78.52           C
ATOM      8  O   MET A   1      13.301  32.365  -4.085  1.00 79.40           O
ATOM      9  CB  MET A   1      12.634  29.724  -2.148  1.00 83.38           C
ATOM     10  HB1 MET A   1      12.600  29.755  -1.059  1.00  0.00           H
ATOM     11  HB2 MET A   1      11.725  30.170  -2.552  1.00  0.00           H
ATOM     12  CG  MET A   1      12.680  28.264  -2.580  1.00 85.71           C
ATOM     13  HG1 MET A   1      12.066  28.145  -3.473  1.00  0.00           H
ATOM     14  HG2 MET A   1      12.274  27.650  -1.776  1.00  0.00           H
ATOM     15  SD  MET A   1      14.354  27.697  -2.951  1.00100.24           S
ATOM     16  CE  MET A   1      15.137  27.675  -1.338  1.00 97.10           C
ATOM     17  HE1 MET A   1      15.118  28.678  -0.913  1.00  0.00           H
ATOM     18  HE2 MET A   1      14.597  26.991  -0.683  1.00  0.00           H
ATOM     19  HE3 MET A   1      16.170  27.342  -1.441  1.00  0.00           H
ATOM     20  N   GLU A   2      13.604  30.306  -5.018  1.00 75.32           N
ATOM     21  H   GLU A   2      13.628  29.291  -4.794  1.00  0.00           H
ATOM     22  CA  GLU A   2      13.626  30.684  -6.444  1.00 72.66           C
ATOM     23  HA  GLU A   2      13.508  29.732  -6.961  1.00  0.00           H
ATOM     24  C   GLU A   2      14.903  31.384  -6.843  1.00 65.44           C
ATOM     25  O   GLU A   2      15.489  31.127  -7.892  1.00 67.46           O
ATOM     26  CB  GLU A   2      12.516  31.672  -6.810  1.00 72.71           C
ATOM     27  HB1 GLU A   2      12.094  32.068  -5.886  1.00  0.00           H
ATOM     28  HB2 GLU A   2      12.960  32.486  -7.383  1.00  0.00           H
ATOM     29  CG  GLU A   2      11.397  31.081  -7.622  1.00 77.72           C
ATOM     30  HG1 GLU A   2      10.901  31.877  -8.178  1.00  0.00           H
ATOM     31  HG2 GLU A   2      11.809  30.354  -8.321  1.00  0.00           H
ATOM     32  CD  GLU A   2      10.385  30.393  -6.731  1.00 79.89           C
ATOM     33  OE1 GLU A   2       9.171  30.574  -6.987  1.00 82.69           O
ATOM     34  OE2 GLU A   2      10.824  29.702  -5.779  1.00 72.73           O
ATOM     35  N   ASN A   3      15.194  32.399  -6.043  1.00 60.12           N
ATOM     36  H   ASN A   3      14.712  32.444  -5.123  1.00  0.00           H
ATOM     37  CA  ASN A   3      16.140  33.452  -6.363  1.00 56.21           C
ATOM     38  HA  ASN A   3      16.072  33.715  -7.419  1.00  0.00           H
ATOM     39  C   ASN A   3      17.548  32.968  -6.047  1.00 50.40           C
ATOM     40  O   ASN A   3      18.513  33.666  -6.331  1.00 50.95           O
ATOM     41  CB  ASN A   3      15.821  34.667  -5.493  1.00 57.10           C
ATOM     42  HB1 ASN A   3      15.913  34.381  -4.445  1.00  0.00           H
ATOM     43  HB2 ASN A   3      16.538  35.457  -5.717  1.00  0.00           H
ATOM     44  CG  ASN A   3      14.425  35.196  -5.732  1.00 62.07           C
ATOM     45  OD1 ASN A   3      14.220  36.001  -6.641  1.00 69.12           O
ATOM     46  ND2 ASN A   3      13.464  34.748  -4.923  1.00 57.21           N
ATOM     47 1HD2 ASN A   3      12.486  35.079  -5.045  1.00  0.00           H
ATOM     48 2HD2 ASN A   3      13.693  34.068  -4.170  1.00  0.00           H
ATOM     49  N   PHE A   4      17.643  31.764  -5.492  1.00 47.14           N
ATOM     50  H   PHE A   4      16.781  31.193  -5.378  1.00  0.00           H
ATOM     51  CA  PHE A   4      18.914  31.219  -5.038  1.00 46.77           C
ATOM     52  HA  PHE A   4      19.696  31.915  -5.340  1.00  0.00           H
ATOM     53  C   PHE A   4      19.239  29.888  -5.695  1.00 48.73           C
ATOM     54  O   PHE A   4      18.490  28.919  -5.594  1.00 49.50           O
ATOM     55  CB  PHE A   4      18.922  31.095  -3.507  1.00 42.32           C
ATOM     56  HB1 PHE A   4      18.091  30.456  -3.207  1.00  0.00           H
ATOM     57  HB2 PHE A   4      19.862  30.635  -3.201  1.00  0.00           H
ATOM     58  CG  PHE A   4      18.787  32.415  -2.806  1.00 39.42           C
ATOM     59  CD1 PHE A   4      17.545  32.878  -2.408  1.00 40.76           C
ATOM     60  HD1 PHE A   4      16.673  32.232  -2.509  1.00  0.00           H
ATOM     61  CD2 PHE A   4      19.881  33.258  -2.648  1.00 35.13           C
ATOM     62  HD2 PHE A   4      20.857  32.941  -3.015  1.00  0.00           H
ATOM     63  CE1 PHE A   4      17.397  34.138  -1.889  1.00 36.49           C
ATOM     64  HE1 PHE A   4      16.400  34.520  -1.669  1.00  0.00           H
ATOM     65  CE2 PHE A   4      19.757  34.474  -2.043  1.00 32.88           C
ATOM     66  HE2 PHE A   4      20.640  35.090  -1.873  1.00  0.00           H
ATOM     67  CZ  PHE A   4      18.507  34.929  -1.643  1.00 43.95           C
ATOM     68  HZ  PHE A   4      18.401  35.892  -1.144  1.00  0.00           H
END



here my code :

from Bio import PDB

from Bio.PDB import PICIO, PDBIO

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser(PERMISSIVE=1, QUIET=1)

structure: PDB.Structure.Structure = parser.get_structure("prova", "example_short.pdb")

structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
d:  Dict[Tuple[PDB.internal_coords.AtomKey,
PDB.internal_coords.AtomKey,
PDB.internal_coords.AtomKey,
PDB.internal_coords.AtomKey],
PDB.internal_coords.Dihedron] = ic_chain.dihedra

cnt = 1
for key in d:

if key[0].akl[3] == 'N':
if key[1].akl[3] == 'CA':
if key[2].akl[3] == 'C':
if key[3].akl[3] == 'N':

print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

d[key].angle += 45

cnt += 1

cnt = 1
for key in d:

if key[0].akl[3] == 'N':
if key[1].akl[3] == 'CA':
if key[2].akl[3] == 'C':
if key[3].akl[3] == 'N':

print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

cnt += 1

structure.internal_to_atom_coordinates(verbose = True)

io = PDBIO()

io.set_structure(structure)

io.save('atom_coord.pdb',  preserve_atom_numbering=True)


here an image of the original vs new quadripeptid (cartoon repr, just backbone) :

original is rainbowed one (blue to red), magenta is the new one

saved as atom_coord.pdb

I believe I changed φ (Phi) adding 45° since

Now, a φ angle is N-CA-C-N

It is worth mentioning that changing

 d[key].angle += 45 to d[key].angle += 6645


works too. Apparently, need confirmation from someone more knowledgeable,

structure.internal_to_atom_coordinates


can cope with more than 360° angle values.

• A nice and simple experiment is making a spaghetti (spaghetto?) into a helix with &phi;=-57.8 and &psi;=-47.0 Aug 25, 2022 at 14:24
• For the record, I wanted to award the bounty to this answer, but it was auto-assigned :( Aug 29, 2022 at 20:13

I do not know if there is a function in Biopython.PDB for dihedral space transformations, but due to the lack of answers I would say no. Biopython.PDB is nice and clean, but I would use something with more molecular mechanics power... or should I say potential (:drum::laughing:).

Personally, I use PyRosetta, but it is not open source (code available to academics), so others may opt for OpenForceField. PyRosetta requires a licence (free for academia), which is different from the Rosetta proper licence.

The code below uses a helper module I wrote, so I'll use the installer therein (which is shorthand for the official pip wheel, which provides different options etc., while the official shorthand is actually pyrosettacolabsetup)

pip3 install pyrosetta-help
install_pyrosetta -u 👾👾👾 -p 👾👾👾


## Ligands?

First lets make sure there's nothing interesting ligand wise.

NB. RCSB has a new shiny API which may be better... Not tried it: currently using European.

import requests

def retrieve_nonpolymer_data(entry_id: str) -> dict:
"""cf. https://www.ebi.ac.uk/pdbe/api/doc/"""
response: requests.Response = requests.get(f' https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{entry_id.lower()}')
response.raise_for_status()
return response.json()[entry_id.lower()]

info = retrieve_nonpolymer_data('1enh')
chem_resns = [c for entity in info for c in entity.get('chem_comp_ids', [])]
print(chem_resns)  # ['HOH']


If is were something that is not water, I made a list here, which is based off Peter Curran's repo for the HotSpots algorithm, which allows the elimination of worthless ligands like glycol.

I mention this because the route differs if there is a ligand.

import pyrosetta
import pyrosetta_help as ph
from types import ModuleType

logger = ph.configure_logger()
pyrosetta.init(extra_options=ph.make_option_string(no_optH=False,
ex1=None,
ex2=None,
#mute='all',
ignore_unrecognized_res=True,
ignore_waters=True)
)

# shorthand for sanity
prc: ModuleType = pyrosetta.rosetta.core

# this route strips ligands
pose: Pyrosetta.Pose = pyrosetta.toolbox.pose_from_rcsb('1enh')


Note ignore_unrecognized_res, load_PDB_components and ignore_waters.

## Angles

This will not happen to an non-Rosetta user, but I get tripped out by the fact that ψ & φ data is not stored in residue, but in the pose.

for residue in pose.residues:  #: prc.conformation.Residue
# residue.get_phi() # ---> Error!
pass


For the residue index in pose, Rosetta was written in Fortran so is one-indexed.

import numpy as np
import numpy.typing as npt
phi:npt.NDArray[float] = np.fromiter(map(pose.phi, range(1, pose.total_residue()+1)), dtype=float)
psi:npt.NDArray[float] = np.fromiter(map(pose.psi, range(1, pose.total_residue()+1)), dtype=float)


## Change φ/ψ angles

Time for some fun in dihedral space, say:

r = 2
phi = 180
pose.set_phi(r, phi)


For modelling, I wrote some code for setting spans with the archetypical angles for different SS: here. I do not recall whence the data came (maybe Wikipedia?), but it may be handy anyway.

## Save

To save a PDB file:

pose.dump_pdb('👾👾👾.pdb')


There is also pose.dump_scored_pdb, but that is only for friends of J. Willard Gibbs.

## Minimisation

It is said in the comments that energy minimisation is not required. But here is a snippet that minimises the sidechains only

scorefxn = pyrosetta.get_fa_scorefxn()  # ref2015 will do.
before: float = scorefxn(pose)
cycles = 5
movemap = pyrosetta.MoveMap()
movemap.set_bb(False)  # backbone
movemap.set_chi(True)  # sidechain
relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, cycles)
relax.set_movemap(movemap)
relax.apply(pose)
after: float = scorefxn(pose)
print(after - before)  # Willard says "Remember free energy is a potential so it should be negative"


If the BB tweaks are minimal, there's also the option of scoring/minimising against the electron density etc.

• Thanks for the answer. In BioPython there is functionality for setting angles in the internal_to_atom_coordinates() method. To test your solution, could you edit to a single script that does Steps 1-5 and sets all the phi/psi angles? I am not familiar with PyRosetta. Something that worries me: are you doing some behind-the-scenes optimization? I just want to make the transformation of going from dihedral to Cartesian space. A test of correctness: get back the modified dihedrals from the reconstructed atom coordinates. Aug 24, 2022 at 9:05
• Uhm... It does strike me as odd that you cannot alter the dihedrals in that case —let me check. Bio.PDB is rather basic so I would assume it might be convoluted hence the lack of Bio.PDB replies. PyRosetta does add protons during reading and strips waters if told to do so. If told to autoload the chemical components (ligands) and the structure has an anomaly in its ligand (i.e. it is not what it ought to be) it will not be distorted on load unless minimised. Aug 24, 2022 at 13:11

I answered about doing it with PyRosetta as I assumed that as BioPython PDB module is basic and the fact it was being asked made me think it was not possible. Turns out I was wrong. So I am making a separate answer —as I do not actually believe doing it with Bio.PDB is the best solution.

The provided code from above with type hints is:

import numpy as np
from Bio import PDB

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser()
structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
dihedrals:  Dict[Tuple[PDB.internal_coords.AtomKey,
PDB.internal_coords.AtomKey,
PDB.internal_coords.AtomKey,
PDB.internal_coords.AtomKey],
PDB.internal_coords.Dihedron] = ic_chain.dihedra


I got the types simply by using type and the attributes of note with dir. The latter helps me explore the Bio.PDB.internal_coords.Dihedron telling me of the angle attribute:

d: PDB.internal_coords.Dihedron = list(dihedrals.values())[0]
d.angle: np.float64


A change here appears it the main PDB.Structure.Structure instance, i.e. it is not a disconnected copy.

The class PDB.internal_coords.AtomKey appears to be hard to instantiate by humans and may not work as expected —help is gibberish so it is a case for inspect.getsource. So iterating across the dihedrals for the Ramachandran dihedrals is the solution. Now, a φ angle is N-CA-C-N and a ψ is C-N-CA-C (in PDB atom names the carboxylic carbon is C ). When the __str__ method of one of these PDB.internal_coords.AtomKey instances is called some values are given, which are available in the tuple attribute akl

atom_key: PDB.internal_coords.AtomKey = list(ic_chain.dihedra.keys())[0][0]
print(atom_key)  #3_R_N
print(atom_key.akl) #('3', None, 'R', 'N', None, None)


The first entry is the residue number and is a string because of insertion codes. The fourth entry is the atom name. So iterating across the keys for the correct sequence is the trick —I know there's bounty asking for a complete solution, but it's simple, so anyone care to fill out the blanks and get the meaningless points, please do!

• Thanks. Indeed I am actually looking for a complete solution that takes care also of the 'meaningless points' and that has a test of correctness. Aug 24, 2022 at 16:43
• I am at a loss here, one question so biopython doesnt give you just the backbone/Ramachandran dihedrals but you need to sort them out ? Aug 24, 2022 at 19:21