In [9]:
import sys 
sys.path.append('/usr/lib/modeller9.12/modlib/')
sys.path.append('/usr/lib/modeller9.12/lib/x86_64-intel8/python2.5/')
import modeller 
import _modeller
import modeller.automodel 
#impo modeller.automodel import *
In [10]:
env=modeller.environ()
env.io.hetatm=True
                         MODELLER 9.12, 2013/06/12, r9480

     PROTEIN STRUCTURE MODELLING BY SATISFACTION OF SPATIAL RESTRAINTS


                     Copyright(c) 1989-2013 Andrej Sali
                            All Rights Reserved

                             Written by A. Sali
                               with help from
           B. Webb, M.S. Madhusudhan, M-Y. Shen, M.A. Marti-Renom,
                N. Eswar, F. Alber, M. Topf, B. Oliva, A. Fiser,
                    R. Sanchez, B. Yerkovich, A. Badretdinov,
                      F. Melo, J.P. Overington, E. Feyfant
                 University of California, San Francisco, USA
                    Rockefeller University, New York, USA
                      Harvard University, Cambridge, USA
                   Imperial Cancer Research Fund, London, UK
              Birkbeck College, University of London, London, UK


Kind, OS, HostName, Kernel, Processor: 4, Linux srv 3.13.0-88-generic x86_64
Date and time of compilation         : 2013/06/11 21:33:06
MODELLER executable type             : x86_64-intel8
Job starting time (YY/MM/DD HH:MM:SS): 2016/11/29 19:02:54

In [1]:
cd /tmp/mod-teach/
/tmp/mod-teach

Download template

In [5]:
%%bash
wget http://www.pdb.org/pdb/files/1lmp.pdb -O 1lmp.pdb
--2016-11-29 18:36:23--  http://www.pdb.org/pdb/files/1lmp.pdb
Resolving www.pdb.org (www.pdb.org)... 132.249.231.77
Connecting to www.pdb.org (www.pdb.org)|132.249.231.77|:80... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: http://www.rcsb.org/pdb/files/1lmp.pdb [following]
--2016-11-29 18:36:23--  http://www.rcsb.org/pdb/files/1lmp.pdb
Resolving www.rcsb.org (www.rcsb.org)... 132.249.231.10
Connecting to www.rcsb.org (www.rcsb.org)|132.249.231.10|:80... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: http://files.rcsb.org/view/1lmp.pdb [following]
--2016-11-29 18:36:23--  http://files.rcsb.org/view/1lmp.pdb
Resolving files.rcsb.org (files.rcsb.org)... 132.249.213.140
Connecting to files.rcsb.org (files.rcsb.org)|132.249.213.140|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: unspecified [text/plain]
Saving to: ‘1lmp.pdb’

     0K .......... .......... .......... .......... ..........  127K
    50K .......... .......... .......... .......... ..........  254K
   100K .......... .......... ........                          279M=0.6s

2016-11-29 18:36:25 (217 KB/s) - ‘1lmp.pdb’ saved [131301]

Download SEQ from uniprot

In [7]:
%%bash
wget http://www.uniprot.org/uniprot/P61626.fasta
--2016-11-29 18:40:00--  http://www.uniprot.org/uniprot/P61626.fasta
Resolving www.uniprot.org (www.uniprot.org)... 128.175.240.211, 193.62.192.81, 193.62.193.81
Connecting to www.uniprot.org (www.uniprot.org)|128.175.240.211|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 217 [text/plain]
Saving to: ‘P61626.fasta’

     0K                                                       100% 42.5M=0s

2016-11-29 18:40:00 (42.5 MB/s) - ‘P61626.fasta’ saved [217/217]

In [8]:
ls
1lmp.pdb  P61626.fasta
In [11]:
alignm=modeller.alignment(env)
alignm.append(file='P61626.fasta', align_codes='all',alignment_format='FASTA')
mdl = modeller.model(env, file='1lmp.pdb', model_segment=('FIRST:'+'A', 'LAST:'+'A'))
alignm.append_model(mdl, atom_files='1lmp.pdb', align_codes='1lmp_A')

alignm[0].code='LYSC_HUMAN'
read_pd_459W> Residue type  NAG 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.
rdpdb___459W> Residue type  NDG 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.
In [14]:
alignm.salign()
s= alignm[0]
for i in range(len(alignm)):
    print "%s identical to wt in %4.1f perc" %(alignm[i].code,s.get_sequence_identity(alignm[i]))
