Учебная страница курса биоинформатики,
год поступления 2011

Семестры Студенты Преподаватели

• 2024
• 2023
• 2022
• 2021
• 2020
• 2019
• 2018
• 2017

• ...

Работа с 3D-структурами в Bio python

Мотивация

Если вы много работаете с 3D-структурами, полезно иметь какое-то средство для того, чтобы автоматизировать работу с ними.

Если вы пользуетесь pymol, в нём есть по-своему замечательная (но по-своему и ужасная) библиотека для работы со структурами, и он позволяет писать на питоне скрипты и исполнять их. К сожалению, пакет pymol для питона не отчуждаем от самой программы pymol, поэтому вы будете при этом всегда вынуждены общаться с их графической оболочкой, да и запуск программы на много файлов будет усложнён.

С другой стороны, если для вашей задачи вам требуется обработать 20, 100 или 1000 PDB-файлов, и пусть даже сделать с ними что-нибудь совсем примитивное: например, повернуть все молекулы одинаковым образом или убрать везде боковые цепи или перенумеровать все остатки таким образом, чтобы ваша программа для молекулярного моделирования смогла с ними справиться (многие такие программы выдвигают иногда довольно странные требования к PBD-файлам, которые выполняются на большинстве файлов в банке PDB – но далеко не на всех), убрать воду, достроить водороды, посчитать вес молекулы, посчитать статистику расстояний интересного вам типа контакта, найти расстояния между остатками активного центра, отмеченными в выравниваниии, вырезать только нужный фрагмент PDB, найти остатки, которые контактируют с лигандом (или наоборот, найти, какие лиганды контактируют с заданным сайтом) – в любом этом и множестве подобных случаев удобнее написать самостоятельный маленький скрипт, который будет это делать. (Надо надеяться, запускать скрипт из шелла так, чтобы он применился к множеству нужных вам файлов, вы уже умеете). И вот для этих целей идеально подходит пакет Bio.PDB из биопитона.

Биопитон достаточно выразителен и для того, чтобы делать свои большие проекты: например, для того, чтобы сделать свой предсказатель вторичной структуры, свою программу для утряски или подбора мутаций и т.п. (Решение таких задач безусловно потребует от вас более глубокого вникания в возможности и устройство этого пакета).

Документация

По биопитону вообще, и по пакету для работы с 3D-структурами в частности, имеется неплохая документация:

• Tutorial

• FAQ

Чтение файла

Первое, что в любом случае нужно делать при работе с трёхмерными структурами в биопитоне – это прочитать их. Делается это так:

   1 from Bio.PDB import PDBParser
   2 p = PDBParser()
   3 structure = p.get_structure("9zyx", "9zyx.pdb")
   4 header = p.get_header()
   5 trailer = p.get_trailer()

Метод get_structure получает два аргумента: логическое название структуры (например, её PDB-код) и имя файла, в котором лежит структура. Название структуры требуется биопитону исключительно для отладки (чтобы, если ваша программа одновременно работает со многими файлами, вы могли отличать, на атом из какого из них вы сейчас смотрите), поэтому здесь можно оставить и пустую строку.

Последние две строчки нужны, если вы хотите получать из PDB информацию заголовков. Эти строчки обязаны идти послe get_structure, иначе биопитон не будет знать, из какого файла мы хотим получить заголовки.

Биопитон разбирает информацию в заголовках так, чтобы ей было удобно пользоваться, и предоставляет её в виде словаря. Этот словарь попадёт в нашем примере в переменную header. Так как биопитон пытается разбирать её с пониманием происходящего, то попадает в этот словарь далеко не всё, что есть в заголовках. Зачастую, если задача связана с использованием информации в заголовках PDB-файла, приходится разбирать их самостоятельно, так, как мы это делали в прошлом блоке.

В trailer (в нашем примере) попадёт в виде списка строк всё, что идёт после раздела с атомами.

Обобщение здесь можно сделать такое: для большинства жизненных случаев PDB-файлы читаются двумя строками:

   1 parser = PDBParser()
   2 structure = parser.get_structure("", "9zyx.pdb")

Иерархия данных

Данные в структуре представляются в виде иерархии: структура -> модель -> цепь -> остаток -> атом.

Структура представляет все трёхмерные данные данные, которые есть в одном PDB-файле.

Модель представляет одно "изображение" из нескольких, которые хранятся в одном PDB-файле. Несколько "изображений" в PDB-файле погут появиться, например, если мы хотим сделать анимацию про то, как молекула движется (с этими возможностями у pymol вас через год познакомит А. В. Головин в курсе, посвящённом молекулярной динамике). Ещё один типичный случай – если структура молекулы появилась в результате эксперимента по ЯМР: в таком случае результаты эксперимента представляют последовательность разных возможных конформаций в порядке убывания их вероятности. Разные модели в одном PDB-файле могут содержать разные наборы молекул (это довольно типично для анимации). Всё же, как правило, в одном PDB-файле лежит одна модель, либо можно считать, что одна модель представляет всё, что нам от PDB-файла нужно.

