Modeller logo
  • Diff for "Mutate model"
Differences between revisions 2 and 6 (spanning 4 versions)
Revision 2 as of 2005-09-28 19:01:08
Size: 5625
Comment: Add the script itself.
Revision 6 as of 2007-05-07 17:36:56
Size: 5551
Comment:
Deletions are marked like this. Additions are marked like this.
Line 8: Line 8:
from modeller import *
from modeller.optimizers import molecular_dynamics, conjugate_gradients
from modeller.automodel import autosched
Line 10: Line 14:
#    
# Usage: modCVS - < mutate_model.py modelname respos resname chain > logfile
#    
# Example: modCVS - < mutate_model.py 1t29 1699 LEU A > 1t29.log
#
# Usage: mod9v1 - modelname respos resname chain < mutate_model.py > logfile
#
# Example: mod9v1 - 1t29 1699 LEU A < mutate_model.py > 1t29.log
Line 23: Line 27:
def optimize(mdl1): def optimize(atmsel, sched):
Line 25: Line 29:
    for step in range(0,len(mdl1.schedule)):
        mdl1.schedule.step = step + 1
        mdl1
.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)
    for step in sched:
        step.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
Line 29: Line 32:
    refine(mdl1)
    mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)
    refine(atmsel)
    cg = conjugate_gradients()
    cg
.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
Line 35: Line 38:
def refine(mdl1): def refine(atmsel):
Line 37: Line 40:
    md = molecular_dynamics(cap_atom_shift=0.39, md_time_step=4.0,
                            md_return='FINAL')
Line 39: Line 44:
                                (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))):
Line 41: Line 47:
            mdl1.optimize(cap_atom_shift=0.39, md_time_step=4.0, optimization_method=3,
                          md_return='FINAL'
, init_velocities=init_vel, temperature=temp,
                          max_iterations=its, equilibrate=equil)
            md.optimize(atmsel, init_velocities=init_vel, temperature=temp,
                         max_iterations=its, equilibrate=equil)
Line 45: Line 50:
Line 52: Line 56:
   s = selection(mdl1)
Line 53: Line 58:
       rsr.make(restraint_type=typ, aln=aln, spline_on_site=True)        rsr.make(s, restraint_type=typ, aln=aln, spline_on_site=True)
Line 55: Line 60:
       rsr.make(restraint_type=typ+'_dihedral', spline_range=4.0, spline_dx=0.3,
                spline_min_points = 5, aln=aln, spline_on_site=True)
       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)
