Answers.Body: <h2>Per Atom</h2>
  <p>You are correct. Rosetta scorefunction does not store any per atom data. The scoring operates at the per residue level. Whereas each atom has its coordinates and properties in full atom mode, in Pyrosetta it is clear that an atom is a just part of a residue and every operation is applied at the residue level. It's a team effort: the functional group of the side chain works as one and in most cases two atoms have the same role (e.g. <code>OE1</code> and <code>OE2</code> in glutamate behave identically).</p>
  <p>So, the real question is whether you really care about the atom level in a classical mechanical setting?</p>
  <h2>Formal charges</h2>
  <p>A wee parenthesis, we need to make a difference between formal charge and Marsilli-Gasteiger partial charge —one is an integer and the other is a float. Both require a stated protonation state in most applications —say <code>aspartate</code> vs. <code>aspartic acid</code>.</p>
  <p><code>PDB2PQR</code> formats to a PQR file for the ligand docking application <code>autodock</code>. It does not assign how a residue is protonated.
  In the case of Rosetta, the partial charge depends on the 4 column of the topology (params file) of the amino acid, i.e. it is fixed for all amino acids —Rosetta does not implement a Drude particle to make something polarisable. I very strongly suspect this is the same in PDB2PQR.</p>
  <p>Therefore, something like the Advance Poisson-Boltzmann Solver may give more interesting data in terms of partial charges, but also requires the uses to specify the protonation state.</p>
  <p>A small parenthesis is histidine, which can adopt four forms. Two tautomers at neutral pH: single proton on either the epsilon (<code>HIE</code> in Amber) or delta nitrogen (<code>HID</code>), while it can be deprotonated (anion) or diprotonated (<code>HID</code>, cation). Rosetta switches between these two automatically.</p>
  <h2>Properties</h2>
  <p>You list a few properties of interest, which all seem to be based on the residue a priori of structure.</p>
  <p>For example:</p>
  <ul>
  <li>hydrophobicity, unless you have a polarisable model (Drude particle) this is structure independent. RSA, relative surface area, is a must better indication (per residue), which then can get classified as core or surface based on this. Per atom would require just a wee tinker.</li>
  <li>aromaticity, this is also structure independent, but π-π stacking or π-cation interactions are noted at the residue level.</li>
  <li>hydrogen bond acceptor/donor. Whether an atom is a donor or acceptor is actually defined in the atom type used in the topology file. Whether that atom fulfils that bonding is a different matter. If you really really wanted to get per atom values of this, in Rosetta you could exclude the backbone atoms you could mutate to alanine without minimising and see if it's maintained. So differentiating between backbone and sidechain is fine, but per atom is a different matter.</li>
  </ul>
  <p>A structure-dependent property is how flexible the residue is, as measured as an RMS. This can be calculated per atom. This is akin to B-factors, which is how big is the &quot;blur&quot; of the atom. However, if you get the <code>NZ</code> atom on a lysine that is switching between two or more positions you get a partial occupancy. If this atom is basically hand-waving along an axis in an arc you get an anisotropic B-factor. In other words it does not vibrate in all 3 cartesian axes equally. So a B-factor is a terrible approximation. An RSMD of an ensemble is even worse.
  In Rosetta, you can a quick vibration using Backrub, which is very quick at telling you how much motion does an atom experience.</p>
  <h2>Here be dragons</h2>
  <blockquote>
  <p>Warning: As mentioned, it takes two atoms to resulting in a force, hence why per atom scores are a bad idea.</p>
  </blockquote>
  <p>However if for some <em>exceptional</em> reason this is required, say
  figuring out which group on a ligand is causing a clash, in pyrosetta it is possible:</p>
  <pre><code>import pandas as pd
  import pyrosetta
  pyrosetta.init()
  pose = pyrosetta.pose_from_sequence(&quot;ELVISISALIVE&quot;)
  target_resi = 9
  residue = pose.residue(target_resi)
  scores = []
  for i in range(1, residue.natoms() + 1):
      iname = residue.atom_name(i)
      for r in range(1, pose.total_residue() + 1):
          other = pose.residue(r)
          for o in range(1, other.natoms() + 1):
              oname = other.atom_name(o)
              score = pyrosetta.toolbox.atom_pair_energy.etable_atom_pair_energies(residue, 
                                                                                   i,
                                                                                   other,
                                                                                   o, 
                                                                                   scorefxn)
          scores.append(dict(zip(['residue_index_A', 'atom_name_A', 'residue_index_A', 'atom_name_B', 'lj_atr', 'lj_rep', 'fa_solv', 'fa_elec'],
                           [target_resi, residue.atom_name(i), r, other.atom_name(o), *score])))
  
  df = pd.DataFrame(scores)
  df.pivot_table(index='atom_name_A',
                 aggfunc='sum')[['lj_atr','lj_rep','fa_elec','fa_solv']]
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-08T21:44:35.443Z
Answers.LastActivityDate: 2021-02-09T11:11:43.587Z
Answers.LastEditDate: 2021-02-09T11:11:43.587Z
Answers.UserDisplayName: Matteo Ferla
Body: <p>Currently, I am implementing a descriptor of protein structure and would like to calculate properties per atoms like charges, hydrophobicity, hydrogen bond donors / acceptors, hydrophobicity, aromaticity etc.</p>
  <p>I tried using PDB2PQR and succeeded with protein structures but it crashed while calculating nucleic acid bound proteins and choosing other force fields did not help.
  Also, as far as I know, Rosetta tool can calculate hydrophobicity however, it assigns the property per residue but not per atom (correct me if I am wrong).</p>
  <p>Any other suitable tool suggestions would be appreciated.</p>
  
