import sys
stdout = sys.stdout

import __main__
__main__.pymol_argv = [ 'pymol']

import pymol
pymol.finish_launching()
from pymol import cmd,stored

Qt not available (pymol.Qt), using GLUT/Tk interface

sys.stdout = stdout

from IPython.display import Image, HTML, Video

"Воспроизведите работу этого скрипта на свой вкус"

cmd.reinitialize()

# скачиваем структуру кальмодулина, зумимся на начало
# не показываю как lines, потому что потом всё равно меняем на cartoon
# создаём видео длиной 1000 кадров

cmd.do('''
fetch 1cll, async=0
bg_color white
mset 1 x1600
mview store''')

cmd.remove('resi -3') # содержат по одному атому без координат (?), портят интерполяцию потом

# добываем номера всех присутсвующих в структуре остатков
stored.r = [] 
cmd.iterate('1cll and n. CA','stored.r.append(int(resi))')

149

import numpy as np

# создаём градиент от голубого (rgb 0.5, 0.75, 1) на N-конце к зелёному (rgb 0.5, 0.75, 0.5) на C-конце
length = len(stored.r)
colors = np.linspace(1, 0.5, length)
for i, resi in enumerate(stored.r):
    cmd.set_color(f'col{i}', [0.5, 0.75, colors[i]])
    cmd.set('cartoon_color',f'col{i}' ,f'resi {resi}')
cmd.show_as('cartoon','all')

# проходим вдоль остова белка, последовательно зумясь на его аминокислоты
for i, resi in enumerate(stored.r):
    cmd.frame(56 + (10*i))
    cmd.zoom( f'n. CA and resi {resi}')
    cmd.mview('store')

cmd.stereo('crosseye') # для более захватывающего опыта

# экспортируем видео и за кадром превращаем в mp3
cmd.movie.produce('movie0.mpeg', quality=100)

using encoder "mpeg_encode"

Итого получаем 3d-фильм о белковом аналоге американских горок:

Video(url='https://kodomo.fbb.msu.ru/~ta.rom/term8/pr1/movie0.mp4', height=300)

Мутация и видео

# mutate

cmd.reinitialize()
cmd.fetch('1LMP')

cmd.wizard("mutagenesis")
cmd.do("refresh_wizard")
cmd.get_wizard().set_mode("PHE")
cmd.get_wizard().do_select("1LMP///101/")
cmd.get_wizard().apply()

cmd.save('1LMP_D101F.pdb')

# movie

cmd.do('''
# setup PyMOL for movies
reinitialize
set movie_auto_interpolate, off
set matrix_mode, 1
set movie_panel, 1
set scene_buttons, 1
set cache_frames, 1
config_mouse three_button_motions, 1

# download the complex and set it up
fetch 1LMP, orig
load 1LMP_D101F.pdb, mut
dist d, /orig/A/A/ASP`101/OD2, /orig/B/B/NAG`3/O6
remove solvent
select prot, polymer.protein
select myres, resi 101
select lig, chain B
color palegreen, mut and prot
color forest, mut and lig
color gray90, orig and prot
color gray50, orig and lig
color palecyan, d
hide label, d
bg_color white
color atomic, not element C
show sticks, myres
deselect


mset 1 x300

   
# move the mut to the left
translate [-80,0,0], object=mut

# overview of the scene
frame 1
mview store
mview store, object=orig
mview store, object=mut


frame 80
translate [+80,0,0], object=mut
mview store
mview store, object=orig
mview store, object=mut


frame 100
mview store
mview store, object=orig
mview store, object=mut


frame 180
set_view (\
     0.655596912,   -0.007433767,   -0.755073845,\
     0.366541237,    0.877375185,    0.309612960,\
     0.660181880,   -0.479747117,    0.577928841,\
     0.000004765,    0.000007735,  -35.530853271,\
     9.126789093,   56.358558655,   29.965265274,\
    27.702697754,   43.356124878,  -20.000000000 )
mview store
mview store, object=orig
mview store, object=mut
mview store, object=d


frame 240
mview store
mview store, object=orig
mview store, object=mut
mview store, object=d


mview interpolate, object=mut
mview reinterpolate

''')

cmd.movie.produce('movie.mpeg', quality=100)

