Практикум 2. Работа в PyMol.

Предварительные ласки.

In [21]:
from IPython.display import Image, display
In [2]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [353]:
cmd.reinitialize()
In [354]:
cmd.do('bg_color white; set antialias, 2; set ray_trace_mode, 3;')

Загрузим данные.

In [355]:
cmd.fetch('1lmp')
cmd.do('remove resn hoh')
In [28]:
cmd.do('png 1lmp.png, width=1080, height=1080')
In [29]:
Image('1lmp.png', width=600, height=600)
Out[29]:

Тут мы ручками проскульптировали леую альфа-цепь белка (на картинке выше) и ниже отрисуем сравнение нескульптированного и скульптированного белка.

In [356]:
cmd.do('zoom')
In [31]:
cmd.do('png 1lmp_sculpted.png, width=1080, height=1080')
In [37]:
display(Image('1lmp.png', width=300, height=300), Image('1lmp_sculpted.png', width=300, height=300))

Связывание с лигандом

Перезагрузим модель.

In [12]:
cmd.delete('1lmp')
In [13]:
cmd.fetch('1lmp')
cmd.do('remove resn hoh')
In [37]:
cmd.do('util.cbag 1lmp')

Посмотрим на связывание лигандов из данных PDB.

In [19]:
display(Image('https://cdn.rcsb.org/etl/poseview/img/lm/1lmp/N/NDG/1lmp_NDG.png', width=600),
        Image('https://cdn.rcsb.org/etl/poseview/img/lm/1lmp/N/NAG/1lmp_NAG.png', width=600))

Нас интересуют именно боковые цепи АК (иначе мутагенез так себе скажется). Для NDG перспективными кажутся ASP52 и GLU35. Для NAG кажется перспективным ASP101. Все связи являются водородными.

Отметим, что в модели два NAG (130 и 131) и один NDG (132). Именно NAG130 взаимодействует с ASP101. Согласно данным Ligand viewer на PDB, NAG131 взаимодействует с TRP63 водородной связью. NDG132 взаимодействует ещё и с боковыми остатками ASN46 и VAL109 (якобы гидрофобно). Все взаимодействия указаны для боковых цепей.

Лиганд Остаток Тип
NAG130 ASP101 Водородная связь
NAG131 TRP63 Водородная связь
NDG132 GLU35 Водоробная связь
NDG132 ASN46 Водородная связь
NDG132 ASP52 Водородная связь
NDG132 VAL109 Гидрофобное взаимодействие

Теперь посмотрим на зону контакта.

In [357]:
cmd.do('select ligand, resn NAG or resn NDG\n'
       'select contact_zone, ligand around 3.5\n'
       'select contact_res, byres contact_zone\n'
       'as wire, not ligand\n'
       'show sticks, contact_res\n'
       'distance contact_dist, ligand, contact_zone, 3.5\n'
       'orient ligand or contact_zone\n'
       'rotate x, 90')
In [358]:
cmd.do('label name CA and contact_res, "%s%s" % (resn, resi)')
In [256]:
cmd.do('png all_contacts.png, 1080, 1080')
In [257]:
cmd.do('rotate y, 90')
In [52]:
cmd.do('png all_contacts_y_90.png, 1080, 1080')
In [258]:
cmd.do('rotate y, -90')
In [54]:
display(Image('all_contacts.png', width=600), Image('all_contacts_y_90.png', width=600))

Мы видим, что у нас много просто "контактов" на расстоянии 3.5Å. Якобы гидрофобное взаимодействие с VAL109 не подтвердилось — этот остаток выходит на поверхность. Правда, белок может димеризоваться с собой (но в PDB нет такой информации) или с другими белками (но это мы смотреть не будем). Таким образом, мы дальше будем смотреть только полярные контакты боковых цепей.

In [359]:
cmd.do('hide everything, contact_dist\n'
       'distance contact_dist_polar, ligand, contact_zone and not (name O or name N), 3.5, 2\n'
       'orient ligand or cont\n')
