import sys 
import modeller 
import _modeller
import modeller.automodel 
import IPython.display

from IPython.display import display,Image

env=modeller.environ()
env.io.hetatm=True

Лизоцим форели, на основе которого будем проводить гомологичное моделирование.

!wget http://www.pdb.org/pdb/files/1lmp.pdb

--2017-12-11 00:30:50--  http://www.pdb.org/pdb/files/1lmp.pdb
Resolving www.pdb.org (www.pdb.org)... 128.6.244.52
Connecting to www.pdb.org (www.pdb.org)|128.6.244.52|:80... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://www.rcsb.org/pdb/files/1lmp.pdb [following]
--2017-12-11 00:30:51--  https://www.rcsb.org/pdb/files/1lmp.pdb
Resolving www.rcsb.org (www.rcsb.org)... 132.249.213.110
Connecting to www.rcsb.org (www.rcsb.org)|132.249.213.110|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: http://files.rcsb.org/view/1lmp.pdb [following]
--2017-12-11 00:30:52--  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'

    [   <=>                                 ] 131,301      137K/s   in 0.9s    

2017-12-11 00:30:53 (137 KB/s) - `1lmp.pdb' saved [131301]

Возьмем последовательность белка лизоцима из организма Crassostrea gigas.

! wget http://www.uniprot.org/uniprot/Q6L6Q6.fasta

--2017-12-11 00:31:04--  http://www.uniprot.org/uniprot/Q6L6Q6.fasta
Resolving www.uniprot.org (www.uniprot.org)... 193.62.193.81, 128.175.240.211
Connecting to www.uniprot.org (www.uniprot.org)|193.62.193.81|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 210 [text/plain]
Saving to: `Q6L6Q6.fasta'

100%[======================================>] 210         --.-K/s   in 0s      

2017-12-11 00:31:04 (25.2 MB/s) - `Q6L6Q6.fasta' saved [210/210]

Создаем объект выравнивание.

alignm=modeller.alignment(env)

alignm.append(file='Q6L6Q6.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')

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.

alignm[0].code='Q6L6Q6_protein' # поправляем идентификаторы

Делаем выравнивание и сохраняем

alignm.salign()
alignm.write(file='all_in_one.ali', alignment_format='PIR')

SALIGN_____> adding the next group to the alignment; iteration    1

Строим модель

s = alignm[0]
pdb = alignm[1]

print s.code, pdb.code
 
## Создаем объект automodel
a = modeller.automodel.automodel(env, alnfile='all_in_one.ali', knowns= pdb.code, sequence = s.code );

a.name='mod'+s.code
a.starting_model = 1
a.ending_model = 2
a.make()

Получилось 2 модели. Нарисуем одну.

import nglview
import ipywidgets

w2=nglview.NGLWidget()
w2.add_component('Q6L6Q6_protein.B99990002.pdb')
display(Image('lig3.png'))

Лиганд не появился. Где он? Проверим список остатков. Там нет лиганда. Добавим его

alignm[0].residues[-3:]

[Residue VAL, Residue ASN, Residue SER]

blk=[i.code for i in alignm[0].residues]
blk=''.join(blk)+'-...' # добавили 3 остатка к последовательности

alignm.append_sequence(blk)

s = alignm[2]
pdb = alignm[1]
s.code = 'with_legand'

alignm.salign()
alignm.write(file='all_in_one.ali', alignment_format='PIR')

print s.code, pdb.code

## Создаем объект automodel
a = modeller.automodel.automodel(env, alnfile='all_in_one.ali', knowns= pdb.code , sequence = s.code )

a.name='new'+s.code
a.starting_model = 1
a.ending_model = 2
a.make()

Лиганд на месте

import nglview
import ipywidgets
w1 = nglview.show_structure_file('with_legand.B99990001.pdb')
display(Image('lig2.png'))

Теперь мы хотим перместить лиганд в другое место: определяем место, меняем базовые рестрейны, добавляем новые.

Нужно поменять рестрейн на расстояние до Ca и на расстояние между отдельными остатками лиганда (иначе он "гнется")

class mymodel(modeller.automodel.automodel):
    def special_restraints(self, aln):#новые рестрейны
        
        rsr = self.restraints
        
        at = self.atoms
        
        for x,y in [('CA:15','O6:139')]:
            rsr.add(modeller.forms.gaussian(group=modeller.physical.xy_distance, 
                                        feature=modeller.features.distance(
                                        at[x],at[y]),mean=3.0, stdev=0.1))
 
            
    def nonstd_restraints(self, aln): # перезадаем функцию рестрейнов
        
       #"""Create restraints on HETATM and BLK residues."""
       # Select all HETATM residues plus any ATOM residues that have
       # no defined topology (generally speaking, BLK residues)
        allatoms = selection(self)
        selhet = allatoms.only_het_residues() | allatoms.only_no_topology()

        rsrgrp = physical.xy_distance
        self.het_std_restraints(aln, selhet, 0.0001, 2.3, rsrgrp)
        self.het_het_restraints(aln, selhet, 10.0, 2.3, rsrgrp)
        self.het_internal_restraints(aln, selhet, rsrgrp)

from modeller import *
from modeller.automodel import * 


a = mymodel(env,  alnfile='all_in_one.ali', knowns= pdb.code , sequence = s.code)
a.make()

Лиганд поменял место посадки (забился вглубь белка).

import nglview
import ipywidgets
w1 = nglview.show_structure_file('with_legand.B99990001.pdb')
display(Image('lig1.png'))