search

Acellera / htmd

HTMD: Programming Environment for Molecular Discovery
https://software.acellera.com/docs/latest/htmd/index.html
Other
259 stars 59 forks source link

Add phosphorylation patch #1095

Closed H-EKE closed 5 days ago

H-EKE commented 1 month ago

Hi,

I wanted to add a phosphorylation patch to my peptide, but I haven't been able to do it successfully using HTMD. I was wondering which is the correct way to add a phosphoserine patch. I have attached the the input files and log

from htmd.ui import *
from moleculekit.config import config
topos  = charmm.defaultTopo()+ ['toppar_all36_prot_na_combined_1.rtf']
params = charmm.defaultParam()+['toppar_all36_prot_na_combined_1.inp']

prot = Molecule('peptide.pdb')

prot = autoSegment(prot, sel='protein')
prot.set('segid', 'A', sel='resname CA')
[build.zip](https://github.com/user-attachments/files/17013514/build.zip)

patch = ['patch SP2 P0:23']
molbuilt = charmm.build(prot, topo=topos, param=params, outdir='./build', saltconc=0.15,patches=patch) 
BuildError                                Traceback (most recent call last)
Cell In[29], line 1
----> 1 molbuilt = charmm.build(prot, topo=topos, param=params, outdir='./build', saltconc=0.15,patches=patch) 

File /nfs/htmd/lib/python3.10/site-packages/htmd/builder/charmm.py:394, in build(mol, topo, param, stream, prefix, outdir, caps, ionize, saltconc, saltanion, saltcation, disulfide, regenerate, patches, noregen, aliasresidues, psfgen, execute, _clean)
    392     molbuilt.read(path.join(outdir, "structure.psf"))
    393 else:
--> 394     raise BuildError(
    395         "No structure pdb/psf file was generated. Check {} for errors in building.".format(
    396             logpath
    397         )
    398     )
    400 if ionize:
    401     os.makedirs(path.join(outdir, "pre-ionize"))

BuildError: 'No structure pdb/psf file was generated. Check /home/build/log.txt for errors in building.'

log.txt 

psfgen) add atom failed in patch SP2
ERROR: failed to apply patch
MOLECULE DESTROYED BY FATAL ERROR!  Use resetpsf to start over.
    while executing
"patch SP2 P0:23"
    (file "./build.vmd" line 169)
P0:28 
P0:22

build_1.zip

stefdoerr commented 1 month ago

https://www.researchgate.net/post/Error_in_Phosphorylation_patch2 I think SEP is the correct patch name, not SP2

H-EKE commented 1 month ago

Hi @stefdoerr

Thanks for your fast reply.

When i tried 

prot = autoSegment(prot, sel='protein')
prot.set('segid', 'A', sel='resname CA')
patch = ['patch SEP P0:23']
molbuilt = charmm.build(prot, topo=topos, param=params, outdir='./build', saltconc=0.15,patches=patch)

I still have the same error

psfgen) applying patch SEP to 1 residues
psfgen) unknown patch type SEP
ERROR: failed to apply patch
MOLECULE DESTROYED BY FATAL ERROR!  Use resetpsf to start over.
    while executing
"patch SEP P0:23"
    (file "./build.vmd" line 169)
