Error with bfactors

[jdmadura]

It appears that when the load_pdb_rcsb() executes it does not grab the bfactor data from the pdb file. When I run the following python code the bf array is zero when in fact it should have values (see the file 4pti at rcsb.org).

import pytraj as pt pdb = pt.load_pdb_rcsb('4pti') bf = bfactors(pdb, '@CA') print(bf)

When I replace 4pti with 1l2y the bf array contains non-zero values. What is strange is that at rcsb.org the bfactor column for 1l2y is 0.0.

[hainm]

hi,

you are right. pt.load_pdb_rcsb only load coordinates. pytraj is mainly for trajectory analysis from simulation.

pt.bfactors computes bfactor-like from trajectory.

If you want to get the bfactor from xray structure, you should use parmed (https://github.com/ParmEd/ParmEd)

pip install git+https://github.com/ParmEd/ParmEd

import parmed as pmd
pdb = pmd.download_PDB('4pti')
print([atom.bfactor for atom in pdb.atoms])

Let us know if this works for you.

[jdmadura]

Dear Hai,

Yes that worked perfectly… Thank-you for the prompt reply.

Best Regards,

 Jeffry

On Feb 23, 2016, at 10:51 PM, Hai Nguyen notifications@github.com wrote:

hi,

you are right. pt.load_pdb_rcsb only load coordinates. pytraj is manly for trajectory analysis from simulation.

pt.bfactors computes bfactor-like from trajectory.

If you want to get the bfactor from pdb bank, you should use parmed (https://github.com/ParmEd/ParmEd https://github.com/ParmEd/ParmEd)

pip install git+https://github.com/ParmEd/ParmEd import parmed as pmd pdb = pmd.download_PDB('4pti') print([atom.bfactor for atom in pdb.atoms]) Let us know if this works for you.

— Reply to this email directly or view it on GitHub https://github.com/Amber-MD/pytraj/issues/1152#issuecomment-188050776.

Professor ACS Fellow Lambert F. Minucci Endowed Chair in Engineering and Computational Sciences Department of Chemistry & Biochemistry Duquesne University 308 Mellon Hall 600 Forbes Avenue Pittsburgh, PA 15282

Phone: 412-396-4129 FAX: 412-396-5683

e-mail: madura@duq.edu

co-Editor: Journal of Molecular Graphics & Modeling Computational Biology and Chemistry

[drroe]

you are right. pt.load_pdb_rcsb only load coordinates. pytraj is mainly for trajectory analysis from simulation.

That is true, but you should probably preserve as much info from the input PDB as you can (for example cpptraj stores things like B-factors, chain IDs, etc). The reason is to make coordinate conversion as complete as possible.

[hainm]

@drroe

behind the scene: with load_pdb_rcsb, pytraj first download pdb file from rcsb and then call cpptraj to read it. So I think all info that cpptraj read is still there.

Does cpptraj keep info like bfactor in 'parm' command or ib 'trajin' command? Do i need to add any keyword?

[hainm]

@jdmadura: i am glad it helps. Feel free to open new issues if you have any question. Cheers.

[drroe]

Does cpptraj keep info like bfactor in 'parm' command or ib 'trajin' command? Do i need to add any keyword?

This info is stored in the topology. The 'trajin' etc commands read coordinates only. You don't need to add any keyword, it is included if present. See e.g. https://github.com/drroe/cpptraj/blob/217c7f756211cdd97707da273432a827a0978e64/src/Traj_PDBfile.cpp#L425-L429

[hainm]

great. So extra info from pdb file is still there but hidden in Topology.

In [1]: pt.fetch_pdb('4pti').save('test.pdb')

In [2]: !head test.pdb
ATOM      1  N   ARG A   1      26.465  27.452  -2.490  1.00 25.18           N
ATOM      2  CA  ARG A   1      25.497  26.862  -1.573  1.00 17.63           C
ATOM      3  C   ARG A   1      26.193  26.179  -0.437  1.00 17.26           C

question: should cpptraj add those info to atom too?

pdb = pt.load_topology(pdbfile)
print([atom.bfactor for atom in pdb.atoms])

[drroe]

question: should cpptraj add those info to atom too?

The reason I put it in a separate array in Topology (extra_) is that this information is really only specific to PDB/CIF files. I suppose it could be added to Atom and just set to zero, but that's not a simple change and not worthwhile this close to code deadline IMO. Something we can maybe chat about at the dev meeting?

[hainm]

I suppose it could be added to Atom and just set to zero, but that's not a simple change and not worthwhile this close to code deadline IMO. Something we can maybe chat about at the dev meeting?

so we just direct users to use ParmEd. It is handy for those cases.

[hainm]

I go ahead closing this issue since it is solved.