I am trying to perform loop optimization on a protein structure which contains a non-standard ligand. The optimization I am attempting requires constrains between the ligand and the loops. Protein is chain A, ligand is a single residue defined as chain B, described as heteroatoms.

I had no trouble performing homology modeling from a structure containing the ligand, and it was correctly carried upon once I define env.io.hetatm = True .However this seems not to work anymore when doing loop optimization only on the same structure. When I try to run the job, I always get:

Traceback (most recent call last):
  File "tas2r463pbl_mod.py", line 42, in ?
  File "/usr/lib/modeller9.10/modlib/modeller/automodel/loopmodel.py", line 36, in make
    self.build_seq(self.inimodel, 1)
  File "/usr/lib/modeller9.10/modlib/modeller/automodel/loopmodel.py", line 167, in build_seq
    self.loop_restraints(atmsel, aln)
  File "/usr/lib/modeller9.10/modlib/modeller/automodel/loopmodel.py", line 361, in loop_restraints
  File "tas2r463pbl_mod.py", line 24, in special_restraints
    rsr.add(forms.gaussian(group=physical.bond,feature=features.distance(at['CA:70:A'],at['C16:309:B']), mean=1.00, stdev=0.10))
  File "/usr/lib/modeller9.10/modlib/modeller/coordinates.py", line 376, in __getitem__
    (self.offset, self.length, self.suffix))
  File "/usr/lib/modeller9.10/modlib/modeller/util/modutil.py", line 23, in handle_seq_indx
    int_indx = lookup_func(*args)
  File "/usr/lib/modeller9.10/modlib/modeller/coordinates.py", line 289, in _indxatm
    raise KeyError("No such atom: %s" % indx)
KeyError: 'No such atom: C16:309:B'

Looking in the log, I think this is the relevant warning:

transfe_404W> At least one template is aligned with model residue   309:B
              but no coordinates could be transferred. This usually
              occurs when your input files do not use the official
              PDBv3 atom names. Please check your templates.

The script follows:


#!/usr/bin/env python

from modeller import *
from modeller.automodel import *
from modeller.scripts import complete_pdb

env = environ()

env.io.atom_files_directory = ['.']
env.io.hetatm = True

class MyModel(loopmodel):

    def select_loop_atoms(self):
        #return selection(self.residue_range('68:A','83:A'),self.residue_range('149:A','176:A'),self.residue_range('249:A','258:A'),self.residue_range('309:B','309:B'))
        return selection(self.residue_range('68:A','83:A'),self.residue_range('149:A','176:A'),self.residue_range('249:A','258:A'))

    def special_restraints(self,aln):
        rsr = self.restraints
        at = self.atoms
        rsr.add(forms.gaussian(group=physical.bond,feature=features.distance(at['CA:70:A'],at['C16:309:B']), mean=1.00, stdev=0.10))
        rsr.add(forms.gaussian(group=physical.bond,feature=features.distance(at['CA:71:A'],at['C16:309:B']), mean=1.00, stdev=0.10))
        rsr.add(forms.gaussian(group=physical.bond,feature=features.distance(at['CA:253:A'],at['C16:309:B']), mean=1.00, stdev=0.10))

a = MyModel(env, inimodel="cluster_A.pdb", sequence='tas2r46', loop_assess_methods=(assess.DOPE, assess.GA341, assess.normalized_dope))
#a = MyModel(env, inimodel="cluster_A.pdb", loop_assess_methods=(assess.DOPE, assess.GA341, assess.normalized_dope))
a.loop.starting_model = 1
a.loop.ending_model = 1
a.loop.md_level = refine.very_fast

env = environ()

#for some reason, this has to be done
env = environ()

#in case we output bad models...
ok_models = filter(lambda x: x['failure'] is None, a.outputs)

for model in ok_models:
    mdl = complete_pdb(env, model['name'])
    s = selection(mdl)
    s.assess_dope(output='ENERGY_PROFILE NO_REPORT', file= model['name']+'.profile', normalize_profile=True, smoothing_window=15)

Any help on what information should I look for to troubleshoot this is very welcome.