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']]