In [121]:
cmd.do('png polar_contacts.png, 1920, 1080')
In [260]:
cmd.do('rotate x, 30')
In [123]:
cmd.do('png polar_contacts_x_30.png, 1920, 1080')
In [261]:
cmd.do('rotate y, -60')
In [125]:
cmd.do('png polar_contacts_y_m_60.png, 1920, 1080')
In [262]:
cmd.do('rotate y, 60')
In [128]:
display(Image('polar_contacts.png', width=720),
        Image('polar_contacts_x_30.png', width=720),
        Image('polar_contacts_y_m_60.png', width=720))

Мы видим два наиболее перспективных взаимодействия: NAG130 с ASP101 и NDG132 с ASP52. Оба взаимодействия — предполагаемые водородные связи с длиной 2.4Å. Остальные взаимодействия больше 3Å.

Ещё посмотрим на красивую поверхность контакта.

In [360]:
cmd.do('show surface, contact_zone\n'
       'hide sticks labels, contact_res and not (resi 101 or resi 52)\n')
In [140]:
cmd.do('png contact_surface.png, 1920, 1080')
In [141]:
display(Image('contact_surface.png', width=720))

Перед тем, чтобы приступить к мутагенезу, надо посмотреть, а как поведёт себя скульптирование без мутаций.

In [264]:
import time
In [265]:
cmd.do('set auto_sculpt, on\n'
       'set sculpting, on\n'
       'sculpt_activate all\n')
time.sleep(60)
cmd.do('set sculpting, off')
In [153]:
cmd.do('png naive_scuplt.png, 1920, 1080')
In [154]:
display(Image('naive_scuplt.png', width=720))

Лиганд-связывающий карман на глаз поменялся незначительно, лиганд не вытолкнулся.

Теперь проведём замену ASP101 на ALA и посмотрим, что получится.

In [266]:
cmd.do('hide everything, contact_dist_pol')
In [267]:
cmd.do('orient ligand or contact_zone\n')
In [270]:
cmd.do('rotate x, 90')
In [159]:
cmd.do('png asp101ala.png, 1920, 1080')
In [160]:
display(Image('asp101ala.png', width=720))

Теперь NAG130 не удерживается водородной связью с ALA101.

In [161]:
cmd.do('set auto_sculpt, on\n'
       'set sculpting, on\n'
       'sculpt_activate all\n')
time.sleep(60)
cmd.do('set sculpting, off')
In [162]:
cmd.do('png asp101ala_after_sculpt.png, 1920, 1080')
In [163]:
display(Image('asp101ala_after_sculpt.png', width=720))

Ничего особенно не изменилось.

Давайте заменим 101 остаток на TRP.

In [164]:
cmd.do('util.cbac resi 101')
In [165]:
cmd.do('orient ligand or contact_zone')
In [166]:
cmd.do('rotate x, 90')
In [167]:
cmd.do('png asp101trp.png, 1920, 1080')
In [168]:
display(Image('asp101trp.png', width=720))

Теперь скульптируем.

In [169]:
cmd.do('set auto_sculpt, on\n'
       'set sculpting, on\n'
       'sculpt_activate all\n')
time.sleep(60)
cmd.do('set sculpting, off')
In [170]:
cmd.do('png asp101trp_after_sculpting.png, 1920, 1080')
In [171]:
display(Image('asp101trp_after_sculpting.png', width=720))

И опять ничего не поменялось!

Вернём ASP на 101 и будем развлекаться с ASP52. Сделаем на этот раз по-красивому, через команды.

In [361]:
cmd.wizard('mutagenesis')
cmd.refresh_wizard()
cmd.get_wizard().do_select('resi 52')
cmd.get_wizard().set_mode('TRP')
cmd.refresh_wizard()
cmd.get_wizard().apply()
cmd.set_wizard()
In [362]:
cmd.do('util.cbag resi 101\n'
       'util.cbac resi 52\n'
       'orient ligand or contact_zone\n')
In [365]:
cmd.rotate('x', 90)

Заменили ASP52 на TRP, далее скульптирование.

In [366]:
cmd.do('set auto_sculpt, on\n'
       'set sculpting, on\n'
       'sculpt_activate all\n')
