Kodomo

Пользователь

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

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

Суть задания построить 3DQSAR модель для ингибиторов тромбина и предсказать активность для трех веществ, активность которых не известна.

Примечание: приводите в отчете команды, которые выполняете, либо прикладывайте скрипт!

  1. Для проведения 3DQSAR анализа мы будем использовать программы Open3DQSAR и Open3DALIGN (open3dqsar.sourceforge.net). Чтобы использовать их на kodomo, добавьте к переменной PATH директорию /home/preps/grishin/open3dtools/bin)

   1 export PATH=$PATH:/home/preps/grishin/open3dtools/bin

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

  1. Итак, вам дан набор из 88 веществ – ингибиторов тромбина (compounds.sdf). Для 85 из них активность известна, для трех – нам предстоит предсказать.


Для начала необходимо построить пространственное выравнивание активных конформаций исследуемых веществ. Будем считать активной конформацией (то есть конформацией, в которой вещество-ингибитор взаимодействует с белком-мишенью) наиболее энергетически выгодную конформацию (часто это вполне соответствует истине). Попробуем сгенерировать эти конформации, используя программу obconformer из пакета OpenBabel:

   1 obconformer 100 100 compounds.sdf > compounds_best_conformer.sdf

К сожалению, это занимает порядка 40 минут, поэтому лучше взять уже готовый файл (compounds_best_conformer.sdf).


Далее необходимо сделать выравнивание полученных конформеров. Попробуем сделать это с помощью программы Open3DALIGN (open3dalign.sourceforge.net). Документацию смотрите на сайте.

Запустите программу:

   1 open3dalign.sh 

Чтобы сделать выравнивание нужно загрузить ваш SDF файл со структурами веществ (команда import), выполнить выравнивание (чтобы использовать в качестве темплэйта, к которому выравниваются все вещества, первое вещество в списке, наберите align object_list=1), и записать выравнивание в файл (save). Синтаксис команд смотрите на сайте в разделе документация.

Technical note: программа Open3DALIGN сохраняет SDF файл c выравненными веществами в кодировке вашей локали. В результате, если ваша локаль, например, юникод, PyMOL такие файлы может отказаться читать.

Посмотрите на выравнивание в PyMOL (попробуйте команду split_states чтобы отобразить все вещества одновременно).Если все работает, переходите к пункту 3.

Если нет, попробуйте следующие команды (в командной строке юникс):

Перекодировать из юникода в ascii:

   1 iconv -c -f utf-8 -t ascii aligned.sdf > aligned_ascii.sdf

Удалить ненужную информацию из заголовков и добавить $$$$ в конец каждой записи:

   1 sed -e 's/.*HEADER.*\([0-9][0-9]\).*/\1/' -e 's/\(.*M  END.*\)/\1\n$$$$/'  aligned_ascii.sdf > temp
   2 sed -n '/^[0-9a-zA-Z \$\.-]*$/ p'  temp > aligned_ok.sdf
   3 rm temp
  1. Теперь мы можем попробовать выполнить непосредственно 3DQSAR анализ и посмотреть, получается ли построить регрессионную модель с помощью такого выравнивания.

Запустите программу:

   1 open3dqsar.sh

Все дальнейшие команды вводятся в командной строке open3dqsar.
Загрузите файл со структурами:

   1 import  type=sdf file=aligned_ok.sdf

Загрузите файл с данными об активности исследуемых соединений (activity.txt):

   1 import type=dependent file=activity.txt

Обратите внимание, активности трех последних соединений нам предстоит предсказать, поэтому для них пока что указана нулевая активность.

Задайте решетку вокруг исследуемых соединений:

   1 box

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

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

   1 set object_list=60-85 attribute=TEST
   2 set object_list=86-88 attribute=EXCLUDED

Рассчитаем значения энергии ван-дер-Ваальсовых взаимодействий в узлах решетки:

   1 calc_field type=VDW force_field=MMFF94 probe_type=CR

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

   1 cutoff type=max level=5.0 field_list=1
   2 cutoff type=min level=-5.0 field_list=1

