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