Answers.Body: <p>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:).</p>
  <p>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.</p>
  <p>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 <code>pyrosettacolabsetup</code>)</p>
  <pre class="lang-bash prettyprint-override"><code>pip3 install pyrosetta-help
  install_pyrosetta -u  -p 
  </code></pre>
  <h2>Ligands?</h2>
  <p>First lets make sure there's nothing interesting ligand wise.</p>
  <p>NB. RCSB has a new shiny API which may be better... Not tried it: currently using European.</p>
  <pre class="lang-py prettyprint-override"><code>import requests
  
  def retrieve_nonpolymer_data(entry_id: str) -&gt; dict:
      &quot;&quot;&quot;cf. https://www.ebi.ac.uk/pdbe/api/doc/&quot;&quot;&quot;
      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']
  </code></pre>
  <p>If is were something that is not water, I made a <a href="https://blog.matteoferla.com/2019/11/go-away-glycerol.html" rel="nofollow noreferrer">list here</a>, which is based off Peter Curran's repo for the HotSpots algorithm, which allows the elimination of worthless ligands like glycol.</p>
  <p>I mention this because the route differs if there is a ligand.</p>
  <h2>Init and read</h2>
  <pre class="lang-py prettyprint-override"><code>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,
                                                  load_PDB_components=False,
                                                  ignore_waters=True)
                                 )
  
  # shorthand for sanity
  prc: ModuleType = pyrosetta.rosetta.core
  
  # this route strips ligands
  pose: Pyrosetta.Pose = pyrosetta.toolbox.pose_from_rcsb('1enh')
  </code></pre>
  <p>Note <code>ignore_unrecognized_res</code>, <code>load_PDB_components</code> and <code>ignore_waters</code>.</p>
  <h2>Angles</h2>
  <p>This will not happen to an non-Rosetta user, but I get tripped out by the fact that ψ &amp; φ data is not stored in residue, but in the pose.</p>
  <pre class="lang-py prettyprint-override"><code>for residue in pose.residues:  #: prc.conformation.Residue
      # residue.get_phi() # ---&gt; Error!
      pass
  </code></pre>
  <p>For the residue index in pose, Rosetta was written in Fortran so is one-indexed.</p>
  <pre class="lang-py prettyprint-override"><code>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)
  </code></pre>
  <h2>Change φ/ψ angles</h2>
  <p>Time for some fun in dihedral space, say:</p>
  <pre class="lang-py prettyprint-override"><code>r = 2
  phi = 180
  pose.set_phi(r, phi)
  </code></pre>
  <p>For modelling, I wrote some code for setting spans with the archetypical angles for different SS: <a href="https://pyrosetta-help.readthedocs.io/en/latest/pyrosetta_help.common_ops.html?highlight=helix#pyrosetta_help.common_ops.ss_changes.make_sheet" rel="nofollow noreferrer">here</a>. I do not recall whence the data came (maybe Wikipedia?), but it may be handy anyway.</p>
  <h2>Save</h2>
  <p>To save a PDB file:</p>
  <pre class="lang-py prettyprint-override"><code>pose.dump_pdb('.pdb')
  </code></pre>
  <p>There is also <code>pose.dump_scored_pdb</code>, but that is only for friends of J. Willard Gibbs.</p>
  <h2>Minimisation</h2>
  <p>It is said in the comments that energy minimisation is not required. But here is a snippet that minimises the sidechains only</p>
  <pre class="lang-py prettyprint-override"><code>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 &quot;Remember free energy is a potential so it should be negative&quot;
  </code></pre>
  <p>If the BB tweaks are minimal, there's also the option of scoring/minimising against the electron density etc.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-24T09:05:02.990Z
