If your template contains a ligand (or other HETATM residue) then MODELLER can transfer this into your generated model. This is done first by setting env.io.hetatm to True, which instructs MODELLER to read HETATM records from your template PDB files, and then by using the BLK ('.') residue type in your alignment (both in the template(s) and the model sequence) to copy the ligand(s) as a rigid body into the model.
# Homology modeling with ligand transfer from the template
from modeller import * # Load standard Modeller classes
from modeller.automodel import * # Load the automodel class
log.verbose() # request verbose output
env = environ() # create a new MODELLER environment to build this model in
# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']
# Read in HETATM records from template PDBs
env.io.hetatm = True
a = automodel(env,
alnfile = 'align-ligand.ali', # alignment filename
knowns = '5fd1', # codes of the templates
sequence = '1fdx') # code of the target
a.starting_model= 4 # index of the first model
a.ending_model = 4 # index of the last model
# (determines how many models to calculate)
a.make() # do the actual homology modeling
C; Similar to alignment.ali, but with ligands included >P1;5fd1 structureX:5fd1:1 :A:108 :A:ferredoxin:Azotobacter vinelandii: 1.90: 0.19 AFVVTDNCIKCKYTDCVEVCPVDCFYEGPNFLVIHPDECIDCALCEPECPAQAIFSEDEVPEDMQEFIQLNAELA EVWPNITEKKDPLPDAEDWDGVKGKLQHLER..* >P1;1fdx sequence:1fdx:1 : :56 : :ferredoxin:Peptococcus aerogenes: 2.00:-1.00 AYVINDSC--IACGACKPECPVNIIQGS--IYAIDADSCIDCGSCASVCPVGAPNPED----------------- -------------------------------..*
Note that by turning on env.io.hetatm, all HETATM records are read from your templates, so all of these must be listed in your alignment. Use a single '.' character for each HETATM residue in the template sequence in your alignment.2.1 MODELLER always reads PDB residues in the order they're written in the PDB file, so if you have a ligand at the end of PDB file, put the '.' residue at the end of the sequence in the alignment too. You will often see a chain break ('/') immediately preceding '.' residues in example alignments. That's only necessary if you want to force the ligands to have a different chain ID to the amino acids. (If you want them in the same chain, leave out the chain break.)
To get the ligand into your model, you must align a residue in the model with the desired HETATM residue in the template. Use a single '.' residue in your model sequence in your alignment for each ligand you want in the model. This must be aligned with a suitable ligand in the template sequence. If you have extra HETATM ligands in the template which you don't want in the model, simply align them with a gap ('-') in the model sequence. If you have multiple templates, you can copy ligands from any suitable template -- just align the '.' residue in the model with the desired template sequence ligand.
automodel builds restraints on these ligands to keep their geometry and environment reasonably similar to the template, by restraining some intra-HETATM, inter-HETATM, and HETATM-protein distances to their template values. See automodel.nonstd_restraints() for more information.
You can also treat ligands flexibly by defining topology and parameter information. See section 5.2.1 for more information, and the example in the advanced modeling tutorial, at http://salilab.org/modeller/tutorial/advanced.html.
If you want to add ligands to your model which are not present in your template, you will need to do some docking studies, which are beyond the scope of the MODELLER program.
To read in water residues, set env.io.water to True and use the 'w' residue type in your alignment.
To read in hydrogen atoms, set env.io.hydrogen to True. This is not generally necessary, as if you want to build an all hydrogen model, it is easiest just to use the allhmodel class, which turns this on for you automatically; see section 2.2.5.