No edit summary
(Add to Examples category)
 
(12 intermediate revisions by the same user not shown)
Line 2: Line 2:
The script below takes a given PDB file, and mutates a single residue. The residue's position is then optimized, and the unoptimized and optimized energies are reported.
The script below takes a given PDB file, and mutates a single residue. The residue's position is then optimized, and the unoptimized and optimized energies are reported.


Note that this script if run multiple times will produce the same model each time, because Modeller is deterministic. If you want to build multiple models, change the value of <code>rand_seed</code> (see comments in the script) each time. This may be useful if some models, for example, cannot be optimized due to steric clashes.


<pre><nowiki>#!python
 
<syntaxhighlight lang="python">
import sys
import sys
import os
import os
from modeller import *
from modeller.optimizers import MolecularDynamics, ConjugateGradients
from modeller.automodel import autosched


#
#
#  mutate_model.py
#  mutate_model.py
#  
#
#    Usage:  modCVS - < mutate_model.py modelname respos resname chain > logfile
#    Usage:  python mutate_model.py modelname respos resname chain > logfile
#  
#
#    Example: modCVS - < mutate_model.py 1t29 1699 LEU A > 1t29.log
#    Example: python mutate_model.py 1t29 1699 LEU A > 1t29.log
#
#
#
#
Line 19: Line 25:
#  The conformation of the mutant sidechain is optimized by conjugate gradient and
#  The conformation of the mutant sidechain is optimized by conjugate gradient and
#  refined using some MD.
#  refined using some MD.
#
#  Note: if the model has no chain identifier, specify "" for the chain argument.
#
#




def optimize(mdl1):
def optimize(atmsel, sched):
     #conjugate gradient
     #conjugate gradient
     for step in range(0,len(mdl1.schedule)):
     for step in sched:
         mdl1.schedule.step = step + 1
         step.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
        mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)
     #md
     #md
     refine(mdl1)
     refine(atmsel)
     mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)
     cg = ConjugateGradients()
 
    cg.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)




#molecular dynamics
#molecular dynamics
def refine(mdl1):
def refine(atmsel):
     # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
     # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
    md = MolecularDynamics(cap_atom_shift=0.39, md_time_step=4.0,
                          md_return='FINAL')
     init_vel = True
     init_vel = True
     for (its, equil, temps) in ((200, 20, (150.0, 250.0, 400.0, 700.0, 1000.0)),
     for (its, equil, temps) in ((200, 20, (150.0, 250.0, 400.0, 700.0, 1000.0)),
                                 (200, 600, (1000.0, 800.0, 600.0, 500.0, 400.0, 300.0))):
                                 (200, 600,
                                (1000.0, 800.0, 600.0, 500.0, 400.0, 300.0))):
         for temp in temps:
         for temp in temps:
             mdl1.optimize(cap_atom_shift=0.39, md_time_step=4.0, optimization_method=3,
             md.optimize(atmsel, init_velocities=init_vel, temperature=temp,
                          md_return='FINAL', init_velocities=init_vel, temperature=temp,
                        max_iterations=its, equilibrate=equil)
                          max_iterations=its, equilibrate=equil)
             init_vel = False
             init_vel = False




Line 51: Line 59:
   rsr = mdl1.restraints
   rsr = mdl1.restraints
   rsr.clear()
   rsr.clear()
  s = Selection(mdl1)
   for typ in ('stereo', 'phi-psi_binormal'):
   for typ in ('stereo', 'phi-psi_binormal'):
       rsr.make(restraint_type=typ, aln=aln, spline_on_site=True)
       rsr.make(s, restraint_type=typ, aln=aln, spline_on_site=True)
   for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
   for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
       rsr.make(restraint_type=typ+'_dihedral', spline_range=4.0, spline_dx=0.3,
       rsr.make(s, restraint_type=typ+'_dihedral', spline_range=4.0,
                spline_min_points = 5, aln=aln, spline_on_site=True)
                spline_dx=0.3, spline_min_points = 5, aln=aln,
                spline_on_site=True)


#first argument
#first argument
Line 63: Line 73:
log.verbose()
log.verbose()


env = environ(rand_seed=-49837)
# Set a different value for rand_seed to get a different final model
env = Environ(rand_seed=-49837)
 
env.io.hetatm = True
env.io.hetatm = True
#soft sphere potential
#soft sphere potential
env.edat.dynamic_sphere=False
env.edat.dynamic_sphere=False
#lennard-jones potential (more accuate)
#lennard-jones potential (more accurate)
env.edat.dynamic_lennard=True
env.edat.dynamic_lennard=True
env.edat.contact_shell = 4.0
env.edat.contact_shell = 4.0
Line 80: Line 92:


# Read the original PDB file and copy its sequence to the alignment array:
# Read the original PDB file and copy its sequence to the alignment array:
mdl1 = model(env, file=modelname)
mdl1 = Model(env, file=modelname)
ali = alignment(env)
ali = Alignment(env)
ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)
ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)


