dear modellers,
I am tring to mutate a particular residue in a protein using the below script but getting the error, could someone help, Error Traceback (most recent call last): File "(stdin)", line 93, in ? File "C:\Program Files\Modeller8v2\modlib\modeller\model.py", line 220, in generate_topology io=io.modpt, libs=libs.modpt, **vars) File "C:\Program Files\Modeller8v2\modlib\modeller\util\top.py", line 197, in generate_topology return _modeller.generate_topology(mdl, aln, io, libs, *args) _modeller.error: gener___479E> Must not use patching residues here. Residue type index, residue type: 35
------------------------------------------------------------------------------------------------------------------------------------------------------------------ command: mod8v2 - < mutate.py 1qhd 2 GLU A > 1qhd.log ------------------------------------------------------------------------------------------------------------------------------------------------------------------ The script
import sys import os
# # mutate_model.py # # Usage: mod8v2 - < mutate_model.py modelname respos resname chain > logfile # # Example: mod8v2 - < 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 ='2', res_types ='GLU', 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 ='2', res_types ='GLU', 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');
----------------------------------------------------------------------------------------------------------------------------------------------------------
Thanks in advance prabhakar
Prabhakar. G. wrote: > I am tring to mutate a particular residue in a protein using the below > script but getting the error, > could someone help, > Error > Traceback (most recent call last): > File "(stdin)", line 93, in ? > File "C:\Program Files\Modeller8v2\modlib\modeller\model.py", line > 220, in generate_topology > io=io.modpt, libs=libs.modpt, **vars) > File "C:\Program Files\Modeller8v2\modlib\modeller\util\top.py", line > 197, in generate_topology > return _modeller.generate_topology(mdl, aln, io, libs, *args) > _modeller.error: gener___479E> Must not use patching residues here. > Residue type index, residue type: 35 > > ------------------------------------------------------------------------------------------------------------------------------------------------------------------ > > command: mod8v2 - < mutate.py 1qhd 2 GLU A > 1qhd.log
Modeller is confused because the 1qhd PDB file has a HETATM connected to the main chain. Modeller knows about the ACE residue type, but thinks it is the CHARMM ACE type, which is a patching residue, not a residue in its own right, and so it doesn't know what to do with the ACE residue.
The simplest solution would be to either manually remove the ACE residue from your PDB file, or set model_segment so that Modeller ignores residue 0. You can always use model.patch() to apply the ACE patch to the first residue, if you want to preserve the acetylated N terminus.
Ben Webb, Modeller Caretaker