openmm / pdbfixer

PDBFixer fixes problems in PDB files
Other
453 stars 114 forks source link

how to handle non standard amino acids e.g. PTR #218

Open michelsanner opened 3 years ago

michelsanner commented 3 years ago

I have a kinase with the phosphorylated tyrosine PTR which I would like to keep

my first attempt was to fixer1 = PDBFixer(filename="mypdb.pdb") fixer1.findMissingResidues() fixer1.findMissingAtoms() fixer1.addMissingAtoms() fixer1.addMissingHydrogens(7.4) fixer1.removeHeterogens(keepWater=False)

but that last call removed the PTR residue (even though the ATOM records are ATOM not HETATM)

so I tried without calling fixer1.removeHeterogens() but then the following code

from simtk.openmm.app import PDBFile, ForceField, GBn2 forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml') system = forcefield.createSystem(fixer1.topology, implicitSolvent=GBn2, soluteDielectric=1.0, solventDielectric=80.0)

fails with: raise ValueError('No template found for residue %d (%s). %s' % (res.index+1, res.name, _findMatchErrors(self, res))) ValueError: No template found for residue 443 (GLN). The set of atoms matches GLN, but the bonds are different.

and I can see that GLN443 is missing the peptide bond to PTR 444 and PTR 44 has only 3 bonds:

for b in ptr443.bonds(): print(b) ## missing C - N (443) Bond(<Atom 7068 (C) of chain 0 residue 443 (PTR)>, <Atom 7082 (N) of chain 0 residue 444 (GLN)>) Bond(<Atom 7068 (C) of chain 0 residue 443 (PTR)>, <Atom 7082 (N) of chain 0 residue 444 (GLN)>) Bond(<Atom 7068 (C) of chain 0 residue 443 (PTR)>, <Atom 7082 (N) of chain 0 residue 444 (GLN)>)

so I tried to use replaceNonstandardResidues() to change PTR to TYR fixer1 = PDBFixer(filename=recFilename) fixer1.findMissingResidues() fixer1.findMissingAtoms() fixer1.findNonstandardResidues() fixer1.addMissingAtoms() fixer1.addMissingHydrogens(7.4) fixer1.replaceNonstandardResidues() with open('PTRtoTYR.pdb', 'w') as outfile: PDBFile.writeFile(fixer1.topology, fixer1.positions, file=outfile, keepIds=True)

inspection of PTRtoTYR.pdb still show residue PTR

What am I missing here? what is the proper way to prepare this while for amber while keeping the PTR residue ?

thanks

-Michel

peastman commented 3 years ago

Do you have CONECT records for the nonstandard residue? A PDB file needs to have them for anything except standard amino acids, nucleotides, and water.

michelsanner commented 3 years ago

Hi Peter

Thanks for the rapid reply. Aafter adding the CONECT records for PTR it get

system = forcefield.createSystem(fixer1.topology, implicitSolvent=GBn2, ... soluteDielectric=1.0, solventDielectric=80.0) Traceback (most recent call last): File "", line 2, in File "/home/sanner/local/miniconda3/envs/py37/lib/python3.7/site-packages/simtk/openmm/app/forcefield.py", line 1177, in createSystem templateForResidue = self._matchAllResiduesToTemplates(data, topology, residueTemplates, ignoreExternalBonds) File "/home/sanner/local/miniconda3/envs/py37/lib/python3.7/site-packages/simtk/openmm/app/forcefield.py", line 1391, in _matchAllResiduesToTemplates raise ValueError('No template found for residue %d (%s). %s' % (res.index+1, res.name, _findMatchErrors(self, res))) ValueError: No template found for residue 444 (PTR). The set of atoms is similar to DC, but it is missing 14 atoms.

using the following code on the attached pdb file

from pdbfixer import PDBFixer from simtk.openmm.app import PDBFile

fixer1 = PDBFixer(filename='../2src_20pad.pdb') fixer1.findMissingResidues() # required to be able to call findMissingAtoms fixer1.findMissingAtoms() fixer1.addMissingAtoms() fixer1.addMissingHydrogens(7.4)

from simtk.openmm.app import PDBFile, ForceField, GBn2 forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml') system = forcefield.createSystem(fixer1.topology, implicitSolvent=GBn2, soluteDielectric=1.0, solventDielectric=80.0)

[-Michel 2src_20pad.pdb.zip

peastman commented 3 years ago

The Amber14 force field only contains parameters for standard amino acids. If you want to simulate a nonstandard one, you'll need to provide your own force field parameters for it. Are you sure you really want to do that?

michelsanner commented 3 years ago

Well, as you know that In kinases the phosphorylation of tyrosine residues dictates activation and deactivation of which entails domain rearrangements. If I mutate the PTR to TYR I'd be running MD on active form using a crystal structure of the inactive form. I am not sure how much value there will be in that. The goal of the MD run is to study the stability of a peptide bound to the inactive form.

I recall reading a post from John Chodera specifically about this /export/export/people/sanner/python/collab/parang/tmp where he mentioned "For example, suppose Im working with an AMBER variant that has the phosphorylated kinase residues PTR, and want to keep these?"

Amber Parameters have been reported for phosphorylated amaino acids https://pubmed.ncbi.nlm.nih.gov/16240095/ http://amber.manchester.ac.uk/ # has the parameters .off files

How hard would it be in extend OpenMM with these templates ?

thanks

jchodera commented 3 years ago

The Amber14 force field only contains parameters for standard amino acids. If you want to simulate a nonstandard one, you'll need to provide your own force field parameters for it. Are you sure you really want to do that?

We just added the amber/phosaa14SB.xml file to the new openmmforcefields release! This should be out on conda-forge in a day or two.

You can also try amber/phosaa10.xml in the meantime if you're using older force fields.