[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[modeller_usage] Modeling multiple chains with two ligands



Hello,
I am having a great problem with modeller 9v1.
My aim is to model 2 chain protein from a 2 chain template with a ligands. I want two ligands: Ca metal ion and SOX (Penicillin G sulphoxide). My template pdb file briefly looks like this:
----------------------------------------------------------------------------------------------------
ATOM 1 N SER A 3 27.496 31.222 28.262 1.00 48.33 N ATOM 2 CA SER A 3 27.785 31.042 26.803 1.00 43.97 C ... etc to the end of chain a
TER    1679      ALA A 209
ATOM 1684 N SER B 1 14.845 1.596 -1.736 1.00 6.93 N ATOM 1685 CA SER B 1 15.548 2.917 -1.659 1.00 6.75 C ATOM 1686 C SER B 1 16.489 3.059 -2.844 1.00 8.04 C
.....
ATOM 6113 NE ARG B 557 26.460 16.023 30.916 1.00 23.74 N ATOM 6114 CZ ARG B 557 27.175 15.276 30.094 1.00 19.20 C ATOM 6115 NH1 ARG B 557 28.016 15.921 29.310 1.00 21.83 N ATOM 6116 NH2 ARG B 557 27.069 13.952 30.055 1.00 20.49 N ATOM 6117 OXT ARG B 557 23.973 18.664 34.656 1.00 29.98 O
TER    6118      ARG B 557
.... Cutting useless HETATM
HETATM 6159 CA CA B1568 23.187 -16.230 6.242 1.00 10.86 CA << Ca ion HETATM 6160 O8 SOX B1569 13.279 0.047 0.037 1.00 37.00 O << SOX ... HETATM 6161 C7 SOX B1569 13.688 -1.062 0.060 1.00 38.52 C
..
HETATM 6166 O12 SOX B1569 13.790 -4.072 -2.049 1.00 54.27 O ...etc to the end of SOX (I guess resid for SOX = 569, chain B1, ?)
-----------------------------------------------------------------------------------------------------

I have made alignment in PIR using chain breaks "/" and block residues "." for Ca and SOX.
My alignment looks like this:
-----------------------------------------------------------------------------------------------------
>P1;1gm9
structureX:1gm9.pdb:3:A:569:B1:::: #Starting from A resid =3(because first two are missing) and finishing with SOX
--SSSEIKIVRDEYG-PHIYANDTWHLFYGYGYVVAQDRLFQMEMARRSTQGTVAEVLGK
DFVKFDKDIRRNYWPDAIRAQIAALSPEDMSILQGYADGMNAWIDKVNTNPETLLPKQFN
TFGFTPKRWEPFDVAMIFVGTMANRFSDSTSEIDNLALLTALKDKYGVSQGMAVFNQLKW
LVNPSAPTTIAVQESNYPLKFNQQNSQTA/  #Chain break = ending chain A
SNMWVIGKSKAQDAKAIMVNGPQFGWYAPAY
TYGIGLHGAGYDVTGNTPFAYPGLVFGHNGVISWGSTAGFGDDVDIFAERLSAEKPGYYL
HNGKWVKMLSREETITVKNGQAETFTVWRTVHGNILQTDQTTQTAYAKSRAWDGKEVASL
LAWTHQMKAKNWQEWTQQAAKQALTINWYYADVNGNIGYVHTGAYPDRQSGHDPRLPVPG
TGKWDWKGLLPFEMNPKVYNPQSGYIANWNNSPQKDYPASDLFAFLWGGADRVTEIDRLL
EQKPRLTADQAWDVIRQTSRQDLNLRLFLPTLQAATSGLTQSDPRRQLVETLTRWDGINL
LNDDGKTWQQPGSAILNVWLTSMLKRTVVAAVPMPFDKWYSASGYETTQDGPTGSLNISV
GAKILYEAVQGDKSPIPQAVDLFAGKPQQEVVLAALEDTWETLSKRYGNNVSNWKTPAMA
LTFRANNFFGVPQAAAEETRHQAEYQNRGTENDMIVFSPTTSDRPVLAWDVVAPGQSGFI
APDGTVDKHYEDQLKMYENFGRKSLWLTKQDVEAHKESQEVLHVQR/..*
#Chain break = ending chain B and ".." for two HETATM

>P1;1klc
sequence:1klc:     : :     : ::: 0.00: 0.00
ASPPTEVKIVRDEYGMPHIYADDTYRLFYGYGYVVAQDRLFQMEMARRSTQGTVSEVLGK
AFVSFDKDIRQNYWPDSIRAQIASLSAEDKSILQGYADGMNAWIDKVNASPDKLLPQQFS
TFGFKPKHWEPFDVAMIFVGTMANRFSDSTSEIDNLALLTAVKDKYGNDEGMAVFNQLKW
LVNPSAPTTIAARESSYPLKFDLQNTQTA/
SNMWVIGKNKAQDAKAIMVNGPQFGWYAPAY
TYGIGLHGAGYDVTGNTPFAYPGLVFGHNGTISWGSTAGFGDDVDIFAEKLSAEKPGYYQ
HNGEWVKMLSRKETIAVKDGQPETFTVWRTLDGNVIKTDTRTQTAYAKARAWAGKEVASL
LAWTHQMKAKNWPEWTQQAAKQALTINWYYADVNGNIGYVHTGAYPDRQPGHDPRLPVP-
DGKWDWKGLLSFDLNPKVYNPQSGYIANWNNSPQKDYPASDLFAFLWGGADRVTEIDTIL
DKQPRFTADQAWDVIRQTSLRDL-LRLFLPALKDATANLAENDPRRQLVDKLASWDGENL
VNDDGKTYQQPGSAILNAWLTSMLKRTVVAAVPAPFGKWYSASGYETTQDGPTGSLNISV
GAKILYEALQGDKSPIPQAVDLFGGKPEQEVILAALDDAWQTLSKRYGNDVTGWKTPAMA
LTFRANNFFGVPQAAAKEARHQAEYQNRGTENDMIVFSPTSGNRPVLAWDVVAPGQSGFI
APDGKADKHYDDQLKMYESFGRKSLWLTPQDVDEHKESQEVLQVQR/--*
---------------------------------------------------------------------------------------------------------


Then I used advanced example as a template to write my input script:
---------------------------------------------------------------------------------------------------------
from modeller import *
from modeller.automodel import *

class mymodel(automodel):
   def special_restraints(self, aln):
       rsr = self.restraints
for ids in (('NH2:471:B', 'O12:766:C'), # That is for SOX, no restraints for Ca
                       ('ND2:449:B', 'O16:766:C'),
                       ('CE2:280:B', 'C5:766:C')):
           atoms = [self.atoms[i] for i in ids]
           rsr.add(forms.upper_bound(group=physical.upper_distance,
                                     feature=features.distance(*atoms),
                                     mean=3.04, stdev=0.1))
env = environ()
env.io.hetatm = True
a = mymodel(env, alnfile='1klc_1gm9.ali',
             knowns='1gm9m, 1klc_model', sequence='1klc')
a.starting_model = 1
a.ending_model = 1
a.make()
------------------------------------------------------------------------------------------------------

Modeller starts with this error: "import site' failed; use -v for traceback", (I guess it is not a big problem because in other cases it does not cause any trouble) , then program starts calculating, and then I always get :
-----------------------------------------------------------------------------------------------------
Traceback (most recent call last):
 File "model-multiple-hetero.py", line 47, in ?
   a.make()
File "C:\Program Files\Modeller9v1\modlib\modeller\automodel\automodel.py", li
ne 108, in make
   self.homcsr(exit_stage)
File "C:\Program Files\Modeller9v1\modlib\modeller\automodel\automodel.py", li
ne 427, in homcsr
   self.mkhomcsr(selection(self), aln)
File "C:\Program Files\Modeller9v1\modlib\modeller\automodel\automodel.py", li
ne 537, in mkhomcsr
   self.special_restraints(aln)
 File "model-multiple-hetero.py", line 10, in special_restraints
   atoms = [self.atoms[i] for i in ids]
File "C:\Program Files\Modeller9v1\modlib\modeller\coordinates.py", line 127,
in __getitem__
   (self.offset, self.length, self.suffix))