Каждый уровень иерархии устроен как объект, который ведёт себя как словарь, но вдобавок имеет несколько полезных методов:

• get_parent() выдаёт словарь более высокого уровня. Например atom.get_parent() выдаст словарь всех атомов в остатке, если в переменной atom был атом.

• get_list() выдаёт содержимое словаря этого уровня в виде списка. Например chain.get_list() выдаст список всех остатков в цепочке.

• get_atoms() выдаёт список всех атомов, входящих в данный кусок структуры.

• get_residues() список всех остатков в данном куске структуры (есть только у уровней иерархии выше цепочки: у модели и структуры)

Дополнительно, есть атрибуты:

• У атомов:

  • name – название, например: CA, N, CB и т.п.

  • element – химический элемент: C, N, P, MG, CA (здесь это значит кальций) ...

  • coord – координаты атома. Про них подробнее ниже

  • bfactor, occupancy, ... – все данные, которые есть в PDB, здесь доступны; обращайтесь к документации за подробностями

• У остатков:

  • resname – имя остатка

Итак, теперь мы можем написать:

   1 from Bio.PDB import PDBParser
   2 p = PDBParser()
   3 s = p.get_structure("9zyx", "9zyx.pdb")
   4 model = s[0] # первая модель
   5 chain = model["A"] # первая цепочка в модели
   6 res = chain[10] # остаток с номером 10 согласно PDB-файлу
   7 atom = res["CA"] # C-alpha атом
   8 print atom.coord # напечатать координаты атома

Или короче:

   1 atom = s[0]["A"][10]["CA"]
   2 print atom.coord

Тонкости: коды вставки и гетероатомы

Когда фазовую проблему для рентгеноструктурного анализа решают методом замещения, зачастую стараются сохранить номера остатков совпадающими с номерами остатков земещающей молекулы. Однако, если в нашей молекуле имеются вставки (в смысле выравнивания: то есть такие остатки нашей молекулы, для которой нет гомологов у замещающей), то с таким подходом к нумерации возникли бы проблемы. Когда разрабатывали стандарт файла PDB, решили, что это очень существенная ситуация, и её нужно учитывать, поэтому добавили ещё одну колонку к номеру остатка: код вставки. В этой колонке может быть любая цифра или буква (что позволяет сохранять нумерацию при вставках длиной до 36 остатков)¹. Это несколько усложняет нам жизнь: до сих пор мы относились к номеру остатка как к числу, а тут внезапно оказалось, что это пара из числа и буквы.

Вторая проблема состоит в следующем: кроме аминокислотных остатков и нуклеотидов в структурах часто встречаются лиганды или вода² – в PDB файле их атомы называются "гетероатомами". С ними есть такая беда, что иногда им дают номера остатка, совпадающие с одним из имеющихся в структуре (например, символизируя этим, что данный лиганд именно с этим остатком взаимодействует).

Чтобы в таких условиях уметь нормально по номеру остатка обращаться к нему, в биопитоне принято следующее решение: номером остатка является тройка из указания гетероатома, номера, и указания кода вставки. Указание гетероатома бывает либо пустое ("" или " "), если это аминокислота или нуклеотид, бывает буква "W", если это вода, а если это лиганд, то здесь будет стоять "H_" и название остатка лиганда. (Иногда лиганды с разными названиями остатка приписывают к одному номеру остатка – к слову сказать, программы для молекулярного моделирования таким PDB-файлам крайне не рады и ругаются на них страшными недобрыми словами).

Итого, правильный способ получить 10-й остаток в цепочке "A" такой: res = model["A"]["", 10, " "]. Но, чтобы не писать таких безумных записей каждый раз, в биопитоне принято сокращение: можно использовать просто число, если мы имеем дело с остатком, у которого нет кода вставки, и он не является гетероостатком: res = model["A"][10] – как мы и писали в примере выше.

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

Тонкости: альтернативные конформации

Другая трудность, связанная с рентгеноструктурным анализом, создана уже не людьми. Состоит она в том, что в асимметрической ячейке кристалла один и тот же белок может лежать несколько раз в разных конформациях. И экспериментаторам иногда удаётся выяснить расположения атомов в разных конформациях. Эту информацию хочется отобразить в PDB. В силу особенностей технологии альтернативные конформации можно вообще как-либо заметить только в тех случаях, когда большая часть белка имеет одинаковую конформацию, поэтому и в PDB-файлах авторы структур предпочитают отображать это так, чтобы было видно, что вот такая часть молекулы неподвижна, а вот у такой части молекулы есть две конформации.

