Modeller logo

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:   mod8v2 - < mutate_model.py modelname respos resname chain > logfile
   8 #    
   9 #     Example: mod8v2 - < 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');