#set up the mutate residue selection segment
#set up the mutate residue selection segment
mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',selection_segment = (respos+':'+chain, respos+':'+chain), atom_types = 'ALL', selection_status = 'INITIALIZE')
s = Selection(mdl1.chains[chain].residues[respos])


#perform the mutate residue operation
#perform the mutate residue operation
mdl1.mutate(residue_type=restyp)
s.mutate(residue_type=restyp)
 
#get two copies of the sequence.  A modeller trick to get things set up
#get two copies of the sequence.  A modeller trick to get things set up
ali.append_model(mdl1, align_codes=modelname)
ali.append_model(mdl1, align_codes=modelname)


# Generate molecular topology for mutant
# Generate molecular topology for mutant
mdl1.generate_topology(ali, sequence=modelname)
mdl1.clear_topology()
mdl1.generate_topology(ali[-1])




# Transfer all the coordinates you can from the template native structure
# Transfer all the coordinates you can from the template native structure
# to the mutant (this works even if the order of atoms in the native PDB  
# to the mutant (this works even if the order of atoms in the native PDB
# file is not standard):
# file is not standard):
#here we are generating the model by reading the template coordinates
#here we are generating the model by reading the template coordinates
Line 107: Line 119:


#yes model2 is the same file as model1.  It's a modeller trick.
#yes model2 is the same file as model1.  It's a modeller trick.
mdl2 = model(env, file=modelname)
mdl2 = Model(env, file=modelname)


#required to do a transfer_res_numb
#required to do a transfer_res_numb
Line 114: Line 126:
mdl1.res_num_from(mdl2,ali)
mdl1.res_num_from(mdl2,ali)


#It is usually necessary to write the mutated sequence out and read it in before proceeding, because not all sequence related information about MODEL is changed by this command (e.g., internal coordinates, charges, and atom types and radii are not updated).
#It is usually necessary to write the mutated sequence out and read it in
#before proceeding, because not all sequence related information about MODEL
#is changed by this command (e.g., internal coordinates, charges, and atom
#types and radii are not updated).


mdl1.write(file=modelname+restyp+respos+'.tmp')
mdl1.write(file=modelname+restyp+respos+'.tmp')
Line 120: Line 135:


#set up restraints before computing energy
#set up restraints before computing energy
#we do this a second time because the model has been written out and read in, clearing the previously set restraints
#we do this a second time because the model has been written out and read in,
#clearing the previously set restraints
make_restraints(mdl1, ali)
make_restraints(mdl1, ali)


Line 126: Line 142:
mdl1.env.edat.nonbonded_sel_atoms=1
mdl1.env.edat.nonbonded_sel_atoms=1


mdl1.schedule.make(library_schedule=6)
sched = autosched.loop.make_for_model(mdl1)


#only optimize the selected residue (in first pass, just atoms in selected residue, in second pass, include nonbonded neighboring atoms)
#only optimize the selected residue (in first pass, just atoms in selected
#residue, in second pass, include nonbonded neighboring atoms)
#set up the mutate residue selection segment
#set up the mutate residue selection segment
mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',selection_segment = (respos+':'+chain, respos+':'+chain), atom_types = 'ALL', selection_status = 'INITIALIZE')
s = Selection(mdl1.chains[chain].residues[respos])


mdl1.restraints.unpick_all()
mdl1.restraints.unpick_all()
mdl1.restraints.pick()
mdl1.restraints.pick(s)


mdl1.energy()
s.energy()


mdl1.randomize_xyz(deviation=4.0)
s.randomize_xyz(deviation=4.0)


mdl1.env.edat.nonbonded_sel_atoms=2
mdl1.env.edat.nonbonded_sel_atoms=2
optimize(mdl1)
optimize(s, sched)


#feels environment (energy computed on pairs that have at least one member in the selected)
#feels environment (energy computed on pairs that have at least one member
#in the selected)
mdl1.env.edat.nonbonded_sel_atoms=1
mdl1.env.edat.nonbonded_sel_atoms=1
optimize(mdl1)
optimize(s, sched)


mdl1.energy()
s.energy()


#give a proper name
#give a proper name
Line 152: Line 170:


#delete the temporary file
#delete the temporary file
os.remove(modelname+restyp+respos+'.tmp');
os.remove(modelname+restyp+respos+'.tmp')
 
</syntaxhighlight>
</nowiki></pre>
 


[[Category:Examples]]

Latest revision as of 21:17, 16 August 2022

The script below takes a given PDB file, and mutates a single residue. The residue's position is then optimized, and the unoptimized and optimized energies are reported.

Note that this script if run multiple times will produce the same model each time, because Modeller is deterministic. If you want to build multiple models, change the value of rand_seed (see comments in the script) each time. This may be useful if some models, for example, cannot be optimized due to steric clashes.


import sys
import os

from modeller import *
from modeller.optimizers import MolecularDynamics, ConjugateGradients
from modeller.automodel import autosched