File "C:\Program Files\Modeller9v1\modlib\modeller\util\modutil.py", line 76,
in handle_seq_indx
   int_indx = lookup_func(*args)
File "C:\Program Files\Modeller9v1\modlib\modeller\coordinates.py", line 42, i
n _indxatm
   raise KeyError, ("No such atom: %s" % indx)
KeyError: 'No such atom: O12:766:C'
----------------------------------------------------------------------------------------------------------
That means, that it can not see atom O12:766:C.
I am sure, that my seq is 764 residues long (I mean residues in seq, not position in alignment). So Ca must be 765 and SOX - 766. I have tried different counts (766,765,767...) and different chains (B1, C) but it did not help. I have tried to model only chain B with SOX, but it did not help. If I remove user restraints and just use automodel, it output good model, but without ligands at all!
And the latter is what I absolutely do no understand.

My log has this note about SOX:
-----------------------------------------------------------------------------------------------------------
read_pd_459W> Residue type  SOX not recognized. 'automodel' model building
             will treat this residue as a rigid body.
To use real parameters, add the residue type to ${LIB}/restyp.lib,
             its topology to ${LIB}/top_*.lib, and suitable forcefield
             parameters to ${LIB}/par.lib.
-----------------------------------------------------------------------------------------------------------
And nothing about CA (may be the program takes CA as Calpha?).


So my questions are:
1. What am I doing wrong?
2. Why modeller outputs model with no ligands if I use default automodel with no restraints (env.io.hetatm =True)?.
What do you think? Thank you for your suggestions.

P.S. I am using windows compilation. May it be some platform specific modeller bug?