Modeller logo
  • Diff for "Mutate model"
Differences between revisions 2 and 3
Revision 2 as of 2005-09-28 19:01:08
Size: 5625
Comment: Add the script itself.
Revision 3 as of 2005-09-28 19:06:20
Size: 5707
Comment: Tidy up formatting.
Deletions are marked like this. Additions are marked like this.
Line 86: Line 86:
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') 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')
Line 115: Line 117:
#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 126:
#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 129: Line 135:
#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 138:
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') 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')
Line 143: Line 152:
#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)

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