Tags: <pdb><mutations><rosetta>
Title: How to calculate various properties of a protein structure per atom using PDB2PQR and Rosetta tools (or any other tools)?
ThreadText: # How to calculate various properties of a protein structure per atom using PDB2PQR and Rosetta tools (or any other tools)?
  
  (Asked by: username_0 on 2021-01-06)
  
  Currently, I am implementing a descriptor of protein structure and would like to calculate properties per atoms like charges, hydrophobicity, hydrogen bond donors / acceptors, hydrophobicity, aromaticity etc.
  
  I tried using PDB2PQR and succeeded with protein structures but it crashed while calculating nucleic acid bound proteins and choosing other force fields did not help.
  Also, as far as I know, Rosetta tool can calculate hydrophobicity however, it assigns the property per residue but not per atom (correct me if I am wrong).
  
  Any other suitable tool suggestions would be appreciated.
  
  **Tags:** pdb, mutations, rosetta
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-01-08. Score: 2 | Accepted answer
  
  ## Per Atom
  
  You are correct. Rosetta scorefunction does not store any per atom data. The scoring operates at the per residue level. Whereas each atom has its coordinates and properties in full atom mode, in Pyrosetta it is clear that an atom is a just part of a residue and every operation is applied at the residue level. It's a team effort: the functional group of the side chain works as one and in most cases two atoms have the same role (e.g. `OE1` and `OE2` in glutamate behave identically).
  
  So, the real question is whether you really care about the atom level in a classical mechanical setting?
  
  ## Formal charges
  
  A wee parenthesis, we need to make a difference between formal charge and Marsilli-Gasteiger partial charge —one is an integer and the other is a float. Both require a stated protonation state in most applications —say `aspartate` vs. `aspartic acid`.
  
  `PDB2PQR` formats to a PQR file for the ligand docking application `autodock`. It does not assign how a residue is protonated.
  In the case of Rosetta, the partial charge depends on the 4 column of the topology (params file) of the amino acid, i.e. it is fixed for all amino acids —Rosetta does not implement a Drude particle to make something polarisable. I very strongly suspect this is the same in PDB2PQR.
  
  Therefore, something like the Advance Poisson-Boltzmann Solver may give more interesting data in terms of partial charges, but also requires the uses to specify the protonation state.
  
  A small parenthesis is histidine, which can adopt four forms. Two tautomers at neutral pH: single proton on either the epsilon (`HIE` in Amber) or delta nitrogen (`HID`), while it can be deprotonated (anion) or diprotonated (`HID`, cation). Rosetta switches between these two automatically.
  
  ## Properties
  
  You list a few properties of interest, which all seem to be based on the residue a priori of structure.
  
  For example:
  
  * hydrophobicity, unless you have a polarisable model (Drude particle) this is structure independent. RSA, relative surface area, is a must better indication (per residue), which then can get classified as core or surface based on this. Per atom would require just a wee tinker.
  * aromaticity, this is also structure independent, but π-π stacking or π-cation interactions are noted at the residue level.
  * hydrogen bond acceptor/donor. Whether an atom is a donor or acceptor is actually defined in the atom type used in the topology file. Whether that atom fulfils that bonding is a different matter. If you really really wanted to get per atom values of this, in Rosetta you could exclude the backbone atoms you could mutate to alanine without minimising and see if it's maintained. So differentiating between backbone and sidechain is fine, but per atom is a different matter.
  
  A structure-dependent property is how flexible the residue is, as measured as an RMS. This can be calculated per atom. This is akin to B-factors, which is how big is the "blur" of the atom. However, if you get the `NZ` atom on a lysine that is switching between two or more positions you get a partial occupancy. If this atom is basically hand-waving along an axis in an arc you get an anisotropic B-factor. In other words it does not vibrate in all 3 cartesian axes equally. So a B-factor is a terrible approximation. An RSMD of an ensemble is even worse.
  In Rosetta, you can a quick vibration using Backrub, which is very quick at telling you how much motion does an atom experience.
  
  ## Here be dragons
  
  > Warning: As mentioned, it takes two atoms to resulting in a force, hence why per atom scores are a bad idea.
  
  However if for some *exceptional* reason this is required, say
  figuring out which group on a ligand is causing a clash, in pyrosetta it is possible:
  
  ```
  import pandas as pd
  import pyrosetta
  pyrosetta.init()
  pose = pyrosetta.pose_from_sequence("ELVISISALIVE")
  target_resi = 9
  residue = pose.residue(target_resi)
  scores = []
  for i in range(1, residue.natoms() + 1):
      iname = residue.atom_name(i)
      for r in range(1, pose.total_residue() + 1):
          other = pose.residue(r)
          for o in range(1, other.natoms() + 1):
              oname = other.atom_name(o)
              score = pyrosetta.toolbox.atom_pair_energy.etable_atom_pair_energies(residue, 
                                                                                   i,
                                                                                   other,
                                                                                   o, 
                                                                                   scorefxn)
          scores.append(dict(zip(['residue_index_A', 'atom_name_A', 'residue_index_A', 'atom_name_B', 'lj_atr', 'lj_rep', 'fa_solv', 'fa_elec'],
                           [target_resi, residue.atom_name(i), r, other.atom_name(o), *score])))
  
  df = pd.DataFrame(scores)
  df.pivot_table(index='atom_name_A',
                 aggfunc='sum')[['lj_atr','lj_rep','fa_elec','fa_solv']]
  
  ```
  
  
  ---
  
  