Modeller logo
  • Diff for "Mutate model"
Differences between revisions 3 and 6 (spanning 3 versions)
Revision 3 as of 2005-09-28 19:06:20
Size: 5707
Comment: Tidy up formatting.
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 91: Line 95:
mdl1.mutate(residue_type=restyp)
s.mutate(residue_type=restyp)
Line 97: Line 100:
mdl1.generate_topology(ali, sequence=modelname) mdl1.clear_topology()
mdl1.generate_topology(ali[-1])
Line 101: 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 133: Line 137:
mdl1.schedule.make(library_schedule=6) sched = autosched.loop.make_for_model(mdl1)
Line 138: 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 143: Line 145:
mdl1.restraints.pick() mdl1.restraints.pick(s)
Line 145: Line 147:
mdl1.energy() s.energy()
Line 147: Line 149:
mdl1.randomize_xyz(deviation=4.0) s.randomize_xyz(deviation=4.0)
Line 150: Line 152:
optimize(mdl1) optimize(s, sched)
Line 155: Line 157:
optimize(mdl1) optimize(s, sched)
Line 157: Line 159:
mdl1.energy() s.energy()
Line 163: 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')