MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.31k stars 648 forks source link

Don't add segID to PDB files unless requested #3505

Open sbhakat opened 2 years ago

sbhakat commented 2 years ago

I want to add only chain ID to my PDB but couldn't manage to add chain ID without adding segment ID. Is there anyway to add only chain? My script looks like following

import MDAnalysis as mda
U = mda.Universe("prot-mass-nochain.pdb")
newseg = U.add_Segment(segid='A')
U.select_atoms('bynum 1:2410').residues.segments = newseg
U.atoms.write("prot-mass-aadedchain.pdb")

The output looks like

ATOM   2407  CD2 HIS A 299      -1.606  84.378 319.924  1.00  0.00      A    C
ATOM   2408  C   HIS A 299       1.108  85.007 322.375  1.00  0.00      A    C
ATOM   2409  OC1 HIS A 299       1.253  85.851 321.410  1.00  0.00      A    O
ATOM   2410  OC2 HIS A 299       1.270  85.194 323.562  1.00  0.00      A    O

The reason being that CHARMM GUI doesn't accept segment ID as a part of RCSB PDB format. Any help will be highly appreciated.

IAlibay commented 2 years ago

@sbhakat can you confirm which version of MDAnalysis you are using? This part of the PDBWriter has changed a lot in recent iterations, so it's worth asking.

sbhakat commented 2 years ago

@IAlibay using version 1.1.1

IAlibay commented 2 years ago

Thanks @sbhakat

So technically it's possible to get around this, but we should also have the option to not add in segments to PDB files - if it's not in the standard, it shouldn't be done by default.

Getting around the issue:

In 2.0.0 the PDBWriter was changed to used chainID instead of segID for setting chains.

You would have to upgrade to 2.0.0, but in your above case, if you did:

# if you have any segments then set them to an empty string to avoid writing them
u.atoms.segements.segids = ''
# now add a chainIDs attribute for your chains (if you don't have one already)
u.add_TopologyAttr('chainIDs')
# set chain for atoms 0 through to 2409
u.atoms[:2410].chainIDs = 'A'

instead of:

newseg = U.add_Segment(segid='A')

It should fix the issue? (at least it does in the test case I'm looking at).

sbhakat commented 2 years ago

@IAlibay thanks a lot. Worked :D

IAlibay commented 2 years ago

I think I'm going to re-open this and change the issue title to "don't write segment by default".

As far as I'm aware it's not in the PDBv3 format: http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM So we shouldn't be enforcing a non-existing standard. Worth having a discussion imho.

sbhakat commented 2 years ago

Well even if I add the chain ID correctly I still can't upload the PDB in CHARMM GUI. Seems like some formatting problem. Any inputs. Attaching the PDB for convenience. https://www.dropbox.com/s/hn30mr4w6lr1sfk/prot-mass-newchain.pdb?dl=0 #

sbhakat commented 2 years ago

Well. I solved this issue by using mdconvert (mdconvert test.pdb -o test-format.pdb) on top of the output of mdanalysis . Output here: https://www.dropbox.com/s/bpeno8t8oan4auk/test.pdb?dl=0 CHARMM GUI can capture this