analyser.densityObj.getPointDensityFromXyz and coordinates from a PDB file

[aliaksei-chareshneu]

First, I would like to thank @hunter-moseley and @syao13 for developing such a nice package for working with electron densities. It is very comfortable to work with.

The question:

Basically, I would like to get the electron density values corresponding to the position of certain atoms from a PDB file. Based on the documentation, parsed the CCP4 file corresponding to the PDB file, and then used this approach:

analyser.densityObj.getPointDensityFromXyz([10.1, 15.2, 24.4])

It seems to produce valid results. However, another similar package I tried (https://github.com/project-gemmi/gemmi) produces quite different results for the same input. I'm not sure which one is correct.

Could you confirm please that using analyser.densityObj.getPointDensityFromXyz in a way I described (providing to it the xyz coordinates from PDB file for some atom) will indeed output the electron density for that atom?

[hunter-moseley]

I do not know about the gemmi project and you did not indicate how you were using their codebase to get a value. However, you must make sure you are analyzing the same density map.
analyzer.densityObj is the 2Fo-Fc map downloaded from PDBe. analyzer.diffDensityObj is the Fo-Fc map downloaded from PDBe.

Which map are you trying to analyze with the gemmi codebase?

Warm regards, Hunter Moseley

[aliaksei-chareshneu]

Thank you very much for the quick response. I have been using the 2Fo-Fc map (analyzer.densityObj). As for the codebase, I attached the both PDBEDA and gemmi solutions just in case.

The actual numbers:

Atom name | PDBEDA | gemmi | difference
C11 | 0.488911033 | 0.525609374 | -0.0367
C21 | 0.569632411 | 0.54194051 | 0.027692
C31 | 0.681316018 | 0.624416947 | 0.056899
C41 | 0.787282884 | 0.749854207 | 0.037429
C51 | 0.649773419 | 0.607901871 | 0.041872
C61 | 0.307238728 | 0.318572879 | -0.01133

I have an idea why the results could differ. Could you tell me please if PDBEDA analyser.densityObj.getPointDensityFromXyz is doing some interpolation, or it just takes the nearest grid point value?

Best wishes, Aliaksei

PDBEDA solution:

import pdb_eda

pdbid = '6fux'
analyser = pdb_eda.densityAnalysis.fromPDBid(pdbid)

dict = {}
for model in analyser.biopdbObj:
    for chain in model:
        for residue in chain:
           for atom in residue:
                dict[atom.get_name()] = atom.get_coord()

atom_names = ["C11", "C21", "C31", "C41", "C51", "C61"]
print("Atom name", "Coordinates", "ED value")
for key in dict:
    if key in atom_names:
        print(key, dict[key], analyser.densityObj.getPointDensityFromXyz(dict[key]))

gemmi solution:

import gemmi
import numpy as np

# Obtained from: https://www.ebi.ac.uk/pdbe/entry/pdb/6fux
str = gemmi.read_pdb('./6fux/pdb6fux.ent')
map = gemmi.read_ccp4_map('./6fux/6fux.ccp4')
map.setup()

dict = {}
for model in str:
    for chain in model:
        for res in chain:
            # ligand residue I am interested in
            # more info: https://www.ebi.ac.uk/pdbe/entry/pdb/6fux/bound/SRY
            if res.name == 'SRY':
                for atom in res:
                    dict[atom.name] = atom.pos            

# names of atoms I am interested in
atom_names = ["C11", "C21", "C31", "C41", "C51", "C61"]

# iterating over atoms and printing interpolated intensity
print("Atom name", "Coordinates", "ED value")
for key in dict:
    if key in atom_names:
        print(key, dict[key], map.grid.interpolate_value(dict[key]))

[hunter-moseley]

Yes, gemmi is calculating a trilinear interpolation value, whereas pdb_eda is returning the closest voxel value without interpolation. https://en.wikipedia.org/wiki/Trilinear_interpolation

This can be seen at line 282 in the gemmi grid.hpp file: https://github.com/project-gemmi/gemmi/blob/master/include/gemmi/grid.hpp

I hope this helps.

[aliaksei-chareshneu]

Thank you so much! It is clear now. I'm closing the issue then.