Для этого в формате PDB для атомов добавили ещё одно поле: "alter code". Для большинства атомов там стоит пробел, который нужно понимать так: альтернативных конформаций в этом месте нет. Для некоторых атомов там стоит буква – чаще всего, A или B.

Когда биопитон загружает такую структуру, он делает вид, что вы имееете дело только с одной конформацией – т.е. выбирает из доступных конформаций одну.

У всех атомов и остатков есть метод is_disordered(), который отвечает, есть ли альтернативные конформации у данного атома.

У тех атомов, у которых альтернативные конформации есть, метод is_disordered() вернёт истину, и вдобавок, у них есть методы:

• atom.disordered_select("B")  – выбрать для этого атома другую конформацию (в зависимости от этого поменяются его координаты, bfactor, occupancy и т.п.)

• atom.disordered_get_id_list() – возвращает список доступных кодов конформаций

Работа с координатами, углами и расстояниями

Самое полезное, что есть в 3D-структуре – это, разумеется, координаты атомов.

Они лежат в атоме в атрибуте coord – там лежит объект типа array из пакета numpy. Полное рассмотрение этого пакета выходит далеко за рамки нашего курса, здесь мы ограничимся простыми операциями:

• массивы ведут себя как списки, вы можете написать: x, y, z = atom.coord или len(atom.coord)³ или atom.coord[2] += 1

• если к массиву применить арифметическую операцию с числом, она применится поэлементно ко всему массиву: atom.coord * 2 это умножение вектора на 2

• если к массиву применить арифметическую операцию с другим массивом, она применяется поэлементно: atom1.coord - atom2.coord это вектор направления от atom2 к atom1.

Для подсчёта расстояний между атомами есть маленькое удобство: если сами атомы (а не их координаты) вычесть друг из друга, то биопитон вернёт их расстояние: dist = atom1 - atom2

А вот подсчёт углов нужно делать либо средтсвами numpy (это просто, но находится за пределами этого обсуждения), либо средствами биопитона, которые странные неуместные и неудобные. Чтобы посчитать угол, нужно координаты превратить сначала в "вектор" в понимании биопитона (для этого есть функция Vector), а уже к векторам можно применять функции для вычисления углов:

   1 from Bio.PDB import calc_angle, calc_dihedral, Vector
   2 a1 = calc_angle(Vector(atom1.coord), Vector(atom2.coord), Vector(atom3.coord))
   3 a2 = calc_dihedral(Vector(atom1.coord), Vector(atom2.coord), Vector(atom3.coord), Vector(atom4.coord))

Трансформации

В биопитоне есть минимальные удобства для движения частей молекулы. У любой части иерархии (начиная от всей структуры целиком) есть метод transform, который получает на вход матрицу преобразования (которое всегда сводится к комбинации из поворота, масштабирования и отражения) и вектор сдвига, и применяет сначала матрицу преобразования, затем вектор сдвига.

Вдобавок к ним прилагаются и функции для создания матриц для типичных преобразований:

• rotmat(p, q) возвращает матрицу поворота, которая переводит вектор p в вектор q. (Векторы p и q могут быть питонскими списками, массивами numpy, или векторами биопитона)

• rotaxis(alpha, axis) возвращает матрицу поворота вокруг оси, заданной вектором axis (вектор, возвращаемый функцией Vector), проходящей через начало координат, на угол alpha (в радианах)

• refmat(p, q) возвращает матрицу зеркального отражения относительно такой плоскости, чтобы вектор p перешёл в вектор q. (Векторы должны получаться из функции Vector).

Например, если мы хотим повернуть остаток res вокруг оси C-alpha – C-beta на 10 градусов, мы можем сделать это так:

   1 from math import radians
   2 from Bio.PDB import rotaxis, rotmat
   3 
   4 a = res["CA"].coord
   5 m = rotaxis(radians(10), res["CA"].coord - res["CB"].coord)
   6 id = rotmat(res["CA"].coord, res["CA"].coord)
   7 res.transform(id, -a) # сдвигаем C-alpha в начало координат
   8 res.transform(m, a) # поворачиваем матрицей m и сдвигаем в прежнее начало координат)

Либо мы можем просто менять вектор coord у атомов.

Сохранение

Поредактированную структуру можно сохранить в файл. Это делается примерно так же, как загрузка:

   1 from Bio.PDB import PDBIO
   2 writer = PDBIO()
   3 writer.set_structure(structure)
   4 writer.save("out.pdb")