Работа с Pymol

from xmlrpclib import ServerProxy
from IPython.display import display,Image
cmd = ServerProxy(uri = "http://localhost:8888/RPC2")

import os, sys, time

# PyMOL launching
import __main__

# Tell PyMOL to launch quiet (-q), fullscreen (-e), without internal GUI (-i)
# and without GUI (-c) and with arguments from command line
#__main__.pymol_argv = [ 'pymol', '-qei' ]
#__main__.pymol_argv = [ 'pymol', '-cp' ]
__main__.pymol_argv = [ 'pymol', '-x' ]  

import pymol
 
# Call the function below before using any PyMOL modules.
pymol.finish_launching()
 
from pymol import cmd, stored

import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import RDConfig
from rdkit.Chem.Draw import ShowMol
from rdkit.Chem.Draw import IPythonConsole 
from rdkit.Chem import Draw
import numpy as np

Sculpting

Изменение структуры белка 1LMC

Wizard->Demo->Sculpting

Sculpting является real-time-оптимизатором геометрии молекулы при ее изменении ее пространственной структуры: данный алгоритм старается сохранить значения длины связей, углов, торсионных углов при изменении координат. Изменение структуры проводилось пошагово:

• Выбираем Sculpting из меню Wizard
• Выбираем центральный атом для модификаций

  Ctrl-middle-click

• Тянем атом в любую сторону

  Ctrl-left-click-and-drag

x = Image(filename='1lmp_before.png')
y = Image(filename='1lmp_after.png')
display(x, y)

Mutagenesis

Средствами Tcl/Tk интерфейса:

Wizard->Mutagenesis

была проведена одиночная мутация в 1lmp, приводящая к потере связывания белка с лигандом. Команды:

Action->find->polar contacts->within selection

показали связи между лигандом и белком. Проводилась мутация ASP-101 на TRP-101. После мутирования связь между 101 остатком белка и лигандом не образовалась.

Подписать названия остатков можно либо в меню объекта/выборки (Object Control Panel):

L->residues,

предварительно выделив необходимый остаток, либо с помощью ввода команд:

label n. CA and 1lmp and i. 101, "%s %s" % (resn, resi) label n. CH2 and /mutagenesis and i. 101, "%s %s" % (resn, resi)

Желтым кругом выделено место мутации.

Image(filename='mutagenesis.jpg')

На двух нижепредставленных картинках изображен остаток 101: сверху до мутации, снизу - после.

x = Image(filename='ASP-101.png')
y = Image(filename='TRP-101.png')
display(x, y)

Далее изображено выравнивание изначальной и мутированной структур: красный - TRP-101, желтый - ASP-101. Структуры накладывались друг на друга с помощью команды align:

• Performs a sequence alignment followed by a structural superposition, and then carries out zero or more cycles of refinement in order to reject structural outliers found during the fit. Its does a good job on proteins with decent sequence similarity (identity >30%). For comparing proteins with lower sequence identity, the super and cealign commands perform better

align

Команды cealign, super, fit, pair_fir в данном случае не подошли, т. к. последовательности двух белков отличались только одной аминокислотой.

Image(filename='alignment.png')

Ролик был создан с использованием команд: mset, mview, align, translate. Получившаяся анимация была сохранена как фильм:

File->Export Movie As->MPEG

Совмещение белков производится командой align по вышеуказанным причинам. Место мутации показано серией set_view положений, перехожы между которыми сглажены интерполяцией. Ниже приведены указания из раздела object motions PyMOLWiki.

Object Motions

Now that we're experts at moving the PyMOL camera around, let's start moving objects around while keeping the camera steady. To do this you must have matrix_mode set to 1, otherwise PyMOL won't save your object's repositioning.

• Use the mouse to get the 'right orientation & zoom'. Then use get_view to get the view matrix. Finally, store that camera-position view matrix in your script. Works every time.
• For object motions, the command translate [-70,70,70], object=AA would be the same as using the mouse and moving the object AA -70 units on the X-axis, +70 units on the Y and 70 on the Z. If you don't use the object= you will not get the desired effect.
• For the above 'explosion' you can get quick motions by interpolating a large change over just a few frames.

