Modeller logo
  • Diff for "Mutate model"
Differences between revisions 5 and 10 (spanning 5 versions)
Revision 5 as of 2007-02-13 19:40:44
Size: 5551
Comment: Update for 9v1.
Revision 10 as of 2008-07-09 01:53:42
Size: 5946
Comment: Add note about changing rand_seed to change models.
Deletions are marked like this. Additions are marked like this.
Line 2: Line 2:

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.
Line 15: Line 17:
# Usage: mod9v1 - < mutate_model.py modelname respos resname chain > logfile # Usage: mod9v2 - modelname respos resname chain < mutate_model.py > logfile
Line 17: Line 19:
# Example: mod9v1 - < mutate_model.py 1t29 1699 LEU A > 1t29.log # Example: mod9v2 - 1t29 1699 LEU A < mutate_model.py > 1t29.log
Line 70: Line 72:
# Set a different value for rand_seed to get a different final model
Line 71: Line 74:
Line 74: Line 78:
#lennard-jones potential (more accuate) #lennard-jones potential (more accurate)

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.

   1 import sys
   2 import os
   3 
   4 from modeller import *
   5 from modeller.optimizers import molecular_dynamics, conjugate_gradients
   6 from modeller.automodel import autosched
   7 
   8 #
   9 #  mutate_model.py
  10 #
  11 #     Usage:   mod9v2 - modelname respos resname chain < mutate_model.py > logfile
  12 #
  13 #     Example: mod9v2 - 1t29 1699 LEU A < mutate_model.py > 1t29.log
  14 #
  15 #
  16 #  Creates a single in silico point mutation to sidechain type and at residue position
  17 #  input by the user, in the structure whose file is modelname.pdb
  18 #  The conformation of the mutant sidechain is optimized by conjugate gradient and
  19 #  refined using some MD.
  20 #
  21 
  22 
  23 def optimize(atmsel, sched):
  24     #conjugate gradient
  25     for step in sched:
  26         step.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
  27     #md
  28     refine(atmsel)
  29     cg = conjugate_gradients()
  30     cg.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
  31 
  32 
  33 #molecular dynamics
  34 def refine(atmsel):
  35     # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
  36     md = molecular_dynamics(cap_atom_shift=0.39, md_time_step=4.0,
  37                             md_return='FINAL')
  38     init_vel = True
  39     for (its, equil, temps) in ((200, 20, (150.0, 250.0, 400.0, 700.0, 1000.0)),
  40                                 (200, 600,
  41                                  (1000.0, 800.0, 600.0, 500.0, 400.0, 300.0))):
  42         for temp in temps:
  43             md.optimize(atmsel, init_velocities=init_vel, temperature=temp,
  44                          max_iterations=its, equilibrate=equil)
  45             init_vel = False
  46 
  47 
  48 #use homologs and dihedral library for dihedral angle restraints
  49 def make_restraints(mdl1, aln):
  50    rsr = mdl1.restraints
  51    rsr.clear()
  52    s = selection(mdl1)
  53    for typ in ('stereo', 'phi-psi_binormal'):
  54        rsr.make(s, restraint_type=typ, aln=aln, spline_on_site=True)
  55    for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
  56        rsr.make(s, restraint_type=typ+'_dihedral', spline_range=4.0,
  57                 spline_dx=0.3, spline_min_points = 5, aln=aln,
  58                 spline_on_site=True)
  59 
  60 #first argument
  61 modelname, respos, restyp, chain, = sys.argv[1:]
  62 
  63 
  64 log.verbose()
  65 
  66 # Set a different value for rand_seed to get a different final model
  67 env = environ(rand_seed=-49837)
  68 
  69 env.io.hetatm = True
  70 #soft sphere potential
  71 env.edat.dynamic_sphere=False
  72 #lennard-jones potential (more accurate)
  73 env.edat.dynamic_lennard=True
  74 env.edat.contact_shell = 4.0
  75 env.edat.update_dynamic = 0.39
  76 
  77 # Read customized topology file with phosphoserines (or standard one)
  78 env.libs.topology.read(file='$(LIB)/top_heav.lib')
  79 
  80 # Read customized CHARMM parameter library with phosphoserines (or standard one)
  81 env.libs.parameters.read(file='$(LIB)/par.lib')
  82 
  83 
  84 # Read the original PDB file and copy its sequence to the alignment array:
  85 mdl1 = model(env, file=modelname)
  86 ali = alignment(env)
  87 ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)
  88 
  89 #set up the mutate residue selection segment
  90 s = selection(mdl1.chains[chain].residues[respos])
  91 
  92 #perform the mutate residue operation
  93 s.mutate(residue_type=restyp)
  94 #get two copies of the sequence.  A modeller trick to get things set up
  95 ali.append_model(mdl1, align_codes=modelname)
  96 
  97 # Generate molecular topology for mutant
  98 mdl1.clear_topology()
  99 mdl1.generate_topology(ali[-1])
 100 
 101 
 102 # Transfer all the coordinates you can from the template native structure
 103 # to the mutant (this works even if the order of atoms in the native PDB
 104 # file is not standard):
 105 #here we are generating the model by reading the template coordinates
 106 mdl1.transfer_xyz(ali)
 107 
 108 # Build the remaining unknown coordinates
 109 mdl1.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
 110 
 111 #yes model2 is the same file as model1.  It's a modeller trick.
 112 mdl2 = model(env, file=modelname)
 113 
 114 #required to do a transfer_res_numb
 115 #ali.append_model(mdl2, atom_files=modelname, align_codes=modelname)
 116 #transfers from "model 2" to "model 1"
 117 mdl1.res_num_from(mdl2,ali)
 118 
 119 #It is usually necessary to write the mutated sequence out and read it in
 120 #before proceeding, because not all sequence related information about MODEL
 121 #is changed by this command (e.g., internal coordinates, charges, and atom
 122 #types and radii are not updated).
 123 
 124 mdl1.write(file=modelname+restyp+respos+'.tmp')
 125 mdl1.read(file=modelname+restyp+respos+'.tmp')
 126 
 127 #set up restraints before computing energy
 128 #we do this a second time because the model has been written out and read in,
 129 #clearing the previously set restraints
 130 make_restraints(mdl1, ali)
 131 
 132 #a non-bonded pair has to have at least as many selected atoms
 133 mdl1.env.edat.nonbonded_sel_atoms=1
 134 
 135 sched = autosched.loop.make_for_model(mdl1)
 136 
 137 #only optimize the selected residue (in first pass, just atoms in selected
 138 #residue, in second pass, include nonbonded neighboring atoms)
 139 #set up the mutate residue selection segment
 140 s = selection(mdl1.chains[chain].residues[respos])
 141 
 142 mdl1.restraints.unpick_all()
 143 mdl1.restraints.pick(s)
 144 
 145 s.energy()
 146 
 147 s.randomize_xyz(deviation=4.0)
 148 
 149 mdl1.env.edat.nonbonded_sel_atoms=2
 150 optimize(s, sched)
 151 
 152 #feels environment (energy computed on pairs that have at least one member
 153 #in the selected)
 154 mdl1.env.edat.nonbonded_sel_atoms=1
 155 optimize(s, sched)
 156 
 157 s.energy()
 158 
 159 #give a proper name
 160 mdl1.write(file=modelname+restyp+respos+'.pdb')
 161 
 162 #delete the temporary file
 163 os.remove(modelname+restyp+respos+'.tmp')