alignm.write(file='all_in_one.ali', alignment_format='PIR')
SALIGN_____> adding the next group to the alignment; iteration    1
LYSC_HUMAN identical to wt in 100.0 perc
1lmp_A identical to wt in 68.2 perc
In [12]:
s=alignm[0]
print s.code
a = modeller.automodel.automodel(env, alnfile='all_in_one.ali', knowns= ('1lmp_A'), sequence = s.code )
a.name='mod'+s.code
a.starting_model = 1
a.ending_model = 1
a.make()
    
LYSC_HUMAN
fndatmi_285W> Only      129 residues out of      132 contain atoms of type  CA
              (This is usually caused by non-standard residues, such
              as ligands, or by PDB files with missing atoms.)
fndatmi_285W> Only      129 residues out of      132 contain atoms of type  CA
              (This is usually caused by non-standard residues, such
              as ligands, or by PDB files with missing atoms.)

check_ali___> Checking the sequence-structure alignment. 

Implied target CA(i)-CA(i+1) distances longer than  8.0 angstroms:

ALN_POS  TMPL  RID1  RID2  NAM1  NAM2     DIST
----------------------------------------------
END OF TABLE
read_to_681_> topology.submodel read from topology file:        3

getf_______W> RTF restraint not found in the atoms list:
              residue type, indices:    18   148
              atom names           : C     +N
              atom indices         :  1155     0

getf_______W> RTF restraint not found in the atoms list:
              residue type, indices:    18   148
              atom names           : C     CA    +N    O
              atom indices         :  1155  1151     0  1156
fndatmi_285W> Only      129 residues out of      132 contain atoms of type  CA
              (This is usually caused by non-standard residues, such
              as ligands, or by PDB files with missing atoms.)
patch_s_522_> Number of disulfides patched in MODEL:        4
mdtrsr__446W> A potential that relies on one protein is used, yet you have at
              least one known structure available. MDT, not library, potential is used.
iup2crm_280W> No topology library in memory or assigning a BLK residue.
              Default CHARMM atom type assigned:  C1 -->  CT2
              This message is written only for the first such atom.
0 atoms in HETATM/BLK residues constrained
to protein atoms within 2.30 angstroms
and protein CA atoms within 10.00 angstroms
0 atoms in residues without defined topology
constrained to be rigid bodies
condens_443_> Restraints marked for deletion were removed.
              Total number of restraints before, now:    13194    12190


>> ENERGY; Differences between the model's features and restraints:
Number of all residues in MODEL                   :      148
Number of all, selected real atoms                :     1157    1157
Number of all, selected pseudo atoms              :        0       0
Number of all static, selected restraints         :    12190   12190
COVALENT_CYS                                      :        F
NONBONDED_SEL_ATOMS                               :        1
Number of non-bonded pairs (excluding 1-2,1-3,1-4):     2351
Dynamic pairs routine                             : 2, NATM x NATM cell sorting
Atomic shift for contacts update (UPDATE_DYNAMIC) :    0.390
LENNARD_JONES_SWITCH                              :    6.500   7.500
COULOMB_JONES_SWITCH                              :    6.500   7.500
RESIDUE_SPAN_RANGE                                :        0   99999
NLOGN_USE                                         :       15
CONTACT_SHELL                                     :    4.000
DYNAMIC_PAIRS,_SPHERE,_COULOMB,_LENNARD,_MODELLER :        T       T       F       F       F
SPHERE_STDV                                       :    0.050
RADII_FACTOR                                      :    0.820
Current energy                                    :         804.2236





Summary of the restraint violations: 

   NUM     ... number of restraints.
   NUMVI   ... number of restraints with RVIOL > VIOL_REPORT_CUT[i].
   RVIOL   ... relative difference from the best value.
   NUMVP   ... number of restraints with -Ln(pdf) > VIOL_REPORT_CUT2[i].
   RMS_1   ... RMS(feature, minimally_violated_basis_restraint, NUMB).
   RMS_2   ... RMS(feature, best_value, NUMB).
   MOL.PDF ... scaled contribution to -Ln(Molecular pdf).

 #                     RESTRAINT_GROUP      NUM   NUMVI  NUMVP   RMS_1   RMS_2         MOL.PDF     S_i
------------------------------------------------------------------------------------------------------
 1 Bond length potential              :    1177       0      0   0.006   0.006      11.387       1.000
 2 Bond angle potential               :    1592       0      7   2.150   2.150      147.53       1.000
 3 Stereochemical cosine torsion poten:     738       0     24  47.751  47.751      266.00       1.000
 4 Stereochemical improper torsion pot:     481       0      0   1.402   1.402      20.918       1.000
 5 Soft-sphere overlap restraints     :    2351       0      0   0.002   0.002     0.65140       1.000
 6 Lennard-Jones 6-12 potential       :       0       0      0   0.000   0.000      0.0000       1.000
 7 Coulomb point-point electrostatic p:       0       0      0   0.000   0.000      0.0000       1.000
 8 H-bonding potential                :       0       0      0   0.000   0.000      0.0000       1.000
 9 Distance restraints 1 (CA-CA)      :    2405       0      0   0.131   0.131      40.957       1.000