cmd.do('''
load 1lmp.pdb
remove solvent
hide everything
show sticks
show spheres, 1lmp
set stick radius, .07, 1lmp
set sphere scale, .18, 1lmp
set sphere scale, .13, elem H
set bg rgb=[1, 1, 1]
set stick quality, 50, 1lmp
set sphere quality, 4, 1lmp
color gray85, elem C
color red, elem O
color slate, elem N
color gray98, elem H
show cartoon
set stick color, black, 1lmp
set ray trace mode, 1
set ray texture, 2
set antialias, 3
set ambient, 0.5
set spec count, 5
set shininess, 50
set specular, 1
set reflect, .1
set dash gap, 0
set dash color, black
set dash gap, .15
set dash length, .05
set dash round ends, 0
set dash radius, .05

load mutagenesis.pdb
preset.publication("mutagenesis", self=cmd)

extract ligands, het 
hide everything, ligands
show sticks, ligands
color deepblue, ligands

set retain order
set ray trace mode, 1
set ray texture, 2
set antialias, 3
set ambient, 0.5
set spec count, 5
set shininess, 50
set specular, 1
set reflect, .1
set dash gap, 0
set dash color, black
set dash gap, .15
set dash length, .05
set dash round ends, 0
set dash radius, .05



set retain order
set matrix mode, 1
set movie panel, 1

mset 1 x1100
frame 1
translate [-50, 0, 0], 1lmp
orient 
mview store
mview store, object=1lmp
mview store, object=mutagenesis

frame 100
align 1lmp, mutagenesis
mview store
mview store, object=1lmp
mview store, object=mutagenesis

frame 200
set view (\
     0.310487568,   -0.351513296,   -0.883196354,\
     0.072738476,    0.935177326,   -0.346630543,\
     0.947790325,    0.043382134,    0.315929443,\
     0.000000000,    0.000000000, -126.023445129,\
    14.233245850,   51.144100189,   20.735301971,\
    99.357849121,  152.689041138,  -20.000000000 )
mview store

frame 250
zoom i. 101
mview store

frame 300
color yellow, /1lmp//A/ASP`101
mview store

frame 350
color tv red, /mutagenesis//A/TRP`101
show sticks, /mutagenesis//A/TRP`101
set view (\
     0.027544985,    0.392354786,   -0.919401050,\
     0.185692057,    0.901733696,    0.390379518,\
     0.982219279,   -0.181478143,   -0.048018824,\
     0.000001287,    0.000002790,  -37.680519104,\
     4.573862553,   56.978332520,   28.689073563,\
    18.598453522,   56.762577057,  -20.000000000 )
mview store

frame 400
select /1lmp//A/ASP`101/OD2 or /ligands/B/A/NAG`130/O6
hide everything, mutagenesis
mview store

frame 450
  set view (\
      0.874671102,    0.465957344,   -0.133546099,\
     -0.245426506,    0.663314044,    0.706953108,\
      0.417990744,   -0.585572362,    0.694536507,\
     -0.000008550,    0.000009350,  -14.178204536,\
      7.696578026,   58.481658936,   32.177547455,\
     -4.169544697,   32.525955200,  -20.000000000 )
mview store

frame 500
show cartoon, mutagenesis
spectrum count, rainbow, mutagenesis
color tv red, /mutagenesis//A/TRP`101
show sticks, /mutagenesis//A/TRP`101
orient /mutagenesis//A/TRP`101
  set view (\
      0.081681184,    0.597173035,   -0.797941208,\
      0.041003622,    0.797925293,    0.601360202,\
      0.995813429,   -0.081837848,    0.040691219,\
      0.000000000,    0.000000000,  -23.100343704,\
      4.899285793,   58.137855530,   30.553213120,\
     18.212486267,   27.988201141,  -20.000000000 )
mview store

frame 550
  set view (\
      0.640822947,    0.750949860,   -0.159429103,\
     -0.407308728,    0.508619845,    0.758554935,\
      0.650723577,   -0.421163052,    0.631805360,\
      0.000000000,    0.000000000,  -23.100343704,\
      4.899285793,   58.137855530,   30.553213120,\
     18.212486267,   27.988201141,  -20.000000000 )
mview store

frame 600
  set view (\
      0.423352420,    0.569107592,    0.704899907,\
     -0.893512607,    0.390838474,    0.221085414,\
     -0.149684191,   -0.723433375,    0.673969448,\
     -0.000007562,    0.000001580,  -23.100343704,\
      5.659269810,   59.242378235,   31.976255417,\
      2.682543993,   43.518146515,  -20.000000000 )
mview store

frame 650
  set view (\
     -0.111630894,    0.917784572,    0.381058097,\
     -0.991070449,   -0.074678324,   -0.110466041,\
     -0.072931081,   -0.389986157,    0.917925060,\
     -0.000007562,    0.000001580,  -23.100343704,\
      5.659269810,   59.242378235,   31.976255417,\
     15.350089073,   30.850589752,  -20.000000000 )
mview store
frame 700
  set view (\
     -0.750151157,    0.503197253,    0.429022849,\
      0.061019193,    0.698701203,   -0.712805569,\
     -0.658440828,   -0.508530557,   -0.554838359,\
     -0.000007562,    0.000001580,  -23.100343704,\
      5.659269810,   59.242378235,   31.976255417,\
     15.350089073,   30.850589752,  -20.000000000 )
mview store

frame 750
  set view (\
     -0.965006828,    0.223388940,    0.137302488,\
      0.242283061,    0.559418619,    0.792683423,\
      0.100267045,    0.798214555,   -0.593964756,\
     -0.000005804,    0.000001609,  -23.100343704,\
      5.487463474,   59.898345947,   30.582403183,\
     10.618468285,   35.582187653,  -20.000000000 )
mview store

frame 800
  set view (\
      0.448935241,    0.174991682,   -0.876257598,\
     -0.879195750,    0.261651903,   -0.398186386,\
      0.159592822,    0.949162602,    0.271315128,\
     -0.000005804,    0.000001609,  -23.100343704,\
      5.487463474,   59.898345947,   30.582403183,\
     10.618468285,   35.582187653,  -20.000000000 )
mview store

frame 850
  set view (\
      0.448935241,    0.174991682,   -0.876257598,\
     -0.879195750,    0.261651903,   -0.398186386,\
      0.159592822,    0.949162602,    0.271315128,\
     -0.000003759,   -0.000015527,  -23.100343704,\
     -6.543584824,   54.431243896,   34.307575226,\
     -1.223415494,   47.423828125,  -20.000000000 )
hide labels
mview store

frame 900
set view (\
     0.310487568,   -0.351513296,   -0.883196354,\
     0.072738476,    0.935177326,   -0.346630543,\
     0.947790325,    0.043382134,    0.315929443,\
     0.000000000,    0.000000000, -126.023445129,\
    14.233245850,   51.144100189,   20.735301971,\
    99.357849121,  152.689041138,  -20.000000000 )
mview store
mview store, object=1lmp
mview store, object=mutagenesis

mview interpolate''')

