Hello, 

I have used modeller to attempt to model a loop region in a tetrameric structure. 

To speed up computational time I've started with just the dimer ( which I think also makes sense biologically ) and ran the following script with an .ali file I created. The region of interest is a disordered loop AA~399-417. Following this I used the evaluate model script to pick the best model, which I then used as the template input for the loop.py script. I run the loop.py script at many iterations and then when it finishes I take the top 20 models (out of 1000) by DOPE score and attempt to analyze them. Any advice on how I could improve my method would be great.. I really am at the edge of my understanding here. 

I've edited the name of the protein. I can include it as soon as I can verify that this can stay anonymous. 
Thanks!!!

Scripts +files :

First )using a template pdb and .ali file  I made myself of the Dimer (BIO_headers in chimera and deletion of one dimer of the tetramer) 
Input : 


>P1;ABCD
structureX:ABCD.pdb:   3 :A:+990:B:::-1.00:-1.00
--TSWSDRLQNAADMPANMDKHALKKYRREAYHRVFVNRSLAMEKIKCFGFNMDYTLAVYKSPEYESLGFELTVE
RLVSIGYPQELLSFAYDSTFPTRGLVFDTLYGNLLKVDAYGNLLVCAHGFNFIRGPETREQYPNKFIQRDDTERF
YILNTLFNLPETYLLACLVDFFTNCPRYTSCETGFKDGDLFMSYRSMFQDVRDAVDWVHYKGSLKEKTVENLEKY
VVKDGKLPLLLSRMKEVGKVFLATNSDYKYTDKIMTYLFDFPHGPKPGSSHRPWQSYFDLILVDARKPLFFGEGT
VLRQVDTKTGKLKIGTYTGPLQHGIVYSGGSSDTICDLLGAKGKDILYIGDHIFGDILKSKKRQGWRTFLVIPEL
AQELHVWTDKSSLFEELQSLDIFLAS----------------SIQRRIKKVTHDMDMCYGMMGSLFRSGSRQTLF
ASQVMRYADLYAASFINLLYYPFSYLFRAAHVLMPHES/--TSWSDRLQNAADMPANMDKHALKKYRREAYHRVFVNRSLAMEKIKCFGFNMDYTLAVYKSPEYESLGFELTVE
RLVSIGYPQELLSFAYDSTFPTRGLVFDTLYGNLLKVDAYGNLLVCAHGFNFIRGPETREQYPNKFIQRDDTERF
YILNTLFNLPETYLLACLVDFFTNCPRYTSCETGFKDGDLFMSYRSMFQDVRDAVDWVHYKGSLKEKTVENLEKY
VVKDGKLPLLLSRMKEVGKVFLATNSDYKYTDKIMTYLFDFPHGPKPGSSHRPWQSYFDLILVDARKPLFFGEGT
VLRQVDTKTGKLKIGTYTGPLQHGIVYSGGSSDTICDLLGAKGKDILYIGDHIFGDILKSKKRQGWRTFLVIPEL
AQELHVWTDKSSLFEELQSLDIFLAS----------------SIQRRIKKVTHDMDMCYGMMGSLFRSGSRQTLF
ASQVMRYADLYAASFINLLYYPFSYLFRAAHVLMPHES*

