ksy141 / mstool

Multiscale Simulation Tool for Backmapping
https://mstool.readthedocs.io/en/latest/#
GNU General Public License v3.0
16 stars 1 forks source link

Making topology files from the Atomistic pdb file #1

Closed maryammajd closed 9 months ago

maryammajd commented 10 months ago

Hi,

I have an issue getting the topology file from the provided Atomistic PDB file using gromacs pdb2gmx. I have a membrane of POPC and POPE, it seems that the pdb file that I get from mstool does not have the same order of atoms for POPE and therefore the resulted top file has mismatched order of atoms. I have attached an image of the top file and the pdb file. It seems that the hydrogen atoms that are attached to the N have different orders in charmm36 forcefield and the provided mapping files in mstools. I was wondering if there is any way to solve this issue.

Bests, Maryam Majd

Screenshot from 2023-12-18 17-53-11 Screenshot from 2023-12-18 17-52-57

ksy141 commented 10 months ago

Hi Maryam,

Thank you for reporting this issue. I think I used the order of atoms to be consistent with the CHARMM36 native force field. It seems gromacs uses a slightly different order.

The order of backmapped atoms in the final system is determined by the order of atoms written in a mapping scheme. In the default file ($mstool/mapping/martini.lipid.c36.dat),

RESI POPE [ NH3 ] N C12 H12A H12B HN1 HN2 HN3

The atom order in the backmapped structure will be N -> C12 -> H12A -> H12B -> HN1 -> HN2 -> HN3. However, gromacs asks the atom order to be N -> HN1 -> HN2 -> HN3 -> C12 -> H12A -> H12B, based on your snapshot.

I made a new mapping file for you that simply has a change in the atom order. Please find it attached. Pass this file to every line. Please let me know if it solves the problem. Otherwise, please send me your structure if you can.

import mstool mstool.Ungroup('cg.pdb', 'aa.pdb', mapping='mapping.dat') mstool.REM('aa.pdb', 'aa_final.pdb', mapping='mapping.dat') mstool.CheckStructure('aa_final.pdb', mapping='mapping.dat')

Best, Siyoung

On Mon, Dec 18, 2023 at 11:52 AM maryammajd @.***> wrote:

Hi,

I have an issue getting the topology file from the provided Atomistic PDB file using gromacs pdb2gmx. I have a membrane of POPC and POPE, it seems that the pdb file that I get from mstool does not have the same order of atoms for POPE and therefore the resulted top file has mismatched order of atoms. I have attached an image of the top file and the pdb file. It seems that the hydrogen atoms that are attached to the N have different orders in charmm36 forcefield and the provided mapping files in mstools. I was wondering if there is any way to solve this issue.

Bests, Maryam Majd

Screenshot.from.2023-12-18.17-51-24.png (view on web) https://github.com/ksy141/mstool/assets/116288311/58562bea-6694-4548-937c-7d93ce2d77ae Screenshot.from.2023-12-18.17-51-10.png (view on web) https://github.com/ksy141/mstool/assets/116288311/144f78e8-3b7d-466d-bb83-92dcb759abfd

— Reply to this email directly, view it on GitHub https://github.com/ksy141/mstool/issues/1, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALSZEC4BQP7K6KV3T2RKKRLYKBYDBAVCNFSM6AAAAABAZ2PLDKVHI2DSMVQWIX3LMV43ASLTON2WKOZSGA2DOMBWGQ3DQNA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

RESI POPC

[ NC3 ] N C12 H12A H12B C13 H13A H13B H13C C14 H14A H14B H14C C15 H15A H15B H15C

[ PO4 ] C11 H11A H11B P O13 O14 O12 O11 C1 HA HB

[ GL1 ] C2 HS O21 C21 O22 C22 H2R H2S

[ GL2 ] C3 HX HY O31 C31 O32 C32 H2X H2Y

[ C1A ] C23 H3R H3S C24 H4R H4S C25 H5R H5S C26 H6R H6S

[ D2A ] C27 H7R H7S C28 H8R H8S C29 H91 C210 H101 C211 H11R H11S

[ C3A ] C212 H12R H12S C213 H13R H13S C214 H14R H14S

[ C4A ] C215 H15R H15S C216 H16R H16S C217 H17R H17S C218 H18R H18S H18T

[ C1B ] C33 H3X H3Y C34 H4X H4Y C35 H5X H5Y C36 H6X H6Y

[ C2B ] C37 H7X H7Y C38 H8X H8Y C39 H9X H9Y C310 H10X H10Y C311 H11X H11Y

[ C3B ] C312 H12X H12Y C313 H13X H13Y C314 H14X H14Y

[ C4B ] C315 H15X H15Y C316 H16X H16Y H16Z

[ cis ] C28 C29 C210 C211

[ chiral ] HS C2 O21 C1 C3 O13 P O11 O14 O12

RESI POPE [ NH3 ] N HN1 HN2 HN3 C12 H12A H12B

