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,

Roland F. Rivera-Santiago
Ph.D. Candidate
Biochemistry & Molecular Biophysics Graduate Group
Speicher Lab, Wistar Institute
E-mail: roland.f.rivera.santiago@gmail.com