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