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