#
#  mutate_model.py
#
#     Usage:   python mutate_model.py modelname respos resname chain > logfile
#
#     Example: python mutate_model.py 1t29 1699 LEU A > 1t29.log
#
#
#  Creates a single in silico point mutation to sidechain type and at residue position
#  input by the user, in the structure whose file is modelname.pdb
#  The conformation of the mutant sidechain is optimized by conjugate gradient and
#  refined using some MD.
#
#  Note: if the model has no chain identifier, specify "" for the chain argument.
#


def optimize(atmsel, sched):
    #conjugate gradient
    for step in sched:
        step.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
    #md
    refine(atmsel)
    cg = ConjugateGradients()
    cg.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)


#molecular dynamics
def refine(atmsel):
    # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
    md = MolecularDynamics(cap_atom_shift=0.39, md_time_step=4.0,
                           md_return='FINAL')
    init_vel = True
    for (its, equil, temps) in ((200, 20, (150.0, 250.0, 400.0, 700.0, 1000.0)),
                                (200, 600,
                                 (1000.0, 800.0, 600.0, 500.0, 400.0, 300.0))):
        for temp in temps:
            md.optimize(atmsel, init_velocities=init_vel, temperature=temp,
                         max_iterations=its, equilibrate=equil)
            init_vel = False


#use homologs and dihedral library for dihedral angle restraints
def make_restraints(mdl1, aln):
   rsr = mdl1.restraints
   rsr.clear()
   s = Selection(mdl1)
   for typ in ('stereo', 'phi-psi_binormal'):
       rsr.make(s, restraint_type=typ, aln=aln, spline_on_site=True)
   for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
       rsr.make(s, restraint_type=typ+'_dihedral', spline_range=4.0,
                spline_dx=0.3, spline_min_points = 5, aln=aln,
                spline_on_site=True)

#first argument
modelname, respos, restyp, chain, = sys.argv[1:]


log.verbose()

# Set a different value for rand_seed to get a different final model
env = Environ(rand_seed=-49837)

env.io.hetatm = True
#soft sphere potential
env.edat.dynamic_sphere=False
#lennard-jones potential (more accurate)
env.edat.dynamic_lennard=True
env.edat.contact_shell = 4.0
env.edat.update_dynamic = 0.39

# Read customized topology file with phosphoserines (or standard one)
env.libs.topology.read(file='$(LIB)/top_heav.lib')

# Read customized CHARMM parameter library with phosphoserines (or standard one)
env.libs.parameters.read(file='$(LIB)/par.lib')


# Read the original PDB file and copy its sequence to the alignment array:
mdl1 = Model(env, file=modelname)
ali = Alignment(env)
ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)

#set up the mutate residue selection segment
s = Selection(mdl1.chains[chain].residues[respos])

#perform the mutate residue operation
s.mutate(residue_type=restyp)
#get two copies of the sequence.  A modeller trick to get things set up
ali.append_model(mdl1, align_codes=modelname)

# Generate molecular topology for mutant
mdl1.clear_topology()
mdl1.generate_topology(ali[-1])


# Transfer all the coordinates you can from the template native structure
# to the mutant (this works even if the order of atoms in the native PDB
# file is not standard):
#here we are generating the model by reading the template coordinates
mdl1.transfer_xyz(ali)

# Build the remaining unknown coordinates
mdl1.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

#yes model2 is the same file as model1.  It's a modeller trick.
mdl2 = Model(env, file=modelname)

#required to do a transfer_res_numb
#ali.append_model(mdl2, atom_files=modelname, align_codes=modelname)
#transfers from "model 2" to "model 1"
mdl1.res_num_from(mdl2,ali)

#It is usually necessary to write the mutated sequence out and read it in
#before proceeding, because not all sequence related information about MODEL
#is changed by this command (e.g., internal coordinates, charges, and atom
#types and radii are not updated).

mdl1.write(file=modelname+restyp+respos+'.tmp')
mdl1.read(file=modelname+restyp+respos+'.tmp')

#set up restraints before computing energy
#we do this a second time because the model has been written out and read in,
#clearing the previously set restraints
make_restraints(mdl1, ali)

#a non-bonded pair has to have at least as many selected atoms
mdl1.env.edat.nonbonded_sel_atoms=1

sched = autosched.loop.make_for_model(mdl1)

#only optimize the selected residue (in first pass, just atoms in selected
#residue, in second pass, include nonbonded neighboring atoms)
#set up the mutate residue selection segment
s = Selection(mdl1.chains[chain].residues[respos])

mdl1.restraints.unpick_all()
mdl1.restraints.pick(s)

s.energy()

s.randomize_xyz(deviation=4.0)

mdl1.env.edat.nonbonded_sel_atoms=2
optimize(s, sched)

#feels environment (energy computed on pairs that have at least one member
#in the selected)
mdl1.env.edat.nonbonded_sel_atoms=1
optimize(s, sched)

s.energy()

#give a proper name
mdl1.write(file=modelname+restyp+respos+'.pdb')

#delete the temporary file
os.remove(modelname+restyp+respos+'.tmp')