Answers.Text: Thanks for the answer. In BioPython there is functionality for setting angles in the [`internal_to_atom_coordinates()`](https://biopython.org/docs/1.77/api/Bio.PDB.internal_coords.html) 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.
Answers.UserDisplayName: epsilone
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-24T13:11:26.057Z
Answers.Text: 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.
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-23T13:28:04.517Z
Answers.LastActivityDate: 2022-08-23T19:10:37.830Z
Answers.LastEditDate: 2022-08-23T19:10:37.830Z
Answers.UserDisplayName: Matteo Ferla
Answers.Body: <p>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.</p>
  <p>The provided code from above with type hints is:</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  from Bio import PDB
  
  from typing import TypedDict, Dict, Tuple
  
  parser = PDB.PDBParser()
  structure: PDB.Structure.Structure = parser.get_structure(&quot;&quot;, &quot;/Users/matteo/Downloads/1enh.pdb&quot;)
  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
  </code></pre>
  <p>I got the types simply by using <code>type</code> and the attributes of note with <code>dir</code>. The latter helps me explore the <code>Bio.PDB.internal_coords.Dihedron</code> telling me of the <code>angle</code> attribute:</p>
  <pre class="lang-py prettyprint-override"><code>d: PDB.internal_coords.Dihedron = list(dihedrals.values())[0]
  d.angle: np.float64
  </code></pre>
  <p>A change here appears it the main <code>PDB.Structure.Structure</code> instance, i.e. it is not a disconnected copy.</p>
  <p>The class <code>PDB.internal_coords.AtomKey</code> appears to be hard to instantiate by humans and may not work as expected —<code>help</code> is gibberish so it is a case for <code>inspect.getsource</code>. 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 <code>C </code>). When the <code>__str__</code> method of one of these <code>PDB.internal_coords.AtomKey</code> instances is called some values are given, which are available in the tuple attribute <code>akl</code></p>
  <pre class="lang-py prettyprint-override"><code>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)
  </code></pre>
  <p>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!</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-24T16:43:58.853Z
Answers.Text: 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.
Answers.UserDisplayName: epsilone
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-24T19:21:46.820Z
Answers.Text: 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 ?
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-24T13:46:36.817Z
Answers.LastActivityDate: 2022-08-24T13:46:36.817Z
Answers.UserDisplayName: Matteo Ferla
Answers.Body: <p>ok, here my input pdb, <code>example_short.pdb</code>  :</p>
  <pre><code>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
  
  </code></pre>
  <p>here my code :</p>
  <pre><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(&quot;prova&quot;, &quot;example_short.pdb&quot;)
  
  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) 
  </code></pre>
  <p>here an image of the original vs new quadripeptid (cartoon repr, just backbone) :</p>
  <p>original is rainbowed one (blue to red), magenta is the new one</p>
  <p>saved as <code>atom_coord.pdb</code></p>
  <p><a href="https://i.sstatic.net/3OtQt.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/3OtQt.png" alt="enter image description here" /></a></p>
  <p>I believe I changed φ (Phi) adding 45° since</p>
  <blockquote>
  <p>Now, a φ angle is N-CA-C-N</p>
  </blockquote>
  <p>It is worth mentioning that changing</p>
  <pre><code> `d[key].angle += 45` to `d[key].angle += 6645`
  </code></pre>
  <p>works too. Apparently, need confirmation from someone more knowledgeable,</p>
  <pre><code>`structure.internal_to_atom_coordinates`
  </code></pre>
  <p>can cope with more than 360° angle values.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-25T14:24:14.390Z