[peptide.zip](https://github.com/user-attachments/files/17042361/peptide.zip)

P0:28 
P0:22

I thought patch should be either SP1 or SP2 , so one can chose between monoanionic and dianionic 


PRES SP1        -1.00  ! convert serine to monoanionic phosphoserine
                       ! use in a patch statement as follows
! patch to convert serine to phosphoserine
!PATCH SP1 segid resnum SETUP WARN
!
DELE ATOM 1HG1
GROUP
ATOM CB   CT2    -0.08  !
ATOM HB1  HA2     0.09  !
ATOM HB2  HA2     0.09  !
ATOM OG   ON2    -0.62  !maintain NA atom type
ATOM P    P       1.50
ATOM O1P  ON3    -0.82
ATOM O2P  ON3    -0.82
ATOM OT   ON4    -0.68
ATOM HT   HN4     0.34
BOND OG   P    P   OT   OT  HT
BOND P    O1P  P   O2P
BILD CA    CB   OG  P      0.0000  000.00  180.00   000.00   0.0000
BILD CB    OG   P   OT     0.0000  000.00  180.00   000.00   0.0000
BILD OG    P    OT  HT     0.0000  000.00  180.00   000.00   0.0000
BILD OG    OT   *P  O1P    0.0000  000.00 -120.00   000.00   0.0000
BILD OG    OT   *P  O2P    0.0000  000.00  120.00   000.00   0.0000

PRES SP2         -2.00  ! convert serine to dianionic phosphoserine
                        ! use in a patch statement
DELE ATOM 1HG1
GROUP
ATOM CB   CT2    -0.18  !
ATOM HB1  HA2     0.09  !
ATOM HB2  HA2     0.09  !
ATOM OG   ON2    -0.40  !maintain NA atom type
ATOM P    P       1.10
ATOM O1P  ON3    -0.90
ATOM O2P  ON3    -0.90
ATOM OT   ON3    -0.90
BOND OG   P    P   OT  
BOND P    O1P  P   O2P
BILD CA    CB   OG  P      0.0000  000.00  180.00   000.00   0.0000
BILD CB    OG   P   OT     0.0000  000.00  180.00   000.00   0.0000
BILD OG    OT   *P  O1P    0.0000  000.00 -120.00   000.00   0.0000
BILD OG    OT   *P  O2P    0.0000  000.00  120.00   000.00   0.0000

RESI SEP         -1.00  !phosphorylated serine, monoanionic
GROUP   
ATOM N    NH1    -0.47  !     |       
ATOM HN   H       0.31  !  HN-N       
ATOM CA   CT1     0.07  !     |   HB1
ATOM HA   HB1     0.09  !     |   |  
GROUP                   !  HA-CA--CB--OG  O1P
ATOM CB   CT2    -0.08  !     |   |     \ //
ATOM HB1  HA2     0.09  !     |   HB2    P
ATOM HB2  HA2     0.09  !   O=C         / \\ 
                        !     |       O3P  O2P   
                        !              |
                        !             H3T
ATOM OG   ON2    -0.62  !maintain NA atom type
ATOM P    P       1.50
ATOM O1P  ON3    -0.82
ATOM O2P  ON3    -0.82
ATOM O3P  ON4    -0.68
ATOM H3T  HN4     0.34

ATOM C    C       0.51
ATOM O    O      -0.51
BOND CB CA   OG CB  N HN  N  CA   
BOND C  CA   C +N  CA HA  CB HB1   
BOND CB HB2 
BOND OG P    P   O3P  O3P  H3T
BOND P  O1P  P   O2P
DOUBLE   O  C      
IMPR N -C CA HN  C CA +N O   
CMAP -C  N  CA  C   N  CA  C  +N
DONOR HN N   
!DONOR HG1 OG   
ACCEPTOR OG   
ACCEPTOR O C   
IC -C   CA   *N   HN    1.3474 124.3700  180.0000 114.1800  0.9999
IC -C   N    CA   C     1.3474 124.3700  180.0000 105.8100  1.5166
IC N    CA   C    +N    1.4579 105.8100  180.0000 117.7200  1.3448
IC +N   CA   *C   O     1.3448 117.7200  180.0000 120.2500  1.2290
IC CA   C    +N   +CA   1.5166 117.7200  180.0000 124.6300  1.4529
IC N    C    *CA  CB    1.4579 105.8100  124.7500 111.4000  1.5585
IC N    C    *CA  HA    1.4579 105.8100 -115.5600 107.3000  1.0821
IC N    CA   CB   OG    1.4579 114.2800  180.0000 112.4500  1.4341
IC OG   CA   *CB  HB1   1.4341 112.4500  119.3200 108.1000  1.1140
IC OG   CA   *CB  HB2   1.4341 112.4500 -123.8600 110.3800  1.1136
!IC CA  CB   OG   HG1   1.5585 112.4500  165.9600 107.0800  0.9655
IC CA   CB   OG   P     0.0000 000.00    180.00   000.00    0.0000
IC CB   OG   P    O3P   0.0000 000.00    180.00   000.00    0.0000
IC OG   P    O3P  H3T   0.0000 000.00    180.00   000.00    0.0000
IC OG   O3P  *P   O1P   0.0000 000.00   -120.00   000.00    0.0000
IC OG   O3P  *P   O2P   0.0000 000.00    120.00   000.00    0.0000
stefdoerr commented 1 month ago

Ah could it be that it expects the residue to be called SEP before applying the SP1/SP2 patches? Can you try mutating it with mol.mutateResidue() and then building the system? That should recreate the sidechain

H-EKE commented 1 month ago

Hi @stefdoerr,

Thanks for your reply.

I tried this and still didn't work

prot = autoSegment(prot, sel='protein')
prot.set('segid', 'A', sel='resname CA')
prot.mutateResidue("resid 23","SEP")
patch = ['patch SP1 A:23']
molbuilt = charmm.build(prot, topo=topos, param=params, outdir='./build', saltconc=0.15) 

2024-10-08 14:46:15,509 - htmd.builder.charmm - WARNING - Segment P0 consists of a peptide with less than 10 residues. It will not be capped by default. If you want to cap it use the caps argument of charmm.build to manually define caps for all segments
2024-10-08 14:46:15,586 - htmd.builder.charmm - INFO - Writing out segments.
2024-10-08 14:46:15,837 - htmd.builder.charmm - INFO - Starting the build.
2024-10-08 14:46:15,905 - htmd.builder.charmm - INFO - Finished building.

---------------------------------------------------------------------------
BuildError                                Traceback (most recent call last)
Cell In[25], line 5
      3 prot.mutateResidue("resid 23","SEP")
      4 patch = ['patch SP2 A:23']
----> 5 molbuilt = charmm.build(prot, topo=topos, param=params, outdir='./build', saltconc=0.15) 

File /nfs/htmd/lib/python3.10/site-packages/htmd/builder/charmm.py:394, in build(mol, topo, param, stream, prefix, outdir, caps, ionize, saltconc, saltanion, saltcation, disulfide, regenerate, patches, noregen, aliasresidues, psfgen, execute, _clean)
    392     molbuilt.read(path.join(outdir, "structure.psf"))
    393 else:
--> 394     raise BuildError(
    395         "No structure pdb/psf file was generated. Check {} for errors in building.".format(
    396             logpath
    397         )
    398     )
    400 if ionize:
    401     os.makedirs(path.join(outdir, "pre-ionize"))

BuildError: 'No structure pdb/psf file was generated. Check /home/build/log.txt for errors in building.'

The error on the log.txt

ERROR: failed on end of segment
MOLECULE DESTROYED BY FATAL ERROR!  Use resetpsf to start over.
    while executing
"segment P0 {
        pdb segments/segmentP0.pdb
}"
    (file "./build.vmd" line 164)
stefdoerr commented 5 days ago

Hi, since you closed this, did you figure out a solution? Sorry I was not more active, I just managed last week to get psfgen working again for me so that I can check these issues.

H-EKE commented 5 days ago

Hi @stefdoerr ,

I just looked how it works on charmm-gui and try to extrapolate it. Apparently there is no need to change residue name. One just need to import the strucutre and apply the patch directly, and at the end psfgen will regenerate the PTM.

#### Importt protein
proteina = Molecule("protein.pdb")

patch = ['patch SP2 R:396']

charmm.build(mol_solvated,topo=topos,param=params,caps=caps,outdir='./pre-build/',ionize=False,patches=patch)

That worked with different systems that I tried. It would be nice if this information would be added to the webpage, as I was not the only one struggling with this issue