Size: 6024
Comment: Use python; add note for chain-less PDBs.
|
← Revision 12 as of 2021-02-10 20:51:44 ⇥
Size: 6019
Comment: Update for 10.0 class names
|
Deletions are marked like this. | Additions are marked like this. |
Line 11: | Line 11: |
from modeller.optimizers import molecular_dynamics, conjugate_gradients | from modeller.optimizers import MolecularDynamics, ConjugateGradients |
Line 37: | Line 37: |
cg = conjugate_gradients() | cg = ConjugateGradients() |
Line 44: | Line 44: |
md = molecular_dynamics(cap_atom_shift=0.39, md_time_step=4.0, md_return='FINAL') |
md = MolecularDynamics(cap_atom_shift=0.39, md_time_step=4.0, md_return='FINAL') |
Line 60: | Line 60: |
s = selection(mdl1) | s = Selection(mdl1) |
Line 75: | Line 75: |
env = environ(rand_seed=-49837) | env = Environ(rand_seed=-49837) |
Line 93: | Line 93: |
mdl1 = model(env, file=modelname) ali = alignment(env) |
mdl1 = Model(env, file=modelname) ali = Alignment(env) |
Line 98: | Line 98: |
s = selection(mdl1.chains[chain].residues[respos]) | s = Selection(mdl1.chains[chain].residues[respos]) |
Line 120: | Line 120: |
mdl2 = model(env, file=modelname) | mdl2 = Model(env, file=modelname) |
Line 148: | Line 148: |
s = selection(mdl1.chains[chain].residues[respos]) | s = Selection(mdl1.chains[chain].residues[respos]) |
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 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')