[ PO4 ] C11 H11A H11B P O13 O14 O12 O11 C1 HA HB

[ GL1 ] C2 HS O21 C21 O22 C22 H2R H2S

[ GL2 ] C3 HX HY O31 C31 O32 C32 H2X H2Y

[ C1A ] C23 H3R H3S C24 H4R H4S C25 H5R H5S C26 H6R H6S

[ D2A ] C27 H7R H7S C28 H8R H8S C29 H91 C210 H101 C211 H11R H11S

[ C3A ] C212 H12R H12S C213 H13R H13S C214 H14R H14S

[ C4A ] C215 H15R H15S C216 H16R H16S C217 H17R H17S C218 H18R H18S H18T

[ C1B ] C33 H3X H3Y C34 H4X H4Y C35 H5X H5Y C36 H6X H6Y

[ C2B ] C37 H7X H7Y C38 H8X H8Y C39 H9X H9Y C310 H10X H10Y C311 H11X H11Y

[ C3B ] C312 H12X H12Y C313 H13X H13Y C314 H14X H14Y

[ C4B ] C315 H15X H15Y C316 H16X H16Y H16Z

[ chiral ] HS C2 O21 C1 C3 O13 P O11 O14 O12

[ cis ] C28 C29 C210 C211

RESI NA [ NA ] NA

RESI CL [ CL ] CL

RESI ION [ NA ] SOD

[ CL ] CLA

RESI SOD [ SOD ] SOD

RESI CLA [ CLA ] CLA

maryammajd commented 10 months ago

Hi,

Thanks a lot for your help. Yes the new mapping works perfectly. I just needed to change the order of C11 and C12 to be consistant with charmm36 forcefield in gromacs.

I have another question, the bilayer I have consists of many lipids some of which do not have a mapping file. I tried to build mapping file for one lipid, and it seems to work. However, I do not know where I can get the data about cis and chiral of the atoms for these lipids and this makes lipids have some weird configurations. I was wondering if there is any source or tutorial about getting these parameters for different lipids. Thanks in advance for your help and consideration.

Bests, Maryam

ksy141 commented 10 months ago

Hi Maryam,

Glad to hear that. Thank you for updating me.

You can find lipid names and chemical structures on the CHARMM-GUI website. https://charmm-gui.org/?doc=archive&lib=lipid

The chirality specified for phospholipids defined in the default mapping scheme (R-enantiomer at the C2 atom) is pretty common throughout phospholipids. I think you can copy and paste this part for chirality if your lipid is a standard or populated lipid in human cells.

[ chiral ] HS C2 O21 C1 C3 O13 P O11 O14 O12 ; This chirality is completely optional. You can remove or keep it.

Please let me know if you have any questions.

Best, Siyoung

On Tue, Dec 19, 2023 at 5:41 AM maryammajd @.***> wrote:

Hi,

Thanks a lot for your help. Yes the new mapping works perfectly. I just needed to change the order of C11 and C12 to be consistant with charmm36 forcefield in gromacs.

I have another question, the bilayer I have consists of many lipids some of which do not have a mapping file. I tried to build mapping file for one lipid, and it seems to work. However, I do not know where I can get the data about cis and chiral of the atoms for these lipids and this makes lipids have some weird configurations. I was wondering if there is any source or tutorial about getting these parameters for different lipids. Thanks in advance for your help and consideration.

Bests, Maryam

— Reply to this email directly, view it on GitHub https://github.com/ksy141/mstool/issues/1#issuecomment-1862521943, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALSZEC2TQVTRKWVZGQAK2ZTYKFVN7AVCNFSM6AAAAABAZ2PLDKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRSGUZDCOJUGM . You are receiving this because you were assigned.Message ID: @.***>

ksy141 commented 10 months ago

I am closing the issue. Please reopen it if necessary.

maryammajd commented 10 months ago

Hi again Siyoung,

Thank you for your help and instructions. I was trying to make mapping files for the lipids that were not by default in the map files. They worked quite fine with SM headgroups but when I tried to map some cerebroside lipids, I got the following error. I was wondering if there is any way I can solve this issue. I can also provide the mapping file and the pdb file if required. However, github does not accept those file formats, and therefore I could not attach them.

Error: No template found for residue 1 (GLPA). The set of atoms is similar to DXPC, but it is missing 2 atoms. For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#template

Bests, Maryam


ksy141 commented 10 months ago

Hi Maryam,

To use mstool backmapping, there are some requirements:

  1. A residue name defined in an input PDB should be consistent with a residue name in your mapping files and a residue name in the openMM force field. I could not find a residue name, GLPA, in the openMM force field, located at $mstool_dir/FF/charmm36/charmm36.xml.