Line 86: Line 92:
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])
Line 89: Line 95:
mdl1.mutate(residue_type=restyp)
s.mutate(residue_type=restyp)
Line 95: Line 100:
mdl1.generate_topology(ali, sequence=modelname) mdl1.clear_topology()
mdl1.generate_topology(ali[-1])
Line 99: Line 105:
# 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
Line 115: Line 121:
#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).
Line 121: Line 130:
#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
Line 127: Line 137:
mdl1.schedule.make(library_schedule=6) sched = autosched.loop.make_for_model(mdl1)
Line 129: Line 139:
#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)
Line 131: Line 142:
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])
Line 134: Line 145:
mdl1.restraints.pick() mdl1.restraints.pick(s)
Line 136: Line 147:
mdl1.energy() s.energy()
Line 138: Line 149:
mdl1.randomize_xyz(deviation=4.0) s.randomize_xyz(deviation=4.0)
Line 141: Line 152:
optimize(mdl1) optimize(s, sched)
Line 143: Line 154:
#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)
Line 145: Line 157:
optimize(mdl1) optimize(s, sched)
Line 147: Line 159:
mdl1.energy() s.energy()
Line 153: Line 165:
os.remove(modelname+restyp+respos+'.tmp');
os.remove(modelname+restyp+respos+'.tmp')

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.

   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:   mod9v1 - modelname respos resname chain < mutate_model.py > logfile
  12 #
  13 #     Example: mod9v1 - 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 env = environ(rand_seed=-49837)
  67 env.io.hetatm = True
  68 #soft sphere potential
  69 env.edat.dynamic_sphere=False
  70 #lennard-jones potential (more accuate)
  71 env.edat.dynamic_lennard=True
  72 env.edat.contact_shell = 4.0
  73 env.edat.update_dynamic = 0.39
  74 
  75 # Read customized topology file with phosphoserines (or standard one)
  76 env.libs.topology.read(file='$(LIB)/top_heav.lib')
  77 
  78 # Read customized CHARMM parameter library with phosphoserines (or standard one)
  79 env.libs.parameters.read(file='$(LIB)/par.lib')
  80 
  81 
  82 # Read the original PDB file and copy its sequence to the alignment array:
  83 mdl1 = model(env, file=modelname)
  84 ali = alignment(env)
  85 ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)
  86 
  87 #set up the mutate residue selection segment
  88 s = selection(mdl1.chains[chain].residues[respos])
  89 
  90 #perform the mutate residue operation
  91 s.mutate(residue_type=restyp)
  92 #get two copies of the sequence.  A modeller trick to get things set up
  93 ali.append_model(mdl1, align_codes=modelname)
  94 
  95 # Generate molecular topology for mutant
  96 mdl1.clear_topology()
  97 mdl1.generate_topology(ali[-1])
  98 
  99 
 100 # Transfer all the coordinates you can from the template native structure
 101 # to the mutant (this works even if the order of atoms in the native PDB
 102 # file is not standard):
 103 #here we are generating the model by reading the template coordinates
 104 mdl1.transfer_xyz(ali)
 105 
 106 # Build the remaining unknown coordinates
 107 mdl1.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
 108 
 109 #yes model2 is the same file as model1.  It's a modeller trick.
 110 mdl2 = model(env, file=modelname)
 111 
 112 #required to do a transfer_res_numb
 113 #ali.append_model(mdl2, atom_files=modelname, align_codes=modelname)
 114 #transfers from "model 2" to "model 1"
 115 mdl1.res_num_from(mdl2,ali)
 116 
 117 #It is usually necessary to write the mutated sequence out and read it in
 118 #before proceeding, because not all sequence related information about MODEL
 119 #is changed by this command (e.g., internal coordinates, charges, and atom
 120 #types and radii are not updated).
 121 
 122 mdl1.write(file=modelname+restyp+respos+'.tmp')
 123 mdl1.read(file=modelname+restyp+respos+'.tmp')
 124 
 125 #set up restraints before computing energy
 126 #we do this a second time because the model has been written out and read in,
 127 #clearing the previously set restraints
 128 make_restraints(mdl1, ali)
 129 
 130 #a non-bonded pair has to have at least as many selected atoms
 131 mdl1.env.edat.nonbonded_sel_atoms=1
 132 
 133 sched = autosched.loop.make_for_model(mdl1)
 134 
 135 #only optimize the selected residue (in first pass, just atoms in selected
 136 #residue, in second pass, include nonbonded neighboring atoms)
 137 #set up the mutate residue selection segment
 138 s = selection(mdl1.chains[chain].residues[respos])
 139 
 140 mdl1.restraints.unpick_all()
 141 mdl1.restraints.pick(s)
 142 
 143 s.energy()
 144 
 145 s.randomize_xyz(deviation=4.0)
 146 
 147 mdl1.env.edat.nonbonded_sel_atoms=2
 148 optimize(s, sched)
 149 
 150 #feels environment (energy computed on pairs that have at least one member
 151 #in the selected)
 152 mdl1.env.edat.nonbonded_sel_atoms=1
 153 optimize(s, sched)
 154 
 155 s.energy()
 156 
 157 #give a proper name
 158 mdl1.write(file=modelname+restyp+respos+'.pdb')
 159 
 160 #delete the temporary file
 161 os.remove(modelname+restyp+respos+'.tmp')