Dear all,

I am trying to fix a protein with 6 missing aa. This protein has a resid from 500 to 650, with three HETATM residues (N-glycan) from 800 to 802 ( discontinuous with the protein resnumber). Here is my script to run modeller:

*First step: to build a model with blank chain ID and resnumber from 1 to 150:*

# Homology modeling by the automodel class 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 env.patch_default = False class MyModel(automodel): def select_atoms(self): # Select residues 41 and 46 (PDB numbering) return selection(self.residue_range('41', '46'))

a = MyModel(env, alnfile = 'xy.pir', # alignment filename knowns = 'x', # sequence = 'y', assess_methods=(assess.DOPE, assess.GA341)) # code of the target a.starting_model= 1 # index of the first model a.ending_model = 10 # index of the last model # Very thorough VTFM optimization: a.library_schedule = autosched.slow a.max_var_iterations = 300 # Thorough MD optimization: a.md_level = refine.slow # Repeat the whole cycle 2 times and do not stop unless obj.func. > 1E6 a.repeat_optimization = 3 a.max_molpdf = 1e6

a.make() # Get a list of all successfully built models from a.outputs ok_models = filter(lambda x: x['failure'] is None, a.outputs)

# Rank the models by DOPE score key = 'DOPE score' ok_models.sort(lambda a,b: cmp(a[key], b[key])) # Get top model m = ok_models[0] print "Top model: %s (DOPE score %.3f)" % (m['name'], m[key]) By running this script I did got a model with missing residues fixed and HETATM keeping, but the chain ID was blank (because only one chain was built), and the resnumber was from 1 to 154. I did hope to keep the original chain ID and the resnum, so I try the rename_segments, like:

*Second step: reassign the chain ID and resnumber:*

# Example for: model.rename_segments() # This will assign new PDB single-character chain id's to all the chains # in the input PDB file (here there are two 'chains': protein and the HETATM # water molecules). from modeller import * # Read the MODEL with all HETATM and water records (so there are two 'chains'): env = environ() env.io.atom_files_directory = ['./', ':../atom_files'] env.io.hetatm = True mdl = model(env, file='y') # Assign new segment names and write out the new model: mdl.rename_segments(segment_ids='A', renumber_residues=[500, 650]) mdl.write(file='y.ini') Now in the final model y.ini the chain ID was defined as "A", and resnum was from 500 to 650. However, the HETATM residues was from 651 to 653, but not 800 to 802. So should I modify my scripts to fix the protein and keeping its original chain ID, resnumber and HETATM residues? Thank you very much for your attention.

P.S. The HETATM residues are also in chain A, the same chain with the protein.