Video(url='https://kodomo.fbb.msu.ru/~ta.rom/term8/pr1/movie.mp4', height=300)
# перевела в mp4, потому что mpeg не воспроизводится в html
# при этом потерялось качество(

флуоресцентная метка TAMRA

cmd.reinitialize()
cmd.fetch('1LMP', 'prot')
cmd.fetch('CID_2762604', 'fluor')
cmd.remove('solvent')

# прикрепим сложноэфирнной связью на какую-нибулдь гидроксильную группу
cmd.show('sticks', '/prot/A/A/SER`81/')

# склеим по кислородам, удалим лишний
cmd.remove('/prot/A/A/SER`81/OG')
cmd.remove('n. H*')
cmd.fuse('fluor and id 4', '/prot/A/A/SER`81/CB') # в метке по id, потому что имена у всех одинаковые
cmd.delete('fluor')

# повернём, чтобы не втыкалась в белок
cmd.edit('prot///UNK`0/O04', '/prot/A/A/SER`81/CB')
cmd.torsion(-30)
cmd.edit('/prot///UNK`0/C', 'prot///UNK`0/O04')
cmd.torsion(-180)

cmd.unpick()

cmd.orient('/prot/A/A/SER`81')
cmd.zoom('resn UNK extend 5')

cmd.bg_color('white')
cmd.png('fluor.png', height=500, ray=1)

1

Image("fluor.png", width=500)

Poly-Ala

%%time

# готовое решение раз
cmd.reinitialize()
cmd.fab('A'*100, ss=1)
cmd.save('polyAla_fab.pdb')

CPU times: total: 2.86 s
Wall time: 2.89 s

%%time

# готовое решение два
cmd.reinitialize()
cmd.set('secondary_structure', 1)
for _ in range(100): cmd._alt('a')
cmd.save('polyAla_alt.pdb')

CPU times: total: 1.8 s
Wall time: 1.82 s

%%time

# готовое решение три
cmd.reinitialize()
cmd.editor.attach_amino_acid(None, 'ALA', ss=1, object='polyAla')
for _ in range(99): cmd.editor.attach_amino_acid('pk1', 'ALA', ss=1)
cmd.save('polyAla_aaa.pdb')

CPU times: total: 3.17 s
Wall time: 3.12 s

С attach_fragment очень муторно возиться, потому что он считает добавленный фрагмент частью исходного остатка, нумерует и переименовывает соответсвенно. С чистым attach ещё больше возни. Если строить "с нуля", то проще тогда тем же fuse:

%%time

# ручное решение через fuse (по мотивам attach_amino_acid)

cmd.reinitialize()
cmd.fragment('ALA', 'polyAla', origin=0)
cmd.alter('polyAla', f"resi=1")


phi=-57.0
psi=-47.0

for res in range(1, 100):

    cmd.fragment('ALA', 'new_res', origin=0)
    cmd.alter('new_res', f"resi={res+1}")
    cmd.set_geometry('new_res and n. N', 3, 3) # make nitrogen planar
    cmd.fuse('new_res and n. N', f'{res}/C', 2) # mode=2: adopt segi/chain
    cmd.delete('new_res')
    

    # omega
    cmd.set_dihedral(f'{res}/O',
                    f'{res}/C',
                    f'{res+1}/N',
                    f'{res+1}/H', 
                    180)
    
    # phi
    cmd.set_dihedral(f'{res}/C',
                    f'{res+1}/N',
                    f'{res+1}/CA', 
                    f'{res+1}/C',
                    phi)

    # psi
    cmd.set_dihedral(f'{res}/N',
                    f'{res}/CA', 
                    f'{res}/C',
                    f'{res+1}/N',
                    psi)
cmd.save('polyAla_man.pdb')

CPU times: total: 3.25 s
Wall time: 3.07 s

# визуализируем итог

cmd.reinitialize()

cmd.load('polyAla_fab.pdb')
cmd.load('polyAla_alt.pdb')
cmd.load('polyAla_aaa.pdb')
cmd.load('polyAla_man.pdb')


cmd.alignto('polyAla_fab')
cmd.orient()

