Hi there,

I'm having some trouble using the model_patch() command to create a disulfide bond. If I use the following alignment file:

>P1;1PRX_edited structure:1PRX_edited.pdb:2 :B:224 :B: : : : -PGGLLLGDVAPNFEANTTVGRIRFHDFLGDSWGILFSHPRDFT--.-----------PE FAKRNVKLIALSIDS---------------SEEPTEKLPFPIIDDRNRELAILLGMLDPA E-----MPVTARVVFVFGPDKKLKLSILYPATTGRNFDEILRVVISLQLTAEKRVATPVD WKDGDSVMVLPTIPEEEAKKLFPKGVFTKELPSGKKYLRYTPQP*

>P1;1PRX_WT sequence:1PRX_WT:2 :A:224 :A: : : : MPGGLLLGDVAPNFEANTTVGRIRFHDFLGDSWGILFSHPRDFTPVCTTELGRAAKLAPE FAKRNVKLIALSIDSVEDHLAWSKDINAYNCEEPTEKLPFPIIDDRNRELAILLGMLDPA EKDEKGMPVTARVVFVFGPDKKLKLSILYPATTGRNFDEILRVVISLQLTAEKRVATPVD WKDGDSVMVLPTIPEEEAKKLFPKGVFTKELPSGKKYLRYTPQP*

in conjunction with the following script:

# Modeling a cysteine bond in a new model from modeller import * from modeller.automodel import * # Load the automodel class from modeller.scripts import complete_pdb

log.verbose() # Verbose error description - helps with troubleshooting! env = environ() # Defines the environment

# directories for input atom files env.io.atom_files_directory = ['.','../atom-files']

#Enables MODELLER to read HETATM entries from a PDB file env.io.hetatm = True

env.libs.topology.read(file='$(LIB)/top_heav.lib')

env.libs.parameters.read(file='$(LIB)/par.lib')

# Make the disulfide bond: def patches (mdl): mdl.patch(residue_type='DISU', residues=(mdl.residues['47:A'], mdl.residues['91:A']))

# Read the sequence: code = '1PRX' mdl = complete_pdb(env, code, special_patches=patches)

# Create the stereochemical restraints: sel = selection(mdl) mdl.restraints.make(sel, restraint_type='stereo', spline_on_site=False)

# Calculate the energy to test the disulfide: sel.energy()

a = automodel(env, alnfile = '1PRX_edited-1PRX_WT.ali', # alignment filename knowns = '1PRX_edited', # codes of the templates sequence = '1PRX_WT', # code of the target assess_methods=(assess.DOPE, assess.GA341)) a.starting_model= 1 # index of the first model a.ending_model = 1 # index of the last model # (together, they determine how many models to calculate) a.make() # do homology modeling

It generates a PDB file that has the desired sequence, but no disulfide bond. If I use the following alignment:

>P1;1PRX structure:1PRX.pdb:2 :B:224 :B: : : : -PGGLLLGDVAPNFEANTTVGRIRFHDFLGDSWGILFSHPRDFTPV.TTELGRAAKLAPE FAKRNVKLIALSIDSVEDHLAWSKDINAYNSEEPTEKLPFPIIDDRNRELAILLGMLDPA E-----MPVTARVVFVFGPDKKLKLSILYPATTGRNFDEILRVVISLQLTAEKRVATPVD WKDGDSVMVLPTIPEEEAKKLFPKGVFTKELPSGKKYLRYTPQP*

>P1;1PRX_WT sequence:1PRX_WT:2 : :224 : : : : : MPGGLLLGDVAPNFEANTTVGRIRFHDFLGDSWGILFSHPRDFTPVCTTELGRAAKLAPE FAKRNVKLIALSIDSVEDHLAWSKDINAYNCEEPTEKLPFPIIDDRNRELAILLGMLDPA EKDEKGMPVTARVVFVFGPDKKLKLSILYPATTGRNFDEILRVVISLQLTAEKRVATPVD WKDGDSVMVLPTIPEEEAKKLFPKGVFTKELPSGKKYLRYTPQP*

Combined with the following script:

# Modeling a cysteine bond in a new model from modeller import * from modeller.automodel import * # Load the automodel class from modeller.scripts import complete_pdb

log.verbose() # Verbose error description - helps with troubleshooting! env = environ() # Defines the environment

# directories for input atom files env.io.atom_files_directory = ['.','../atom-files']

#Enables MODELLER to read HETATM entries from a PDB file env.io.hetatm = True

env.libs.topology.read(file='$(LIB)/top_heav.lib')

env.libs.parameters.read(file='$(LIB)/par.lib')

# Special distance restraints class MyModel (automodel): # Make the disulfide bond: def patches (mdl): mdl.patch(residue_type='DISU', residues=(mdl.residues['47:'], mdl.residues['91:']))

# Read the sequence: code = '1PRX' mdl = complete_pdb(env, code, special_patches=patches)

# Create the stereochemical restraints: sel = selection(mdl) mdl.restraints.make(sel, restraint_type='stereo', spline_on_site=False)

# Calculate the energy to test the disulfide: sel.energy()

a = MyModel(env, alnfile = '1PRX-1PRX_WT.ali', # alignment filename knowns = '1PRX', # codes of the templates sequence = '1PRX_WT', # code of the target assess_methods=(assess.DOPE, assess.GA341)) a.starting_model= 1 # index of the first model a.ending_model = 1 # index of the last model # (together, they determine how many models to calculate) a.make() # do homology modeling

It returns a KeyError: 'No such residue'. I researched this in the answer archives and found the suggestion to create a .ini file to determine the model residue numbers, which are exactly correspondent to those found in the script. I would appreciate any help deciphering what exactly I'm doing wrong. Thanks,