[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?