Слишком маленькие значения энергии приравняем к 0:

   1 zero type=all level=0.05

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

   1 sdcut level=0.1
   2 nlevel
   3 remove_x_vars type=nlevel

Наконец, построим регрессионную модель:

   1 pls 

В результате выполнения этой программы вы получите коэффициенты корреляции для разного количества компонент, выделенных PLS.
Какие коэффициенты вы получили? Можно что-то предсказать с помощью такой модели?

Если вас устраивает модель, попробуйте выполнить кросс-валидацию:

   1 cv type=loo runs=20

И предсказать активность для тестовой выборки:

   1 predict

Отметьте разницу в коэффициентах корреляции, полученных при построении модели, при кросс-валидации и при анализе тестовой выборки.

  1. Теперь попробуем выполнить тот же анализ, но используя выравнивание и конформации, полученные с учетом структуры активного центра белка-мишени (на самом деле они находятся в исходном файле compounds.sdf). Посмотрите на выравнивание в PyMOL.

Повторите 3DQSAR анализ с этим выравниванием. Отметьте разницу в полученных коэффициентах корреляции (полученных при построении модели, при кросс-валидации и при анализе тестовой выборки). Попытайтесь объяснить разницу.

  1. Давайте попробуем использовать получившуюся модель (даже если она не очень хороша) для предсказания активности.

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

   1 set object_list=60-85 attribute=TRAINING
   2 set object_list=86-88 attribute=TEST

Построим модель и предскажем активность трех веществ:

   1 pls
   2 predict

Запишите в отчет эти активности. Чем больше они будут совпадать в реальными (которые мне на самом деле известны), тем выше балл :) .

6*) Давайте теперь посмотрим на вашу модель в PyMOL. После того, как вы построили модель командой pls, наберите команду:

   1 export type=coefficients pc=5 format=insight file=coefs

Откройте PyMOL, загрузите туда файл со структурами, а также получившийся в результате экспорта файл coefs_fld-01_y-01.grd. В этом файле находится информация о коэффициентах модели в каждой точки решетки. Такого рода данные можно визуализовать в PyMOL с помощью команды isosurface, которая рисует поверхность, построенную из всех точек, в которых коэффициенты равны заданному значению:

   1 isosurface positive, coefs_fld-01_y-01, 0.00001

для отображения области, соответствующей коэффициентам 0.00001, и

   1 isosurface negative, coefs_fld-01_y-01, -0.00001 

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

7*) Попробуйте добавить в модель электростатические взаимодействия (calc_field type=MM_ELE, остальные параметры смотрите в документации). Помогают ли они улучшить модель? К сожалению, в Open3DQSAR не реализованы водородные связи и гидрофобные взаимодействия.

8*) Модель «обучается» только тому, что вы ей дали в качестве training set. Поэтому очень важно, чтобы вещества, используемые при построении модели, охватывали максимальный спектр модификаций основной структуры, то есть были максимально разнообразными. Если в тестовой выборке будут вещества, содержащие модификации, которых (или сходных) не было в обучающей выборке, хорошего предсказания не получится! Попробуйте построить несколько моделей, используя в качестве обучающей выборки разные наборы веществ, и следите за качеством модели. Свои наблюдения запишите в отчет.

Примечание: на самом деле OpenBabel не слишком хорош по части генерации конформаций (да и не только по этой части), и я бы не рекомендовал его использовать для реальных задач (разве что когда вы действительно знаете, что делаете :) ). Кроме того, он выдает только одну конформацию, в то время как обычно есть смысл использовать набор (50-200) различных низкоэнергетических конформаций. Во многих коммерческих пакетах программ такая возможность реализована. Одна из лучших (и самая быстрая) программа для генерации конформации - OMEGA от OpenEye (если захотеть, ее можно получить бесплатно). Есть еще смысл попробовать бесплатный (но не открытый) Balloon.