Answers.Text: A nice and simple experiment is making a spaghetti (spaghetto?) into a helix with &phi;=-57.8 and &psi;=-47.0
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-29T20:13:07.747Z
Answers.Text: For the record, I wanted to award the bounty to this answer, but it was auto-assigned :(
Answers.UserDisplayName: epsilone
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-25T09:46:34.193Z
Answers.LastActivityDate: 2022-08-25T13:06:26.193Z
Answers.LastEditDate: 2022-08-25T13:06:26.193Z
Answers.UserDisplayName: pippo1980
Body: <p>As the title summarizes, I am trying to:</p>
  <ol>
  <li>Read a PDB file (for example, <a href="https://files.rcsb.org/download/1ENH.pdb" rel="nofollow noreferrer"><code>1enh.pdb</code></a>).</li>
  <li>Extract the backbone dihedral angles (phi, psi).</li>
  <li>Modify the dihedral angles (phi, psi) arbitrarily (for example, replace them with <code>np.random.uniform(-np.pi, np.pi, size = (56, 2))</code> in 1enh, since there are 56 pairs of dihedrals).</li>
  <li>Reconstruct the Cartesian coordinates of the modified backbone.</li>
  <li>Write a PDB file based on the modified backbone.</li>
  </ol>
  <p>From <a href="https://bioinformatics.stackexchange.com/questions/18849/help-using-biopythons-ic-chain-functionality-for-getting-protein-3d-coordinates">this question</a> it seems the tools for what I want are available in BioPython. There is also a <code>set_angle()</code> 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.</p>
  <p>This is what I have tried so far:</p>
  <pre><code>#!/usr/bin/python
  
  import numpy as np
  from Bio.PDB import PDBParser, PICIO, PDBIO
  
  ## Step 1
  
  parser = PDBParser()
  structure = parser.get_structure(&quot;&quot;, &quot;1enh.pdb&quot;)
  
  ## 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 = PICIO.read_PIC('internal')
  modified.internal_to_atom_coordinates()
  
  ## Step 5
  
  io = PDBIO()
  io.set_structure(modified)
  io.save('1enh_modified.pdb', preserve_atom_numbering = True) 
  </code></pre>
  <p>Many thanks in advance.</p>
  <p><strong>EDIT:</strong> 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.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-20T09:59:16.163Z
Comments.Text: 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 ?
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-20T10:57:25.340Z
Comments.Text: 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.
Comments.UserDisplayName: epsilone
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-20T11:12:07.563Z
Comments.Text: Any Python package that checks PDB geometry  just for clashes ??
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-20T11:50:18.360Z
Comments.Text: 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.
Comments.UserDisplayName: epsilone
Tags: <python><protein-structure><pdb><3d-structure><biopython>
Title: Read PDB file, extract dihedral angles, modify dihedral angles, reconstruct Cartesian coordinates, and write PDB file
ThreadText: # Read PDB file, extract dihedral angles, modify dihedral angles, reconstruct Cartesian coordinates, and write PDB file
  
  (Asked by: username_0 on 2022-08-19)
  
  As the title summarizes, I am trying to:
  
  1. Read a PDB file (for example, [`1enh.pdb`](https://files.rcsb.org/download/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](https://bioinformatics.stackexchange.com/questions/18849/help-using-biopythons-ic-chain-functionality-for-getting-protein-3d-coordinates) 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 = PICIO.read_PIC('internal')
  modified.internal_to_atom_coordinates()
  
  ## Step 5
  
  io = PDBIO()
  io.set_structure(modified)
  io.save('1enh_modified.pdb', preserve_atom_numbering = True) 
  
  ```
  
  Many thanks in advance.
  
  **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.
  
  **Tags:** python, protein-structure, pdb, 3d-structure, biopython
  
  > **Comments:**
  > **username_1**: 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 ?  
  > **username_0**: 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.  
  > **username_1**: Any Python package that checks PDB geometry  just for clashes ??  
  > **username_0**: 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.  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_1 on 2022-08-25. Score: 5 | Accepted answer
  
  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`
  
  [![enter image description here](https://i.sstatic.net/3OtQt.png)](https://i.sstatic.net/3OtQt.png)
  
  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.
  
  > **Comments:**
  > **username_2**: A nice and simple experiment is making a spaghetti (spaghetto?) into a helix with &phi;=-57.8 and &psi;=-47.0  
  > **username_0**: For the record, I wanted to award the bounty to this answer, but it was auto-assigned :(  
  
  ---
  
  ### Answer by username_2 on 2022-08-23. Score: 5
  
  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`)
  
  ```bash
  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.
  
  ```py
  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](https://blog.username_2.com/2019/11/go-away-glycerol.html), 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.
  
  ## Init and read
  
  ```py
  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,
                                                  load_PDB_components=False,
                                                  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.
  
  ```py
  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.
  
  ```py
  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:
  
  ```py
  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](https://pyrosetta-help.readthedocs.io/en/latest/pyrosetta_help.common_ops.html?highlight=helix#pyrosetta_help.common_ops.ss_changes.make_sheet). I do not recall whence the data came (maybe Wikipedia?), but it may be handy anyway.
  
  ## Save
  
  To save a PDB file:
  
  ```py
  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
  
  ```py
  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.
  
  > **Comments:**
  > **username_0**: Thanks for the answer. In BioPython there is functionality for setting angles in the [`internal_to_atom_coordinates()`](https://biopython.org/docs/1.77/api/Bio.PDB.internal_coords.html) 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.  
  > **username_2**: 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.  
  
  ---
  
  ### Answer by username_2 on 2022-08-24. Score: 4
  
  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:
  
  ```py
  import numpy as np
  from Bio import PDB
  
  from typing import TypedDict, Dict, Tuple
  
  parser = PDB.PDBParser()
  structure: PDB.Structure.Structure = parser.get_structure("", "/Users/username_2/Downloads/1enh.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
  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:
  
  ```py
  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`
  
  ```py
  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!
  
  > **Comments:**
  > **username_0**: 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.  
  > **username_1**: 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 ?  
  
  ---
  
  