>P1;X
sequence:ABCD:     : :     : ::: 0.00: 0.00
MSTSWSDRLQNAADMPANMDKHALKKYRREAYHRVFVNRSLAMEKIKCFGFDMDYTLAVYKSPEYESLGFELTVE
RLVSIGYPQELLSFAYDSTFPTRGLVFDTLYGNLLKVDAYGNLLVCAHGFNFIRGPETREQYPNKFIQRDDTERF
YILNTLFNLPETYLLACLVDFFTNCPRYTSCETGFKDGDLFMSYRSMFQDVRDAVDWVHYKGSLKEKTVENLEKY
VVKDGKLPLLLSRMKEVGKVFLATNSDYKYTDKIMTYLFDFPHGPKPGSSHRPWQSYFDLILVDARKPLFFGEGT
VLRQVDTKTGKLKIGTYTGPLQHGIVYSGGSSDTICDLLGAKGKDILYIGDHIFGDILKSKKRQGWRTFLVIPEF
AQELHVWTDKSSLFEELQSLDIFLAELYKHLDSSSNERPDISSIQRRIKKVTHDMDMCYGMMGSLFRSGSRQTLF
ASQVMRYADLYAASFINLLYYPFSYLFRAAHVLMPHES/MSTSWSDRLQNAADMPANMDKHALKKYRREAYHRVFVNRSLAMEKIKCFGFDMDYTLAVYKSPEYESLGFELTVE
RLVSIGYPQELLSFAYDSTFPTRGLVFDTLYGNLLKVDAYGNLLVCAHGFNFIRGPETREQYPNKFIQRDDTERF
YILNTLFNLPETYLLACLVDFFTNCPRYTSCETGFKDGDLFMSYRSMFQDVRDAVDWVHYKGSLKEKTVENLEKY
VVKDGKLPLLLSRMKEVGKVFLATNSDYKYTDKIMTYLFDFPHGPKPGSSHRPWQSYFDLILVDARKPLFFGEGT
VLRQVDTKTGKLKIGTYTGPLQHGIVYSGGSSDTICDLLGAKGKDILYIGDHIFGDILKSKKRQGWRTFLVIPEF
AQELHVWTDKSSLFEELQSLDIFLAELYKHLDSSSNERPDISSIQRRIKKVTHDMDMCYGMMGSLFRSGSRQTLF
ASQVMRYADLYAASFINLLYYPFSYLFRAAHVLMPHES*


# Homology modeling by the automodel class
#
# Demonstrates how to build multi-chain models, and symmetry restraints
#
from modeller import *
from modeller.automodel import * # Load the automodel class
 
log.verbose()
 
# Override the 'special_restraints' and 'user_after_single_model' methods:
class MyModel(automodel):
        def special_restraints(self, aln):
          # Constrain the A, B, C and D chains to be identical
          s1 = selection(self.chains['A']).only_atom_types('CA')
          s2 = selection(self.chains['B']).only_atom_types('CA')
          self.restraints.symmetry.append(symmetry(s1, s2, 1))
                   
        def user_after_single_model(self):
          # Report on symmetry violations greater than 1A after building
          # each model:
          self.restraints.symmetry.report(1)
         
env = environ()
# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']
 
# Be sure to use 'MyModel' rather than 'automodel' here!
a = MyModel(env,
            alnfile = 'WT_dimer_NoCterm.ali' ,       # alignment filename
            knowns = 'ABCD',             # codes of the templates
            sequence = 'ABCD')           # code of the target
           
a.starting_model= 1                      # index of the first model
a.ending_model = 20                      # index of the last model
                                        # (determines how many models to calculate)
a.make()                                 # do homology modeling
                                         #
                                         # class MyModel(automodel):

2nd) Run eval model and pick top model to use as input for the loop.py script below --->

# Loop refinement of an existing model
from modeller import *
from modeller.automodel import *

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = './:../atom_files'

# Create a new class based on 'loopmodel' so that we can redefine
# select_loop_atoms (necessary)
class MyLoop(loopmodel):
    # This routine picks the residues to be refined by loop modeling
    def select_loop_atoms(self):
        # 10 residue insertion 
        return selection(self.residue_range('399:A', '417:A'))
        
m = MyLoop(env,
           inimodel='ABCD_Top_Model.pdb', # initial model of the target
           sequence='ABCD')          # code of the target

m.loop.starting_model= 1          # index of the first loop model 
m.loop.ending_model  = 1000       # index of the last loop model
m.loop.md_level = refine.slow # loop refinement method; this yields
                                   # models quickly but of low quality;
                                   # use refine.slow for better models

m.make()


3rd) Run model_energies script

Then I rank them using this perl one liner ( just cats and displays the dope and the file ) can run in any unix terminal. 

cat model*.log |perl -lane' print if /DOPE\ score/ || /ABCD.BL/'| perl -pe 's/ +/\t/g'|perl -pe 's/:\t//g' |cut -f3,4|perl -pe 's/\n/\t/'|perl -pe 's/(\d)\t(ABCD)/$1\n$2/g'|sort -k2,2n|head -20


4) Take top 20 models and analyze 

Any help with where I'm doing things wrong or how I can improve would be awesome.  How many models in loop.py would you recommend to get a good answer? Is there away to interpret the models created as a whole ....for example say something like "60% of the time  residue A is next to residue B and 20% of the time its over in this pocket , therefor".... etc 

My major aim is to compare different ligand bound crystal states, or mutation states, and show that the models created are different between them. Then use this to design wet lab experiments to test. 


Thanks again