I could not find GLPA on the CHARMM-GUI website either. I wonder whether it is a new lipid that is not supported in the conventional CHARMM36 force field? If it is a new lipid that the CHARMM36 force field does not support, you should create a new all-atom force field in XML format, and provide it through ff_add='GLPA.xml' I could not find GLPA on the martini force field. Can I ask where you get this coarse-grained structure?

  1. Atom names defined in an input PDB should be consistent with atom names in your mapping files and atom names in the openMM force field. (The order does not matter).

From an error message, I think you likely hit the first issue. You can provide me with a PDB structure that contains only one GLPA lipid and your mapping file to reduce the size of the file. Best, Siyoung

On Fri, Dec 22, 2023 at 9:33 AM maryammajd @.***> wrote:

Hi again Siyoung,

Thank you for your help and instructions. I was trying to make mapping files for the lipids that were not by default in the map files. They worked quite fine with SM headgroups but when I tried to map some cerebroside lipids, I got the following error. I was wondering if there is any way I can solve this issue. I can also provide the mapping file and the pdb file if required. However, github does not accept those file formats, and therefore I could not attach them.

Error: No template found for residue 1 (GLPA). The set of atoms is similar to DXPC, but it is missing 2 atoms. For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#template

Bests, Maryam

— Reply to this email directly, view it on GitHub https://github.com/ksy141/mstool/issues/1#issuecomment-1867755256, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALSZEC7J5YUZZMFVZKQVRNLYKWK5FAVCNFSM6AAAAABAZ2PLDKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRXG42TKMRVGY . You are receiving this because you modified the open/close state.Message ID: @.***>

maryammajd commented 10 months ago

Hi Siyoung,

A quick explanation about the system. The name of the lipid I wanted to model is DPGS in Martini and I wanted to model it as Glycolipid in the charmm-gui membrane builder method. I checked the .itp file for this lipid, it seems the tails are called CER16 and the head is called BGAL.

About the first requirement, I have changed the name of the residues in the PDB file to GLPA and have added the mapping file that corresponds with this lipid. But as you mentioned, the forcefield cannot detect the model in the openMM forcefield. I would appreciate it if you could guide me to create the xml file for this lipid, as there probably are more lipids in my system that need the xml file. I have also added the PDB (both AA and CG) and Mapping file for this lipid.

GLPA.zip

Bests, Maryam

ksy141 commented 10 months ago

Hi Maryam,

Thank you for explaining your system in detail and sending me your files. I made a new openMM force field for GLPA by stitching ADDG and CER160 parameters in $mstool/FF/charmm36.xml. Please find it attached. Unfortunately, stitching two molecule topologies to make a new molecule is manual and time-consuming. What I did was I made a new file, copied an ADDG topology, copied a CER160 topology, removed atoms that should not be present in GLPA, and connected (added a bond) between one atom of ADDG and one atom of CER160. I wish there were a better way of doing this. I hope the attached file can serve as a template for other lipids in your system.

Also, please review the attached mapping file.

  1. Your file contains H4 C4 C5 C3 O4 chirality; But, I think this one should be revised to H4 C4 C5 O4 C3 chriality.
  2. Your file contains H62 C6 H61 C5 O6 chirality; However, this is not required as H61 and H62 are equivalent, and C6 is not a chiral center.
  3. I added two trans terms (to keep trans isomerism) and two dihedral terms (to keep a chair conformation).

I confirmed the below worked. import mstool mstool.Ungroup('glpa_cg.pdb', 'glpa_aa.pdb', mapping_add='glpa.txt') mstool.REM('glpa_aa.pdb', 'glpa_aa_final.pdb', mapping_add='glpa.txt', ff_add='glpa.xml') mstool.CheckStructure('glpa_aa_final.pdb', mapping_add='glpa.txt')

Best, Siyoung

On Tue, Dec 26, 2023 at 7:20 AM maryammajd @.***> wrote:

Hi Siyoung,

A quick explanation about the system. The name of the lipid I wanted to model is DPGS in Martini and I wanted to model it as Glycolipid in the charmm-gui membrane builder method. I checked the .itp file for this lipid, it seems the tails are called CER16 and the head is called BGAL.

About the first requirement, I have changed the name of the residues in the PDB file to GLPA and have added the mapping file that corresponds with this lipid. But as you mentioned, the forcefield cannot detect the model in the openMM forcefield. I would appreciate it if you could guide me to create the xml file for this lipid, as there probably are more lipids in my system that need the xml file. I have also added the PDB (both AA and CG) and Mapping file for this lipid.

GLPA.zip https://github.com/ksy141/mstool/files/13771416/GLPA.zip

Bests, Maryam

— Reply to this email directly, view it on GitHub https://github.com/ksy141/mstool/issues/1#issuecomment-1869505747, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALSZEC3DYAS4DZOAQTEOA2DYLK6IJAVCNFSM6AAAAABAZ2PLDKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRZGUYDKNZUG4 . You are receiving this because you modified the open/close state.Message ID: @.***>

ksy141 commented 10 months ago

GLPA_new.zip