time.sleep(60)
cmd.do('set sculpting, off')
In [367]:
cmd.do('png asp52trp_after_scuplting.png, 1920, 1080')
In [368]:
display(Image('asp52trp_after_scuplting.png', width=720))

Да тоже ничего не поменялось, в общем-то.

Ну и ладно, оставим так.

Мультик с совмещением мутанта и нативного белка.

Сначала подготовим мутант.

In [276]:
cmd.do('hide surface\n'
       'as cartoon, not resi 52 and not hetatm')
In [277]:
cmd.do('zoom')

Теперь загрузим нативный белок.

In [278]:
cmd.fetch('1lmp', 'native')
cmd.do('remove resn hoh')
In [279]:
cmd.do('util.cbaw native and resi 52')
In [280]:
cmd.do('show sticks, native and resi 52')
In [281]:
cmd.center('1lmp or native')
In [282]:
cmd.zoom('1lmp or native')
In [283]:
cmd.super('1lmp', 'native')
Out[283]:
(0.21096858382225037, 808, 5, 0.42899811267852783, 997, 592.786376953125, 129)
In [284]:
cmd.translate('[-30, 0, 0]', 'native')
In [285]:
cmd.center('1lmp or native')
cmd.zoom('1lmp or native')

Рисуем кадры.

In [286]:
# Начальный кадр
cmd.mset('1 x300')
cmd.frame('1')
cmd.mview('store')
cmd.mview('store', object='native')
In [287]:
# Совмещаем структуры
cmd.frame('100')
cmd.mview('store')
cmd.translate('[30, 0, 0]', object='native')
cmd.mview('store', object='native')
In [288]:
# Центрируем по совмещённым структурам
cmd.frame('150')
cmd.center('native')
cmd.zoom('native')
cmd.mview('store', object='native')
cmd.mview('store')
In [289]:
# Приближаем к лиганду
cmd.frame('200')
cmd.mview('store', object='native')
cmd.center('1lmp and resi 52')
cmd.zoom('1lmp and resi 52')
cmd.mview('store')
In [290]:
# Вращаемся вокруг лиганда
cmd.frame('300')
cmd.mview('store', object='native')
cmd.center('1lmp and resi 52')
cmd.zoom('1lmp and resi 52')
cmd.turn('y', 90)
cmd.mview('store')
In [291]:
# Заполняем кадры
cmd.mview('reinterpolate')
cmd.mview('reinterpolate', object='native_protein')

Сохраняем картинки и собираем гифку, ибо неохота маяться с установкой кодека.

In [295]:
cmd.mpng('movie\\', mode=1, width=480, height=360)
In [297]:
import imageio
In [300]:
with imageio.get_writer('movie.gif', mode='I', fps=30) as writer:
    for i in range(300):
        writer.append_data(imageio.imread(f'movie\\{i + 1:04}.png'))
In [301]:
display(Image('movie.gif', format='png'))

Ура! Положение Cβ почти не изменилось у 52-ого остатка, а дальше наблюдаются изменения. Помимо этого, структуры неполностью совпадают друг с другом, особенно лиганды. Возможно, скульптирование внесло некие изменения, или же мутация повлияла на RMSD.

TAMRA

In [302]:
display(Image('6-Tamra_500.png'))

В метке есть две карбоксильные группы, значит, мы должны цеплять за серин или треонин.

In [303]:
cmd.delete('1lmp')
cmd.delete('ligand')
cmd.delete('contact*')
In [304]:
cmd.zoom()
In [305]:
cmd.show('surface')
cmd.select('potent_zones', 'resn SER+THR and elem O and sidechain')
cmd.color('green', 'potent_zones')
In [306]:
cmd.turn('y', 90)
In [313]:
cmd.label('potent_zones', '"%s%s" % (resn, resi)')
In [320]:
cmd.set('label_position', '(0, 0, 3)')
In [321]:
cmd.do('png potent_zones.png, 1920, 1080')
In [322]:
display(Image('potent_zones.png'))

Посмотрели на белок с противоположной от лиганда стороны и нашли несколько атомов спиртовых групп, выходящих на поверхность.

Мне понравился THR89, будем к нему присоединять.

In [323]:
cmd.hide('surface label', 'all')
In [325]:
cmd.show('sticks', 'resi 89')
cmd.util.cbam('all')
cmd.util.cbag('resi 89')
In [326]:
cmd.zoom('resi 89', buffer=5)

Теперь загрузим метку.

In [327]:
cmd.load('tamra.sdf')

Начнём слияние.

In [329]:
cmd.remove('resi 89 and elem O and sidechain')
cmd.fuse('tamra and donor', 'resi 89 and name CB')
In [333]:
cmd.zoom('resi 89', buffer=10)
In [335]:
cmd.torsion(200)
In [336]:
cmd.unpick()
cmd.delete('tamra')
cmd.remove('elem H')
In [350]:
cmd.zoom('resi 89', buffer=10)
cmd.rotate('y', 90)
In [351]:
cmd.do('png fused.png, 1920, 1080')
In [352]:
display(Image('fused.png'))

Прекрасно: выбрали удачное место, метка ничему не мешает.

Полиаланиновая α-цепь.

Перезапустим среду.

In [395]:
cmd.reinitialize()
In [396]:
cmd.do('bg_color white; set antialias, 2; set ray_trace_mode, 3;')

Начало скрипта.

In [397]:
aa = 'ala'
length = 100
In [398]:
cmd.fragment(aa)
cmd.alter('ala', 'resi=1')
Out[398]:
10
In [399]:
for i in range(length):
    cmd.edit(f'resi {i + 1} and name C')
    cmd.editor.attach_amino_acid('pk1', aa, ss=1)

Конец скрипта.

Нарисуем красивые картинки.

In [400]:
cmd.center('all')
cmd.orient('all')
cmd.zoom('all')
In [402]:
cmd.hide()
cmd.show('sticks')
In [406]:
cmd.do('png poly_ala.png, 1920, 1080')
In [407]:
display(Image('poly_ala.png'))
In [417]:
cmd.rotate('y', 90)
cmd.zoom('ala', buffer=-50)
In [418]:
cmd.do('png poly_ala_inside.png, 1920, 1080')
In [419]:
display(Image('poly_ala_inside.png'))

Получилось красиво. Мы, конечно, верим тому, как pymol строит цепь, но проверим торсионные углы.

In [420]:
def get_phi(start_resi):
    return cmd.get_dihedral(f'resi {start_resi} and name C',
                            f'resi {start_resi + 1} and name N',
                            f'resi {start_resi + 1} and name CA',
                            f'resi {start_resi + 1} and name C')


def get_psi(start_resi):
    return cmd.get_dihedral(f'resi {start_resi} and name N',
                            f'resi {start_resi} and name CA',
                            f'resi {start_resi} and name C',
                            f' resi {start_resi + 1} and name N')
In [421]:
phis = [get_phi(i) for i in range(1, 99)]
psis = [get_psi(i) for i in range(1, 99)]
In [426]:
import matplotlib.pyplot as plt
In [430]:
plt.plot(list(range(1, 99)), phis, label='$\phi$')
plt.plot(list(range(1, 99)), psis, label='$\psi$')
plt.xlabel('Start residue')
plt.ylabel('Dihedral angle')
plt.legend()
plt.show()
In [432]:
phis[0], psis[0]
Out[432]:
(-56.999664306640625, -47.00018310546875)

Все углы φ примерно равны -57°, а все углы ψ примерно равны -47°. Эти данные попадают в зоны стабильной α-цепи на карте Рамачандрана.

B-форма ДНК.

Перезапустим среду.

In [512]:
cmd.reinitialize()
In [513]:
cmd.do('bg_color white; set antialias, 2; set ray_trace_mode, 3;')

Начало скрипта.

  1. Определить наборы нуклеотидов в конкретной модели.
  2. Выбрать две последовательные пары нуклеотидов.
  3. Удаляем старую ДНК.
  4. Совместить пары.
  5. Найти матрицу поворота.
  6. Создаём копии первой пары и поворачиваем их.
  7. Переименовываем номера остатков.
  8. Образуем связи.
  9. Красиво отрисовываем.