from IPython.display import HTML
from base64 import b64encode


video = open('mutagenesis_movie.mp4', 'rb').read()
video_encoded = b64encode(video)
video_tag = '<video controls width="960" height="540" alt = "PyMol Movie" src ="data:video/mp4;base64,{0}" type = "video/mp4">'.format(video_encoded)
HTML(data = video_tag)

Метка TAMRA

Присоединение флуоресцентной метки TAMRA к белку через сложноэфирную связь при помощи команд:

• Joins two objects into one by forming a bond: a copy of the object containing the first atom is moved so as to form an approximately reasonable bond with the second, and is then merged with the first object

  fuse

• Rotates the torsion on the bond currently picked for editing (the rotated fragment will correspond to the first atom specified when picking the bond or the nearest atom, if picked using the mouse)

  torsion

Для выбора оптимальноого расположения метки могла использоваться команда torsion 15, однако эього не потребовалось делать.

Сложные эфиры - производные оксокислот (гетерофункциональные соединения, содержащие карбоксильную и карбонильную, т. е. альдегидную или кетонную, группы) с общей формулой RkE(=O)l(OH)m, где l ≠ 0, формально являющиеся продуктами замещения атомов водорода в гидроксилах —OH кислотной функции на углеводородный остаток (алифатический, алкенильный, ароматический или гетероароматический); рассматриваются также как ацилпроизводные спиртов. Следовательно, сложноэфироная связь между белком и TAMRA образуется через OH-группы β-оксоаминокислот (серина и треонина). Присоеденим метку к SER-122.

Получим PDB-файл метки из SMILES:

TAMRA = "CN(C)C1=CC2=C(C=C1)C(=C3C=CC(=[N+](C)C)C=C3O2)C4=C(C=CC(=C4)C(=O)[O-])C(=O)O"
mol = Chem.MolFromSmiles(TAMRA)
AllChem.Compute2DCoords(mol)
img = Chem.Draw.MolToImageFile(mol, 'TAMRA.png', size=(400, 400))

pdb_block = rdkit.Chem.rdmolfiles.MolToPDBFile(mol, 'tamra.pdb')
tamra_file = open('tamra.pdb', 'r')
print(tamra_file.read())