10 Distance restraints 2 (N-O)        :    2564       0      0   0.155   0.155      68.747       1.000
11 Mainchain Phi dihedral restraints  :       0       0      0   0.000   0.000      0.0000       1.000
12 Mainchain Psi dihedral restraints  :       0       0      0   0.000   0.000      0.0000       1.000
13 Mainchain Omega dihedral restraints:     147       0      3   4.636   4.636      37.258       1.000
14 Sidechain Chi_1 dihedral restraints:     120       0      0  68.491  68.491      25.477       1.000
15 Sidechain Chi_2 dihedral restraints:      86       0      0  63.198  63.198      33.297       1.000
16 Sidechain Chi_3 dihedral restraints:      35       0      0  89.980  89.980      22.231       1.000
17 Sidechain Chi_4 dihedral restraints:      20       0      0 102.640 102.640      12.968       1.000
18 Disulfide distance restraints      :       4       0      0   0.006   0.006     0.27561E-01   1.000
19 Disulfide angle restraints         :       8       0      0   2.789   2.789      1.3737       1.000
20 Disulfide dihedral angle restraints:       4       0      0  23.141  23.141      2.0451       1.000
21 Lower bound distance restraints    :       0       0      0   0.000   0.000      0.0000       1.000
22 Upper bound distance restraints    :       0       0      0   0.000   0.000      0.0000       1.000
23 Distance restraints 3 (SDCH-MNCH)  :    1691       0      0   0.425   0.425      34.838       1.000
24 Sidechain Chi_5 dihedral restraints:       0       0      0   0.000   0.000      0.0000       1.000
25 Phi/Psi pair of dihedral restraints:     146      16     13  23.998  62.053      26.471       1.000
26 Distance restraints 4 (SDCH-SDCH)  :     972       0      0   0.637   0.637      52.045       1.000
27 Distance restraints 5 (X-Y)        :       0       0      0   0.000   0.000      0.0000       1.000
28 NMR distance restraints 6 (X-Y)    :       0       0      0   0.000   0.000      0.0000       1.000
29 NMR distance restraints 7 (X-Y)    :       0       0      0   0.000   0.000      0.0000       1.000
30 Minimal distance restraints        :       0       0      0   0.000   0.000      0.0000       1.000
31 Non-bonded restraints              :       0       0      0   0.000   0.000      0.0000       1.000
32 Atomic accessibility restraints    :       0       0      0   0.000   0.000      0.0000       1.000
33 Atomic density restraints          :       0       0      0   0.000   0.000      0.0000       1.000
34 Absolute position restraints       :       0       0      0   0.000   0.000      0.0000       1.000
35 Dihedral angle difference restraint:       0       0      0   0.000   0.000      0.0000       1.000
36 GBSA implicit solvent potential    :       0       0      0   0.000   0.000      0.0000       1.000
37 EM density fitting potential       :       0       0      0   0.000   0.000      0.0000       1.000
38 SAXS restraints                    :       0       0      0   0.000   0.000      0.0000       1.000
39 Symmetry restraints                :       0       0      0   0.000   0.000      0.0000       1.000



# Heavy relative violation of each residue is written to: LYSC_HUMAN.V99990001
# The profile is NOT normalized by the number of restraints.
# The profiles are smoothed over a window of residues:    1
# The sum of all numbers in the file:   14675.1133



List of the violated restraints:
   A restraint is violated when the relative difference
   from the best value (RVIOL) is larger than CUTOFF.

   ICSR   ... index of a restraint in the current set.
   RESNO  ... residue numbers of the first two atoms.
   ATM    ... IUPAC atom names of the first two atoms.
   FEAT   ... the value of the feature in the model.
   restr  ... the mean of the basis restraint with the smallest
              difference from the model (local minimum).
   viol   ... difference from the local minimum.
   rviol  ... relative difference from the local minimum.
   RESTR  ... the best value (global minimum).
   VIOL   ... difference from the best value.
   RVIOL  ... relative difference from the best value.


-------------------------------------------------------------------------------------------------