Зададим начальные данные.

In [514]:
length = 100

Загрузим структуру B-DNA, чтобы оттуда взять пары и матрицу перехода.

In [515]:
cmd.fetch('355d')
cmd.remove('hetatm')

Посмотрим, какие цепи есть в структуре.

In [516]:
stored.chains = []
cmd.iterate('all', 'stored.chains.append(chain)')
set(stored.chains)
Out[516]:
{'A', 'B'}

Посмотрим, сколько и какие пары есть в структуре.

In [517]:
stored.a, stored.b = [], []
cmd.iterate("chain A and name C1'", 'stored.a.append([resi, resn])')
cmd.iterate("chain B and name C1'", 'stored.b.append([resi, resn])')
Out[517]:
12
In [518]:
pairs = [[stored.a[_], stored.b[11 - _]]
         for _ in range(len(stored.a))]
In [519]:
pairs
Out[519]:
[[['1', 'DC'], ['24', 'DG']],
 [['2', 'DG'], ['23', 'DC']],
 [['3', 'DC'], ['22', 'DG']],
 [['4', 'DG'], ['21', 'DC']],
 [['5', 'DA'], ['20', 'DT']],
 [['6', 'DA'], ['19', 'DT']],
 [['7', 'DT'], ['18', 'DA']],
 [['8', 'DT'], ['17', 'DA']],
 [['9', 'DC'], ['16', 'DG']],
 [['10', 'DG'], ['15', 'DC']],
 [['11', 'DC'], ['14', 'DG']],
 [['12', 'DG'], ['13', 'DC']]]

Возьмём пары по середине (6 и 7) и будем считать их 1 и 2 парой. Пара 0 будет как пара 1, но её используем, чтобы найти матрицу перехода.

In [520]:
cmd.create('0', 'chain A and resi 6 or chain B and resi 19')
cmd.create('1', '0')
cmd.create('2', 'chain A and resi 7 or chain B and resi 18')
In [521]:
cmd.pair_fit('0', '2')
trans_matrix = cmd.get_object_matrix('0')
cmd.delete('0')

Создадим оставшиеся пары нуклеотидов простым копированием. Добавим свойство resi к каждой паре.

In [522]:
for i in range(3, length + 1):
    cmd.create(str(i), str(i - 1))
    cmd.transform_selection(str(i), trans_matrix)
for i in range(length):
    cmd.alter(str(i + 1), f'resi={i + 1}')

Удалим модель ДНК и создадим новую из наших пар.

In [530]:
cmd.delete('355d')
cmd.create('dna', 'all')

Удалим пары и образуем связи сахаро-фосфатного остова.

In [534]:
for i in range(length):
    cmd.delete(str(i + 1))
    cmd.bond("/dna/A/A/*`{i + 1}/O3'", "/dna/A/A/*`{i + 2}/P")
    cmd.bond("/dna/B/B/*`{i + 2}/O3'", "/dna/A/A/*`{i + 1}/P")

PyMol не хочет красиво отрисовывать картонки, поэтому приходится сохранять модель в pdb-файл и снова загружать.

In [545]:
cmd.save('bdna.pdb')
cmd.reinitialize()
cmd.do('bg_color white; set antialias, 2; set ray_trace_mode, 3;')
cmd.load('bdna.pdb')
In [559]:
cmd.hide('all')
cmd.set('cartoon_ring_mode', 2)
cmd.show('cartoon')
cmd.center('all')
cmd.orient('all')
cmd.zoom('all')
In [560]:
cmd.do('png bdna_1.png, 1920, 1080')
In [561]:
cmd.rotate('y', 45)
In [562]:
cmd.zoom('bdna and resi 1-10')
In [565]:
cmd.do('png bdna_2.png, 1920, 1080')
In [566]:
display(Image('bdna_1.png', width=720),
        Image('bdna_2.png', width=720))

Получилось красиво.