x = Image(filename='TAMRA_SMILES.png')
y = Image(filename='TAMRA_PDB.png', width='400')
display(x, y)

Присоединим метку к остатку лизоцима.

cmd.do('''
select 1lmp_O, /1lmp/A/A/SER`122/OG
select tamra_C, /tamra///UNL`1/C14
fuse 1lmp_O, tamra_C''')

SER-122 (желтый) с меткой TAMRA (сиреневая).

Image(filename='fuse.png')

Поли-цепь

x = Image(filename='ramachandran_plot.jpg')
y = Image(filename='ramachandran_amino_acids.jpg')
z = Image(filename='dihedral_angles.jpg')
display (x, y, z)

Необходимые значения углов были взяты из вышепреведенных изображений. Первый скрипт строит полаланиновую цепь длинной 100 аминокислотных остатков; второй - цепь из чередующихся серинов и треонинов.

Потребовавшиеся команды:

• Picks an atom or bond for editing

  edit

• Adds a whole residue (amino acid) from the library

  editor.attach_amino_acid

residues =  "ALA"

cmd.reinitialize()
cmd.bg_color('white')
cmd.fragment(residues)


phi = '-57'
psi = '-47'
n = 100

def polyhelix(n, residues, phi=None, psi=None):
# At first, just add all the aminoacids
   # for res in residues:
  
    for i in range(2, n+1):
        # Choose the atom to attach a new residue to
            cmd.edit("i. %i & n. C" % i)  
            cmd.do('editor.attach_amino_acid("pk1", "%s")' % residues)
            time.sleep(1)  

# Choose the atom to attach a new residue to
     
    for i in range(2, n+1):
        if i % 2 == 0:
            cmd.set_dihedral("i. %i and n. N"  % i, 
                             "i. %i and n. CA" % i,
                             "i. %i and n. C"  % i, 
                             "i. %i and n. N"  % (i+1),
                             phi)
            cmd.set_dihedral("i. %i & n. C"  % i,     
                             "i. %i & n. N"  % (i+1),
                             "i. %i & n. CA" % (i+1), 
                             "i. %i & n. C"  % (i+1),
                             psi)
    cmd.show('cartoon')

polyhelix(100, residues, phi=phi, psi=psi)

Image(filename='ala_all.png')

Image(filename='ala_zoom.png')

cmd.reinitialize()
cmd.bg_color('white')
cmd.fragment('THR')   # The first residue of the future helix

residues =  ("THR", "SER")   # These two residues 
phi = ('-69.6', '-66.5')
psi = ('-35.0', '-33.5')
n = 100

def polyhelix(n, residues, phi=None, psi=None):
# At first, just add all the aminoacids
   # for res in residues:
  
    for i in range(2, n+1, 1):
        if i % 2 == 0:
        # Choose the atom to attach a new residue to
            cmd.edit("/THR///THR`%i/C" %i)  
            cmd.do('editor.attach_amino_acid("pk1", "%s")' % residues[1])
            time.sleep(1)  
        if not i % 2 == 0:
            cmd.edit("/THR///SER`%i/C" %i) 
            cmd.do('editor.attach_amino_acid("pk1", "%s")' % residues[0])
            time.sleep(1)

#    for i in range(3, n+1, 2):
       # Choose the atom to attach a new residue to
     
    for i in range(2, n+1, 1):
        if i % 2 == 0:
            cmd.set_dihedral("i. %i and n. N"  % i, 
                             "i. %i and n. CA" % i,
                             "i. %i and n. C"  % i, 
                             "i. %i and n. N"  % (i+1),
                             phi[0])
            cmd.set_dihedral("i. %i & n. C"  % i,     
                             "i. %i & n. N"  % (i+1),
                             "i. %i & n. CA" % (i+1), 
                             "i. %i & n. C"  % (i+1),
                             psi[0])
        if not i % 2 == 0:
            cmd.set_dihedral("i. %i and n. N"  % i, 
                             "i. %i and n. CA" % i,
                             "i. %i and n. C"  % i, 
                             "i. %i and n. N"  % (i+1),
                             phi[1])
            cmd.set_dihedral("i. %i & n. C"  % i,     
                             "i. %i & n. N"  % (i+1),
                             "i. %i & n. CA" % (i+1), 
                             "i. %i & n. C"  % (i+1),
                             psi[1])

polyhelix(100, residues, phi=phi, psi=psi)

Image(filename='polyhelix_all.png')

Image(filename='polyhelix_zoom.png')

B-форма ДНК

В связи с тем, что в PyMOL нет заготовок для нуклеиновых кислот нет, невозможно пострить B-формы ДНК длинной 100 пар нуклеотидов теми же командами, что и полипептидную цепь. Для решения поставленной задачи нужно определить матрицу превращения при совмещении одной пары нуклеотидов с последующей. После этого можно множить эту пару и применять к ней матрицу превращения. Стартовую пару пару возьмем из середины цепи. Код был написан на основе псевдокода:

fetch perfectdna
create pair create nextpair
create 1, pair
python cmd.pair_fit("pair","nextpair")
trans=cmd.get_object_matrix("pair")
for x in range(2,500):

cmd.create x,x-1
cmd.transform_selection( str(x), trans, homogenous=0 )
python end

cmd.do('''
reinitialize 
bg_color white
fetch 1bna
remove solvent
create pair, c. A and i. 5 + c. B and i. 20
create nextpair, c. A and i. 4 + c. B and i. 21
create 1, pair''')
python
cmd.pair_fit('pair', 'nextpair')
trans = cmd.get_object_matrix('pair')
for x in range(2, 101):
    cmd.create (str(x), str(x-1))
    cmd.transform_selection(str(x), trans, homogenous = 0)
    time.sleep(1)
# python end

x = Image(filename='100bna_zoom.png')
y = Image(filename='100bna_tor.png')
display(x, y)

Image(filename='100bna_all.png')

G-квадруплекс и пентаплексная ДНК

Для построения G-квадруплексной и пентаплексную ДНК длинной 100 пар нуклеотидов использовалься тот же принцип, что и при построении B-формы ДНК. Необходимо, чтобы количество цепей равнялось "плетности" молекулы. В качестве основы для построения G-квадруплекса бы взят четырехтяжевый параллельный квдруплекс, фрагмент вирусного генома SV40 - 1EVO. Другой структурой являлась 2LUJ - параллельная олигоизогуаниновая пентаплексная ДНК, сформированная d(T(iG)4T) мотивомю

cmd.do('''
reinitialize 
bg_color white
fetch 1evo
remove solvent
create tetrada, c. P and i. 104 + c. R and i. 404 + c. O and i. 304 + c. Q and i. 204
create nexttetrada, c. P and i. 103 + c. R and i. 403 + c. O and i. 303 + c. Q and i. 203
create 1, tetrada''')
python
cmd.pair_fit('tetrada', 'nexttetrada')
trans = cmd.get_object_matrix('tetrada')
for x in range(2, 101):
    cmd.create (str(x), str(x-1))
    cmd.transform_selection(str(x), trans, homogenous = 0)
    time.sleep(1)
python end

x = Image(filename='quadr_zoom.png')
y = Image(filename='quadr_trs.png')
display(x, y)

cmd.do('''
reinitialize 
bg_color white
fetch 2luj
remove solvent
create penta, c. A and i. 5 + c. E and i. 29 + c. B and i. 11  +  c. D and i. 23 + c. C and i. 17
create nextpenta, c. A and i. 4 + c. E and i. 28 + c. B and i. 10  +  c. D and i. 22 + c. C and i. 16
create 1, penta''')
python
cmd.pair_fit('penta', 'nextpenta')
trans = cmd.get_object_matrix('penta')
for x in range(2, 101):
    cmd.create (str(x), str(x-1))
    cmd.transform_selection(str(x), trans, homogenous = 0)
    time.sleep(1)
python end

x = Image(filename='penta_zoom.png')
y = Image(filename='penta_trs.png')
display(x, y)

Ссылки

• Визуализация: PyMOLWiki

  https://pymolwiki.org/index.php/Main_Page

• Анимация: MovieSchool

  https://www.pymolwiki.org/index.php/MovieSchool

• Логические операторы и алгебра выбора: Selection Algebra

  https://pymolwiki.org/index.php/Selection_Algebra

• Карты Рамачандрана все аминокислот

  Hovmöller, S., Zhou, T. and Ohlson, T. (2002) Conformations of amino acids in proteins

• Значения углов для серина и треонина:

  https://www.sciencedirect.com/science/article/pii/S0006349500765143#tbl1

• Моделирование белков. Принципы организации структуры белков

  Höltje H.-D., Sippl W., Folkers G., Rognan D. (2003) Molecular Modelling. Basic Principles and Applications
  Leach Andrew R. (1996) Molecular Modelling. Priciples and Applications

• Псевдокод

  https://kodomo.fbb.msu.ru/wiki/2009/6/help2