Modeller logo
  • Diff for "Mutate model"
Differences between revisions 1 and 2
Revision 1 as of 2005-09-28 18:58:00
Size: 171
Comment:
Revision 2 as of 2005-09-28 19:01:08
Size: 5625
Comment: Add the script itself.
Deletions are marked like this. Additions are marked like this.
Line 2: Line 2:


{{{#!python
import sys
import os

#
# mutate_model.py
#
# Usage: modCVS - < mutate_model.py modelname respos resname chain > logfile
#
# Example: modCVS - < 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.
#


def optimize(mdl1):
    #conjugate gradient
    for step in range(0,len(mdl1.schedule)):
        mdl1.schedule.step = step + 1
        mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)
    #md
    refine(mdl1)
    mdl1.optimize(optimization_method=1,max_iterations=200,min_atom_shift=0.001)



#molecular dynamics
def refine(mdl1):
    # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
    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:
            mdl1.optimize(cap_atom_shift=0.39, md_time_step=4.0, optimization_method=3,
                          md_return='FINAL', 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()
   for typ in ('stereo', 'phi-psi_binormal'):
       rsr.make(restraint_type=typ, aln=aln, spline_on_site=True)
   for typ in ('omega', 'chi1', 'chi2', 'chi3', 'chi4'):
       rsr.make(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()

env = environ(rand_seed=-49837)
env.io.hetatm = True
#soft sphere potential
env.edat.dynamic_sphere=False
#lennard-jones potential (more accuate)
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
mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',selection_segment = (respos+':'+chain, respos+':'+chain), atom_types = 'ALL', selection_status = 'INITIALIZE')

#perform the mutate residue operation
mdl1.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.generate_topology(ali, sequence=modelname)


# 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

mdl1.schedule.make(library_schedule=6)

#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
mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',selection_segment = (respos+':'+chain, respos+':'+chain), atom_types = 'ALL', selection_status = 'INITIALIZE')

mdl1.restraints.unpick_all()
mdl1.restraints.pick()

mdl1.energy()

mdl1.randomize_xyz(deviation=4.0)

mdl1.env.edat.nonbonded_sel_atoms=2
optimize(mdl1)

#feels environment (energy computed on pairs that have at least one member in the selected)
mdl1.env.edat.nonbonded_sel_atoms=1
optimize(mdl1)

mdl1.energy()

#give a proper name
mdl1.write(file=modelname+restyp+respos+'.pdb')

#delete the temporary file
os.remove(modelname+restyp+respos+'.tmp');

}}}

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:   modCVS - < mutate_model.py modelname respos resname chain > logfile
   8 #    
   9 #     Example: modCVS - < 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',selection_segment = (respos+':'+chain, respos+':'+chain), atom_types = 'ALL', selection_status = 'INITIALIZE')
  83 
  84 #perform the mutate residue operation
  85 mdl1.mutate(residue_type=restyp)
  86 
  87 #get two copies of the sequence.  A modeller trick to get things set up
  88 ali.append_model(mdl1, align_codes=modelname)
  89 
  90 # Generate molecular topology for mutant
  91 mdl1.generate_topology(ali, sequence=modelname)
  92 
  93 
  94 # Transfer all the coordinates you can from the template native structure
  95 # to the mutant (this works even if the order of atoms in the native PDB 
  96 # file is not standard):
  97 #here we are generating the model by reading the template coordinates
  98 mdl1.transfer_xyz(ali)
  99 
 100 # Build the remaining unknown coordinates
 101 mdl1.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
 102 
 103 #yes model2 is the same file as model1.  It's a modeller trick.
 104 mdl2 = model(env, file=modelname)
 105 
 106 #required to do a transfer_res_numb
 107 #ali.append_model(mdl2, atom_files=modelname, align_codes=modelname)
 108 #transfers from "model 2" to "model 1"
 109 mdl1.res_num_from(mdl2,ali)
 110 
 111 #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).
 112 
 113 mdl1.write(file=modelname+restyp+respos+'.tmp')
 114 mdl1.read(file=modelname+restyp+respos+'.tmp')
 115 
 116 #set up restraints before computing energy
 117 #we do this a second time because the model has been written out and read in, clearing the previously set restraints
 118 make_restraints(mdl1, ali)
 119 
 120 #a non-bonded pair has to have at least as many selected atoms
 121 mdl1.env.edat.nonbonded_sel_atoms=1
 122 
 123 mdl1.schedule.make(library_schedule=6)
 124 
 125 #only optimize the selected residue (in first pass, just atoms in selected residue, in second pass, include nonbonded neighboring atoms)
 126 #set up the mutate residue selection segment
 127 mdl1.pick_atoms(pick_atoms_set=1, selection_from ='ALL', res_types ='ALL',selection_segment = (respos+':'+chain, respos+':'+chain), atom_types = 'ALL', selection_status = 'INITIALIZE')
 128 
 129 mdl1.restraints.unpick_all()
 130 mdl1.restraints.pick()
 131 
 132 mdl1.energy()
 133 
 134 mdl1.randomize_xyz(deviation=4.0)
 135 
 136 mdl1.env.edat.nonbonded_sel_atoms=2
 137 optimize(mdl1)
 138 
 139 #feels environment (energy computed on pairs that have at least one member in the selected)
 140 mdl1.env.edat.nonbonded_sel_atoms=1
 141 optimize(mdl1)
 142 
 143 mdl1.energy()
 144 
 145 #give a proper name
 146 mdl1.write(file=modelname+restyp+respos+'.pdb')
 147 
 148 #delete the temporary file
 149 os.remove(modelname+restyp+respos+'.tmp');