cmd.translate([0, 20, 0], 'polyAla_alt')
cmd.translate([0, 40, 0], 'polyAla_aaa')
cmd.translate([0, 60, 0], 'polyAla_man')

cmd.zoom()

cmd.set('cartoon_side_chain_helper', 1)
cmd.show('sticks')

cmd.bg_color('white')
cmd.png('polyAla.png', height=500, ray=1)

Image('polyAla.png', width=600)

В ручном воспроизведении (верхнее=жёлтое), похоже, не учла что-то с концами, в остальном работает.

Из любопытства можно сравнить времена работы. "Интерактивный" вариант с альтом (который работает при alt+A) получается самым быстрым.

Poly-A

Немного читерства: воспользуемся работой, сделанной до нас другими людьми и извлечём нужную матрицу из идеализированной ДНК, построенной 3DNA и одиночных выравненных в референсную систему отсчёта пар нуклеотидов оттуда же.

cmd.reinitialize()
cmd.load('gatc-b.pdb')  # остался с третьего семестра (fiber -b последовательности gatc пять раз)
cmd.load('a.pdb')  # fiber -b одиночного a
cmd.align('a', 'gatc-b and resi 2')
transformation = cmd.get_object_matrix('a', state=1)

transformation

(0.8090222477912903,
 0.5877780318260193,
 -1.01000059657963e-05,
 0.0001397366397384303,
 -0.5877780318260193,
 0.8090222477912903,
 1.589298881299328e-05,
 8.710888340379697e-05,
 1.7512678823550232e-05,
 -6.921220574440667e-06,
 1.0,
 -3.3749778410082887,
 0.0,
 0.0,
 0.0,
 1.0)

transformation = (0.8090222477912903,
 0.5877780318260193,
 -1.01000059657963e-05,
 0.0001397366397384303,
 -0.5877780318260193,
 0.8090222477912903,
 1.589298881299328e-05,
 8.710888340379697e-05,
 1.7512678823550232e-05,
 -6.921220574440667e-06,
 1.0,
 -3.3749778410082887,
 0.0,
 0.0,
 0.0,
 1.0)

def gen_dna(seq):
    cmd.reinitialize()
    cmd.create('dna', None)
    
    for resi, res in enumerate(seq):
        cmd.load(f'{res}.pdb', 'new_bp')
        
        for _ in range(resi):
            cmd.transform_selection('new_bp', transformation)
            # будь это нормальная матрица трансформации, можно было бы воспользоваться
            # матричным умножением для ускорения, но это не она, а особый pymol'овский объект        

        cmd.alter('new_bp and chain A', f'resi={resi+1}')
        cmd.alter('new_bp and chain B', f'resi={len(seq)-resi}')
        cmd.create('dna', 'dna or new_bp')
        cmd.delete('new_bp')

        if resi>0:
            cmd.bond(f"A/{resi}/O3'", f"A/{resi+1}/P")
            cmd.bond(f"B/{resi}/P", f"B/{resi+1}/O3'")

seq = 'gatc' * 5
gen_dna(seq)
cmd.save('dna20.pdb')

cmd.reinitialize()
cmd.bg_color('white')
cmd.load('gatc-b.pdb', 'ref')
cmd.color('gray', 'ref')
cmd.load('dna20.pdb')
cmd.orient()

cmd.show('sticks')
cmd.hide('cartoon')

cmd.png('ref_dna.png', height=300, ray=1)

1

Image('ref_dna.png')

cmd.rms_cur('dna20', 'ref')

0.0018551647663116455

Референстная и построенная структура сопадают практически идально.

Теперь к заданию:

seq = 'a' * 100
gen_dna(seq)
cmd.save('dna100.pdb')

cmd.reinitialize()
cmd.bg_color('white')
cmd.load('dna100.pdb')
cmd.orient()

cmd.png('dna100.png', height=500, width=1600, ray=1)

1

Image('dna100.png', width=800)

ДНК кажется кривой, но это просто особенности перспективы, проявляющиеся из-за большой длины. Если перспективу выключить, всё становится в порядке:

cmd.set('orthoscopic', 'on')
cmd.png('dna100_ortho.png', height=500, width=1600, ray=1)

1

Image('dna100_ortho.png', width=800)

Blender

Image('str2.png', width=600)