Содержание
Повествование предполагает, что читатель уже знаком с pymol с точки зрения пользователя.
Программирование для Pymol
Если рассуждать теоретически, то pymol – это прекрасное средство. Однако, всякий раз, когда я начинаю изучать его поближе – например, с тем, чтобы рассказать вам о том, как для него программировать, – я всякий раз прихожу к мнению, что pymol – замечательный пример того, как не надо программировать.
И всё же, многие вещи для него сделать довольно просто.
Итак.
Запуск скриптов
Первый способ запуска скриптов в pymol – командой @. Такие скрипты содержат в себе просто последовательность команд pymol так же, как если бы мы их писали в его командной строке.
Полагаю, с этим видом скриптов вы уже знакомы.
Эти скрипты, как правило, довольно понятные, но в них нельзя ни использовать какую-нибудь информацию из структуры, чтобы принять решение, что изображать, ни использовать более компактную запись за счёт циклов или функций.
Соответственно, чтобы дать пользователю pymol в руки более адекватные средства, в pymol есть второй вид скриптов – они запускаются командой run. Эти скрипты пишутся на питоне и дальше пойдёт речь о них.
Модуль cmd
В модуле cmd находятся функции, которые, реализуют все команды командной строки pymol. Фактически, всё взаимодействие с pymol из скриптов происходит через вызов функций из этого модуля.
Пример:
Эти строки питона эквивалентны pymol-команде:
color red, all
Нелирическое отступление: from ... import ...
До сих пор про работу с модулями я рассказывал только про одну форму команды import: import имя_модуля. И хорошо бы так всегда и оставалось, только не всегда так удобно писать. Теперь мы наткнулись на ещё одну форму from имя_модуля import что_нибудь.
Команда from pymol import cmd эквивалентна на самом деле трём командам подряд:
То есть мы загружаем модуль, сразу берём что-то из его внутренностей, а имя самого модуля выбрасываем. Понятно, для чего это полезно: чтобы обращение внутрь модуля было более коротким.
Этим способом не стоит злоупотреблять: если вы пользуетесь конструкциями from, то глядя на какой-нибудь вызов функции в середине программы, сложнее становится угадать, откуда она пришла и как про неё добыть хелпы.
cmd.iterate & cmd.alter
Зададимся задачей, чтобы на её примере продемонстрировать возможности pymol. Здача будет такая: мы хотим в pymol выделить все точки излома полипептидной цепи.
DSSP определяет точку излома полипептидной цепи следующим образом: построим вектора направления полипептидной цепи – это вектора, соединяющие центры Cα-атомов, стоящих в цепи через один; если между двумя смежными векторами направления острый угол больше 70 градусов, то в их смежной точке есть излом цепи.
Чтобы решить эту задачу, для начала нужно научиться получать из pymol данные о структуре. Это можно делать несколькими способами, основные из которых можно назвать по реализующим их функциям: cmd.iterate, cmd.get_model и cmd.identify.
Начнём с cmd.iterate.
В pymol есть семейство из четырёх очень похожих команд: iterate, alter, iterate_state и alter_state.
iterate и alter используются одинаковым образом: они получают два аргумента, выделение и строку питонского кода, и для каждого атома выделения выполняют строку питонского кода, при этом подсовывая в неё в виде переменных свойства очередного атома. Разница между ними в том, что alter дополнительно, после исполнения питонской строки для атома, собирает данные из этих же переменных и возвращает их во внутреннее состояние структуры в pymol. По сути это значит, что при помощи alter можно менять структуру, а при помощи iterate нельзя (но зато она чуть-чуть быстрее работает1, яснее выражает намерение "только посмотреть" и защищает от ошибки, если вы случайно всё-таки ещё и "потрогали"). Маленькое замечание: даже после того, как alter изменил структуру, pymol не отобразит изменения; чтобы изменения вступили в силу, после alter нужно сказать rebuild.
Команды iterate_state и alter_state полностью аналогичны, но при этом позволяют выбирать, на каком "состоянии" команды работают. (iterate и alter работают на всех "состояниях" pymol). Что такое "состояние", я вам внятно объяснить не могу – вроде бы предлагается ими представлять разные конформации структуры. (Почему состояния не сделали одним из критериев алгебры выделений – одному Аллаху только и известно).
Можно сказать, например, в pymol:
iterate n. ca, print([x,y,z])
И pymol распечатает в командном окне координаты всех Cα атомов.
Самого по себе iterate нам недостаточно, чтобы собрать данные из структуры. Нам нужно ещё какое-нибудь место, через которое можно было бы передать данные, полученные строками питонского кода в iterate наружу. Поясню. Если мы сделаем скрипт bend.py с текстом:
или даже разумнее:
то мы ничего, кроме сообщения об ошибке, что переменной coords не существует, не получим. Код внутри iterate исполняется в специальном замкнутом пространстве, откуда переменные нашего bend.py не видны. Чтобы это обойти, в pymol сделали специальную штуку по имени stored. Изнутри pymol оно выглядит как просто питонский объект, в который можно дописывать по мере желания свои поля. Из питонских скриптов, которые мы пишем, оно выглядит как модуль. Соответственно, пользоваться этим так:
Вот после таких слов мы увидим, наконец, список координат Cα-атомов.
Нелирическое отступление: Pymol, Chempy
Чтобы искать точки излома, нам нужно научиться делать всяческие операции над векторами в pymol: считать разность, считать углы между векторами. Чтобы их научиться делать, полезно иметь немножко представления, как устроен pymol.
А устроен он по большому счёту так:
Есть для питона маленькая и немогучая (сравнительно) библиотека для работы со всякой молекулярной биологией, зовётся chempy. (Пишут её ровно те же ребята, кто делают pymol). В этой библиотечке лежат всякие средства для работы со всем трёхмерным: там есть чтение-запись форматов PDB и MDL, там есть классы Atom (в котором есть всё про атом и даже его координаты) и Bond (в котором хотя бы есть номера атомов, которые pymol считает связанными), там есть библиотека функций для работы 3-мерными (и только 3-мерными) векторами по имени cpv (и не спрашивайте меня, как это имя расшифровывается, я не знаю), и много чего ещё. Вектора в этой библиотеке представляются не иначе как список или тупль из трёх чисел. (И это я считаю определённым идиотизмом, потому, что это заставляет нас писать гораздо более некрасивый код, чем можно было бы в питоне). Эта библиотека является частью pymol.
Рядом с ней, как-то с ней взаимодействующее есть ядро pymol, написанное зачем-то на славном языке C. Когда из питона мы смотрим на части кода, написанные на си, мы не видим про них, как правило, документации, и ими, как правило, довольно неудобно пользоваться. Соответственно, выглядит это ядро pymol большим и непостижимым. На мой скромный вкус, это пример нарушения авторами pymol первого правила оптимизации, правильнее было писать весь pymol на питоне и только проседающие по производительности места (т.н. "узкие горлышки") переносить в си.
И вокруг этого сишного ядра написана обёртка на питоне – именно эту обёртку мы в первую очередь и видим.
find_bends
Всё это я говорил к тому, где же именно нам искать способы борьбы с векторами.
Кажется осмысленным отделить часть, которая собственно ищет точки излома от всяческой работы с pymol, поэтому мы опишем функцию, которая получает на вход список координат Cα-атомов и ищет среди них точки излома...
Вооружившись новыми знаниями, пишем:
И тут нам приходится задуматься, а как же мы будем из этой функции возвращать, в каких атомах произошёл излом. Для этого нужно решить, как мы будем делать выделение из интересующих нас атомов. Я предлагаю для этого использовать поле index, которое, покуда структура не изменяется, однозначно идентифицирует атом.
Соответственно, теперь считаем, что в списке atoms находятся не просто координаты, а пары из индекса атома и его координат. И для найденных точек излома мы будем возвращать индексы атомов, в которых излом происходит:
1 import chempy.cpv
2 import math
3
4 def find_bends(atoms):
5 bends = []
6 for i in range(len(atoms) - 4):
7 idx1, ca1 = atoms[i]
8 idx2, ca2 = atoms[i+2]
9 idx3, ca3 = atoms[i+4]
10 dir1 = chempy.cpv.sub(ca2, ca1)
11 dir2 = chempy.cpv.sub(ca3, ca2)
12 angle = chempy.cpv.angle(dir1, dir2)
13 if angle > math.radians(70):
14 bends.append(idx2)
15 return bends
cmd.extend
Теперь дело осталось за малым: собрать из получившегося списка выделение и отправить его в pymol:
from pymol import cmd, stored
... # та часть, которую мы описали выше
stored.ca = []
cmd.iterate("n. ca", "stored.ca.append((index, (x,y,z)))")
bends = find_bends(stored.ca)
selections = []
for bend in bends:
selections.append("(idx. %s)" % bend)
selection = " or ".join(selections)
cmd.select("bends", selection)Теперь мы можем кучу раз говорить run bends.py и для каждого излома и, скорее всего, разрыва цепи мы получим по атому в выделении под именем bends.
В структуре 1m5x, которую я везде тыкаю в качестве примера, нуклеазу, насколько я понимаю, дезактивировали, присоединив к ней несколько йонов кальция, и угадайте, как соответствующие атомы в PDB называются! Правильно, ca! И убедить наш скрипт их не выделать довольно трудно.
Поэтому лучше было бы сделать, чтобы наш скрипт получал в качестве аргумента выделение, на котором он действует. Только вот pymol не разрешает передавать скриптам параметры (вместо этого он позволяет запускать одной командой run несколько скриптов). Зато pymol разрешает доопределить новую команду командной строки pymol. (Ужасная тавтология, да). Для этого нам нужно сложить всё, что мы назовём одной командой, в функцию, и сказать pymol, чтобы он проассоциировал какую-нибудь команду с этой функцией, делается это так:
У cmd.extend два аргумента. Первый – строка, под каким именем наша функция будет видна в командной строке pymol. Второй – функция (без кавычек, без скобочек, только имя функции), под каким именем мы её описали в питонском модуле.
В итоге получаем:
1 from pymol import cmd, stored
2 from chempy import cpv
3
4 def find_bends(atoms):
5 bends = []
6 for i in range(len(atoms) - 4):
7 idx1, ca1 = atoms[i]
8 idx2, ca2 = atoms[i+2]
9 idx3, ca3 = atoms[i+4]
10 dir1 = chempy.cpv.sub(ca2, ca1)
11 dir2 = chempy.cpv.sub(ca3, ca2)
12 angle = chempy.cpv.angle(dir1, dir2)
13 if angle > math.radians(70):
14 bends.append(idx2)
15 return bends
16
17 def bend(selection):
18 # Выбираем подмножество selection, содержащее только CA-атомы.
19 selection = "(%s) and n. ca" % selection
20 # Ищем изломы.
21 stored.ca = []
22 cmd.iterate(selection, "stored.ca.append((index, (x,y,z)))")
23 bends = find_bends(stored.ca)
24 # Строим новое выделение.
25 selections = []
26 for bend in bends:
27 selections.append("(idx. %s)" % bend)
28 selection = " or ".join(selections)
29 cmd.select("bends", selection)
30
31 cmd.extend("bend", bend)
Не шибко изящно, но работает2. Теперь в pymol мы можем сказать run bend.py и bend all или bend chain b или bend chain b and not resn ca – и получить выделение с нашими точками излома.
Переменные для iterate
Тут такой вопрос: а откуда узнать, какие же именно переменные pymol делает доступными для питонской строки в cmd.iterate? Мы пока увидели только index, x, y, z, что мало.
Ответ: это не задокументировано.
Вот такая прекрасная программа pymol.
Но это компенсируется тем, какой прекрасный язык питон.
Вспоминаем функцию vars из рассказа про объекты и классы (а точнее про словари и исключения): она возвращает словарь, которым реализованы поля объекта. Если функцию vars() вызвать без аргументов, то она вернёт словарь, которым реализованы локальные переменные, с которыми мы можем в питоне работать. То есть ровно то, что нам хочется узнать.
Поэтому мы говорим pymol:
PyMOL>iterate id 1, print sorted(list(vars())) ['ID', 'alt', 'b', 'cartoon', 'chain', 'color', 'elec_radius', 'elem', 'flags', 'formal_charge', 'index', 'label', 'model', 'name', 'numeric_type', 'partial_charge', 'q', 'rank', 'resi', 'resn', 'resv', 'segi', 'ss', 'state', 'text_type', 'type', 'vdw', 'x', 'y', 'z'] Iterate: iterated over 1 atoms.
и получаем список имён переменных, которые нам доступны.
(Конструкции, в которых много скобочек, нужно начинать читать изнутри: vars() – словарь локальных переменных, т.е. его ключ – имя переменной, значение – значение переменной; list(словарь) – когда словарь оказывается там, где питон ожидает список, словарь воспринимается как список ключей; таким образом, list(vars()) – список имён локальных переменных; sorted(list(vars())) – отсортированный в алфавитном порядке список имён локальных переменных – так читать удобнее; iterate id 1, ... – применить ... к атому с id == 1 – я надеюсь, с id == 1 может быть только один атом, нам важно сделать какую-нибудь такую команду iterate, которая делает ровно один шаг, чтобы не забивать себе экран всякой мурой).
Можно сократить до: iterate id 1, print sorted(vars()) – list здесь на самом деле не нужен – или вовсе до iterate id 1, print vars(), чтобы посмотреть и значения переменных тоже.
Смысл этих переменных местами очевиден, а местами можно догадаться о нём, сравнивая имя переменной с чем-нибудь похожим в списке Property Selectors.
cmd.get_model
iterate и alter – довольно шаманские команды. (Например, я видел версию pymol, которая питонским строкам внутри iterate не дают координаты атома ни в каком виде).
Можно попробовать сделать то же самое несколько красивее. Вспоминаем, что pymol написан в тесной связке с библиотекой chempy, а в ней есть довольно внятное объектное представление атомов. Можно попросить pymol вернуть нам часть структуры в виде объектов chempy. Делается это командой get_model.
cmd.get_model(выделение) возвращает для выделения объект "модели" – это объект, который содержит четыре поля:
atom – список атомов (объектов класса chempy.Atom)
bond – список валентных связей, которые рисует pymol (объекты класса chempy.Bond).
molecule – какая-то непонятная фигня
index – по идее, это должен быть словарь, ключи которого – целое число (поле index, которое мы как раз использовали выше), а значения – атомы, но в том виде, как модель возвращает cmd.get_model, это поле всегда остаётся со значением None.
Разумеется, в нашем случае (как и во всех, какие были у меня до сих пор), нас интересует только список атомов.
Класс chempy.Atom содержит примерно3 те же поля, что и локальные переменные в iterate – самое значимое отличие в том, что координаты атома представлены полем coords, в котором лежит список из трёх чисел – в точности вектор с точки зрения chempy.
Как вы можете догадаться из такого описания, chempy слегка недодокументирован. (Иначе бы я обошёлся вместо этого описания словами RTFM – Read That... эээ... Fascinating Manual).
Соответственно, поиск информации про него аналогичен томе же про iterate. Здесь мы пользуемся тем, что встроенная команда print в pymol на самом деле получает аргументом питонский код (а не просто строку) точно так же, как iterate, и тем позволяет нам весьма безболезненно подглядывать в то, как pymol устроен:
PyMol> print help(chempy.Atom)
... -- большая куча весьма бесполезного
PyMol> print sorted(vars(cmd.get_model("all").atom[0]))
['b', 'chain', 'coord', 'elec_radius', 'flags', 'formal_charge', 'hetatm', 'id', 'index', 'name', 'numeric_type', 'partial_charge', 'q', 'resi', 'resi_number', 'resn', 'segi', 'ss', 'stereo', 'symbol', 'u_aniso', 'vdw']get_model всякий раз, когда мы его вызываем, создаёт все объекты, котороые он возвращает, заново, поэтому сколько бы мы ни меняли модель, которую мы таким образом получили, мы ничего не испортим (но и повлиять на pymol'ское мировоззрение не сможем).
Теперь, когда мы всё про него знаем, мы можем попытаться им воспользоваться. Мы замечаем, что атом у нас сам по себе сразу несёт в себе и свою идентификацию, и координаты, поэтому нам не нужно передавать в find_bends пары атом, координата, а достаточно передавать туда атомы и возврващать оттуда тоже атомы.
1 from pymol import cmd
2 from chempy import cpv
3 import math
4
5 def find_bends(atoms):
6 bends = []
7 for i in range(len(atoms) - 4):
8 dir1 = cpv.sub(atoms[i+2].coord, atoms[i].coord)
9 dir2 = cpv.sub(atoms[i+4].coord, atoms[i+2].coord)
10 angle = cpv.angle(dir1, dir2)
11 if angle > math.radians(70):
12 bends.append(atoms[i+2])
13 return bends
14
15 defind bend(selection):
16 selection = "(%s) and n. ca" % selection
17 atoms = cmd.get_model(selection).atom
18 bends = find_bends(atoms)
19
20 selections = []
21 for atom in bends:
22 selections.append("(idx. %s)" % atom.index)
23 selection = " or ".join(selections)
24 cmd.select("bend", selection)
25
26 cmd.extend("bend", bend)
get_*
Есть и ещё один способ добиться того же самого. В pymol есть набор полезных встроенных команд, которые считают разные штуки про структуру. Все они называются get_что_нибудь. Мне из них наиболее интересны три:
get_distance sel1, sel1 – считает расстояние между двумя атомами
get_angle sel1, sel2, sel3 – считает угол по трём атомам
get_dihedral sel1, sel2, sel3, sel4 – считает торсионный угол по четырём атомам
Соответственно, чтобы суметь этими командами воспользоваться, нам нужно получить список индексов атомов в выделении. Это делается командой: identify. Правда возвращает она не список значений index, а список значений id.
Собирая вместе эти кусочки знаний, получаем:
1 from pymol import cmd
2 import math
3
4 def bend(selection):
5 selection = "(%s) and n. ca" % selection
6 bends = []
7 ids = cmd.identify(selection)
8 for i in range(len(ids) - 4):
9 ca1 = "id %s" % ids[i]
10 ca2 = "id %s" % ids[i + 2]
11 ca3 = "id %s" % ids[i + 4]
12 angle = cmd.get_angle(ca1, ca2, ca3)
13 if angle > math.radians(70):
14 bends.append(ca2)
15 selection = " or ".join(bends)
16 cmd.select("bends", selection)
17
18 cmd.extend("bend", bend)
При всей красивости третьего варианта, для программы, которой я буду пользоваться больше одного раза, и которая будет достаточно сложной, я бы предпочёл второй подход: самая ключевая логика в нём завязана на chempy, а не на pymol – следовательно, можно было бы чуть-чуть постараться и сделать, чтобы этими функциями можно было пользоваться и в других случаях. Впрочем, может быть, это просто ещё один вид преждевременной оптимизации: оптимизация для повторного использования кода.
index vs. id
Про многие свойства атома можно понять из документации, что они значат, либо понять по их названию, что они не интересны. А вот index и id – это явно полезные вещи, которые так и тянет путать.
Поэтому проговариваю отдельно и внятно.
id – номер атома согласно PDB.
index – это номер атома во внутре-pymol-ской организации. Он не совпадает с id и назначается по какому-то странному принципу.
Если в структуре удалить какой-нибудь атом, то index многих атомов сдвинутся, а id останутся неизменными (если мы удаляем строчку в PDB, то поле с номером атома в каждой строке всё равно остаётся тем же).
Покуда структура остаётся достаточно неизменной (то есть не добавляются/удаляются новые атомы или связи), пользоваться можно что id, что index – чем проще.
alter
Напоследок вернёмся к нашему примеру.
На самом деле, конечно же, правильно было бы делать, чтобы наша программа не создавала выделение с точками излома, а помечала информацию о вторичной структуре. Для этого в pymol имеется поле ss (secondary structure).
Единственный способ как-то повлиять на внутреннее представление структуры в pymol – используя команду alter.
Соответственно:
1 from pymol import cmd
2 import math
3
4 def bends(selection):
5 selection = "(%s) and n. ca" % selection
6 bends = []
7 ids = cmd.identify(selection)
8 for i in range(len(ids) - 4):
9 ca1 = "id %s" % ids[i]
10 ca2 = "id %s" % ids[i + 2]
11 ca3 = "id %s" % ids[i + 4]
12 angle = cmd.get_angle(ca1, ca2, ca3)
13 if angle > math.radians(70):
14 bends.append(ca2)
15 selection = " or ".join(bends)
16 cmd.alter(selection, "ss = 'B'")
17 cmd.select("bend", "ss B")
18
19 cmd.extend("mark_bends", bends)
Команда alter велика и непостижима, поэтому по первому проверяйте результаты её выполнения руками. Например, когда мы меняем поле ss, она у нас за спиной меняет регистр букв с маленьких на большие. (Поэтому в этом примере я сразу сделал букву 'B' заглавной — чтобы наша программа по принципу наименьшего удивления).
Документация
Теперь о самом главном.
Так как pymol написан замечательными людьми, источников документации о том, как для него программировать, приходится использовать довольно много:
- К pymol прилагается хорошая и подробная документация. Наверное. Я её никогда не видел в глаза потому, что авторы продают её исключительно за деньги.
Зато вокруг pymol сложилось большое сообщество пользователей, которые совместными усилиями выдирают из него куски полезной информации. Все такие усилия собираются вместе на сайте http://pymolwiki.org – набор самых полезных ссылок на документацию на нём на мой вкус я сделал на странице занятия.
Про всё содержимое модуля cmd, т.е. про все команды командной строки pymol (пардон), имеется прекрасная встроенная документация, доступная из командной строки командой help
К сожалению, эта команда закрывает нам возможность использовать встроенный питонский help() для изучения внутренностей pymol. Но это можно обойти: команда print командной строки pymol исполняет питонский код и пишет его результат в консоль – то, что нам нужно. Соответственно: print help(chempy), или print help(cmd) или print help(pymol).
Задокументировано в pymol далеко не всё из тех внутренностей, которыми хотелось бы пользоваться для программирования. Поэтому приодится пользоваться просто возможностями питона по изучению самого себя4, два самых важных примера тому я вам показывал, если вы поняли, откуда взялись они, вы сможете сочинить и свои подобные: print sorted(vars(cmd.get_model("all").atom[0])) и iterate id 1, print sorted(vars())
Полезные мелочи
И напоследок, чуточку полезных команд pymol:
cd – перейти в другую директорию
ls – посмотреть содержимое директории
pwd – какая директория является текущей?
Очень полезно, чтобы не писать каждый раз полный путь к каждому скрипту и каждому pdb-файлу, с которыми работаем.
По правде сказать, я вам рассказал даже не все подходы к тому, как с pymol можно работать. Я рассказал только то из этого, что на мой вкус достойно внимания. Авторы заложили в pymol довольно много весьма дурацких конструкций сверх этого.