1o9u.
pdb (renum_1, water & ligand remove)
1. Extract the residues sequence by using the following script -
from modeller import *
# Get the sequence of the 1o9u PDB file, and write to an alignment file
code = '1o9u'
e = Environ()
m = Model(e, file=code)
aln = Alignment(e)
aln.append_model(m, align_codes=code)
aln.write(file=code+'.pir')
Result-
2. Generate the alignment in between target and templete residues
from modeller import *
env = Environ()
aln = Alignment(env)
mdl = Model(env, file='1o9u', model_segment=('FIRST:@','END:'))
aln.append_model(mdl, align_codes='1o9u', atom_files='1o9u.pdb')
aln.append(file='1o9u_bjoin.pir', align_codes='1o9u_bjoin')
aln.align2d()
aln.write(file='alignment.ali', alignment_format='PIR')
Result-
3. Model only the break region of the protein chain
# Comparative modeling by the AutoModel class
#
# Demonstrates how to refine only a part of the model.
#
# You may want to use the more exhaustive "loop" modeling routines instead.
#
from modeller import *
from modeller.automodel import * # Load the AutoModel class
log.verbose()
# Override the 'select_atoms' routine in the 'AutoModel' class:
# (To build an all-hydrogen model, derive from AllHModel rather than AutoModel
# here.)
class MyModel(AutoModel):
def select_atoms(self):
# Select residues 1 and 2 in chain A (PDB numbering)
return Selection(self.residue_range('181:A', '183:A'))
# Residues 4, 6, 10 in chain A:
# return Selection(self.residues['4:A'], self.residues['6:A'],
# self.residues['10:A'])
# All residues except 1-5 in chain A:
# return Selection(self) - Selection(self.residue_range('1:A', '5:A'))
env = Environ()
# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']
# selected atoms do not feel the neighborhood
env.edat.nonbonded_sel_atoms = 2
# Be sure to use 'MyModel' rather than 'AutoModel' here!
a = MyModel(env,
alnfile = 'alignment.ali', # alignment filename
knowns = '1o9u', # codes of the templates
sequence = '1o9u_bjoin') # code of the target
a.starting_model= 3 # index of the first model
a.ending_model = 3 # index of the last model
# (determines how many models to calculate)
a.make() # do comparative modeling
Result-
1o9u_bjoin.B99990003.pdb
4. Mutate selected residue (V232G)
import sys
import os
from modeller import *
from modeller.optimizers import MolecularDynamics, ConjugateGradients
from modeller.automodel import autosched
#
# mutate_model.py
#
# Usage: python mutate_model.py modelname respos resname chain > logfile
#
# Example: python mutate_model.py 1t29 1699 LEU A > 1t29.log
#
#
# Creates a single in silico point mutation to sidechain type and at residue position
# input by the user, in the structure whose file is modelname.pdb
# The conformation of the mutant sidechain is optimized by conjugate gradient and
# refined using some MD.
#
# Note: if the model has no chain identifier, specify "" for the chain argument.
#
def optimize(atmsel, sched):
#conjugate gradient
for step in sched:
step.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
#md
refine(atmsel)
cg = ConjugateGradients()
cg.optimize(atmsel, max_iterations=200, min_atom_shift=0.001)
#molecular dynamics
def refine(atmsel):
# at T=1000, max_atom_shift for 4fs is cca 0.15 A.
md = MolecularDynamics(cap_atom_shift=0.39, md_time_step=4.0,
md_return='FINAL')
init_vel = True
for (its, equil, temps) in ((200, 20, (150.0, 250.0, 400.0, 700.0, 1000.0)),
(200, 600,
(1000.0, 800.0, 600.0, 500.0, 400.0, 300.0))):
for temp in temps:
md.optimize(atmsel, init_velocities=init_vel, temperature=temp,
max_iterations=its, equilibrate=equil)
init_vel = False
#use homologs and dihedral library for dihedral angle restraints
def make_restraints(mdl1, aln):
rsr = mdl1.restraints
rsr.clear()
s = Selection(mdl1)
for typ in ('stereo', 'phi-psi_binormal'):
rsr.make(s, restraint_type=typ, aln=aln, spline_on_site=True)
for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
rsr.make(s, restraint_type=typ+'_dihedral', spline_range=4.0,
spline_dx=0.3, spline_min_points = 5, aln=aln,
spline_on_site=True)
#first argument
modelname, respos, restyp, chain, = sys.argv[1:]
log.verbose()
# Set a different value for rand_seed to get a different final model
env = Environ(rand_seed=-49837)
env.io.hetatm = True
#soft sphere potential
env.edat.dynamic_sphere=False
#lennard-jones potential (more accurate)
env.edat.dynamic_lennard=True
env.edat.contact_shell = 4.0
env.edat.update_dynamic = 0.39
# Read customized topology file with phosphoserines (or standard one)
env.libs.topology.read(file='$(LIB)/top_heav.lib')
# Read customized CHARMM parameter library with phosphoserines (or standard one)
env.libs.parameters.read(file='$(LIB)/par.lib')
# Read the original PDB file and copy its sequence to the alignment array:
mdl1 = Model(env, file=modelname)
ali = Alignment(env)
ali.append_model(mdl1, atom_files=modelname, align_codes=modelname)
#set up the mutate residue selection segment
s = Selection(mdl1.chains[chain].residues[respos])
#perform the mutate residue operation
s.mutate(residue_type=restyp)
#get two copies of the sequence. A modeller trick to get things set up
ali.append_model(mdl1, align_codes=modelname)
# Generate molecular topology for mutant
mdl1.clear_topology()
mdl1.generate_topology(ali[-1])
# Transfer all the coordinates you can from the template native structure
# to the mutant (this works even if the order of atoms in the native PDB
# file is not standard):
#here we are generating the model by reading the template coordinates
mdl1.transfer_xyz(ali)
# Build the remaining unknown coordinates
mdl1.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
#yes model2 is the same file as model1. It's a modeller trick.
mdl2 = Model(env, file=modelname)
#required to do a transfer_res_numb
#ali.append_model(mdl2, atom_files=modelname, align_codes=modelname)
#transfers from "model 2" to "model 1"
mdl1.res_num_from(mdl2,ali)
#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).
mdl1.write(file=modelname+restyp+respos+'.tmp')
mdl1.read(file=modelname+restyp+respos+'.tmp')
#set up restraints before computing energy
#we do this a second time because the model has been written out and read in,
#clearing the previously set restraints
make_restraints(mdl1, ali)
#a non-bonded pair has to have at least as many selected atoms
mdl1.env.edat.nonbonded_sel_atoms=1
sched = autosched.loop.make_for_model(mdl1)
#only optimize the selected residue (in first pass, just atoms in selected
#residue, in second pass, include nonbonded neighboring atoms)
#set up the mutate residue selection segment
s = Selection(mdl1.chains[chain].residues[respos])
mdl1.restraints.unpick_all()
mdl1.restraints.pick(s)
s.energy()
s.randomize_xyz(deviation=4.0)
mdl1.env.edat.nonbonded_sel_atoms=2
optimize(s, sched)
#feels environment (energy computed on pairs that have at least one member
#in the selected)
mdl1.env.edat.nonbonded_sel_atoms=1
optimize(s, sched)
s.energy()
#give a proper name
mdl1.write(file=modelname+restyp+respos+'.pdb')
#delete the temporary file
os.remove(modelname+restyp+respos+'.tmp')
Result-
1o9u_bjoin.B99990003GLY232.pdb
5. Loop refinement
# Loop refinement of an existing model
from modeller import *
from modeller.automodel import *
#from modeller import soap_loop
log.verbose()
env = Environ()
# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']
# Create a new class based on 'LoopModel' so that we can redefine
# select_loop_atoms (necessary)
class MyLoop(LoopModel):
# This routine picks the residues to be refined by loop modeling
def select_loop_atoms(self):
# One loop in chain A from residue 19 to 28 inclusive
#return Selection(self.residue_range('19:A', '28:A'))
# Two loops simultaneously
return Selection(self.residue_range('181:A', '183:A'),
self.residue_range('231:A', '233:A'))
m = MyLoop(env,
inimodel='1o9u_bjoin.B99990003GLY232.pdb', # initial model of the target
sequence='1o9u_bjoin', # code of the target
loop_assess_methods=assess.DOPE) # assess loops with DOPE
# loop_assess_methods=soap_loop.Scorer()) # assess with SOAP-Loop
m.loop.starting_model= 20 # index of the first loop model
m.loop.ending_model = 23 # index of the last loop model
m.loop.md_level = refine.very_fast # loop refinement method
m.make()
Result-
1o9u_bjoin.BL00230001.pdb
DOPE Score = -1349.11682
RMSD = 0.00