Feature 25                           : Phi/Psi pair of dihedral restraints     
List of the RVIOL violations larger than   :       6.5000

    #   ICSR  RESNO1/2 ATM1/2   INDATM1/2    FEAT   restr    viol   rviol   RESTR    VIOL   RVIOL
    1   4005   1M   2K C   N       7    9  -64.57  -70.20    5.95    0.39  -62.90  179.26   23.19
    1          2K   2K N   CA      9   10  138.46  140.40                  -40.80
    2   4006   2K   3A C   N      16   18  -74.64  -68.20    9.27    0.58  -62.50  167.58   28.07
    2          3A   3A N   CA     18   19  151.96  145.30                  -40.90
    3   4007   3A   4L C   N      21   23  -98.45 -108.50   13.80    0.77  -63.50 -179.75   27.15
    3          4L   4L N   CA     23   24  141.97  132.50                  -41.20
    4   4010   6V   7L C   N      44   46 -103.37 -108.50    6.06    0.27  -63.50  175.07   22.35
    4          7L   7L N   CA     46   47  129.27  132.50                  -41.20
    5   4011   7L   8G C   N      52   54  -90.16  -80.20   68.15    3.19   82.20 -161.64    8.57
    5          8G   8G N   CA     54   55  106.68  174.10                    8.50
    6   4012   8G   9L C   N      56   58 -125.41 -108.50   41.06    2.01  -63.50  161.23   26.22
    6          9L   9L N   CA     58   59  169.93  132.50                  -41.20
    7   4013   9L  10V C   N      64   66  -58.86  -62.40    4.38    0.54 -125.40 -175.85   10.37
    7         10V  10V N   CA     66   67  -44.98  -42.40                  143.30
    8   4014  10V  11L C   N      71   73 -111.10 -108.50   18.88    0.99  -63.50  174.23   27.12
    8         11L  11L N   CA     73   74  151.20  132.50                  -41.20
    9   4015  11L  12L C   N      79   81  -72.91  -70.70   13.00    1.07  -63.50  170.25   23.12
    9         12L  12L N   CA     81   82  128.79  141.60                  -41.20
   10   4016  12L  13S C   N      87   89  163.67 -136.60   83.64    2.69  -64.10  175.41   22.09
   10         13S  13S N   CA     89   90 -150.26  151.20                  -35.00
   11   4018  14V  15T C   N     100  102 -104.81 -124.80   23.58    0.88  -63.20  178.03   20.35
   11         15T  15T N   CA    102  103  131.00  143.50                  -42.10
   12   4020  16V  17Q C   N     114  116   88.80 -121.10  150.10    5.65  -63.80 -124.45   42.00
   12         17Q  17Q N   CA    116  117  139.14  139.70                  -40.30
   13   4021  17Q  18G C   N     123  125 -101.84  -80.20   44.39    1.04   82.20 -125.09    7.89
   13         18G  18G N   CA    125  126 -147.14  174.10                    8.50
   14   4022  18G  19K C   N     127  129  -77.16  -70.20   13.93    1.14  -62.90  169.73   21.30
   14         19K  19K N   CA    129  130  128.33  140.40                  -40.80
   15   4057  53E  54S C   N     415  417 -136.30  -64.10   78.23    8.47  -64.10   78.23    8.47
   15         54S  54S N   CA    417  418   -4.87  -35.00                  -35.00
   16   4069  65A  66G C   N     508  510  -66.26  -62.40    4.87    0.92   82.20  157.53   12.01
   16         66G  66G N   CA    510  511  -44.18  -41.20                    8.50


report______> Distribution of short non-bonded contacts:


DISTANCE1:  0.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40
DISTANCE2:  2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50
FREQUENCY:     0    0    0    0    0    4   10   37   67  104  113  140  158  187  183


<< end of ENERGY.

>> Summary of successfully produced models:
Filename                          molpdf
----------------------------------------
LYSC_HUMAN.B99990001.pdb       804.22363

In [2]:
import nglview
import ipywidgets
In [5]:
w1 = nglview.show_structure_file('LYSC_HUMAN.B99990001.pdb')
In [7]:
w2 = nglview.show_pdbid('1lmp')
In [14]:
class mymodel(modeller.automodel.automodel):
    def special_restraints(self, aln):
        rsr = self.restraints
        for ids in (('OD1:98:A', 'O6A:131:A'),
                    ('N:65:A', 'O7B:132:A'),
                    ('OD2:73:A', 'O1C:133:A')):
                    atoms = [self.atoms[i] for i in ids]
                    rsr.add(forms.upper_bound(group=physical.upper_distance,
                      feature=features.distance(*atoms), mean=3.5, stdev=0.1)) 
In [ ]: