QSAR | Наталья Ланина

Для проведения 3DQSAR анализа мы будем использовать программы Open3DQSAR и Open3DALIGN (open3dqsar.sourceforge.net).

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

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

wget https://kodomo.fbb.msu.ru/~golovin/qsar/compounds.sdf --no-check-certificate

1. Построение пространственных выравниваний исследуемых веществ.

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

obconformer 100 100 compounds.sdf > compounds_best_conformer.sdf

Далее необходимо сделать выравнивание полученных конформеров с помощью программы Open3DALIGN (open3dalign.sourceforge.net). Для этого нужно загрузить SDF файл со структурами веществ (import), выполнить выравнивание (образец - первое вещество в списке, align object_list=1), и записать выравнивание в файл (save):

open3dalign.sh
import type=SDF file=compounds_best_conformer.sdf
align object_list=1
save file=aligned.sdf

Перекодируем из UNICODE в ASCII:

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

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

sed -e 's/.*HEADER.*\([0-9][0-9]\).*/\1/' -e 's/\(.*M  END.*\)/\1\n$$$$/'  aligned_ascii.sdf > temp
sed -n '/^[0-9a-zA-Z \$\.-]*$/ p'  temp > aligned_ok.sdf
rm temp

Получили файл выравнивания aligned_ok.sdf.

Рисунок 1. Пространственное выравнивание исследуемых молекул.

2. 3DQSAR анализ.

Загружаем файл с данными об активности исследуемых соединений (activity.txt) и пробуем построить регрессионную модель с помощью полученного выравнивания:

wget https://kodomo.fbb.msu.ru/~golovin/qsar/activity.txt --no-check-certificate
open3dqsar.sh
import  type=sdf file=aligned_ok.sdf
import type=dependent file=activity.txt
box

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

set object_list=60-85 attribute=TEST
set object_list=86-88 attribute=EXCLUDED
calc_field type=VDW force_field=MMFF94 probe_type=CR

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

cutoff type=max level=5.0 field_list=1
cutoff type=min level=-5.0 field_list=1
zero type=all level=0.05

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

sdcut level=0.1
nlevel
remove_x_vars type=nlevel

Строим регрессионную модель:

pls

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

          Exp.   Cum. exp.        Exp.   Cum. exp.
PC    var. X %    var. X %    var. Y %    var. Y %        SDEC          r2
--------------------------------------------------------------------------
 0      0.0000      0.0000      0.0000      0.0000      0.9494      0.0000
 1     15.9480     15.9480     32.8386     32.8386      0.7780      0.3284
 2      5.1333     21.0813     36.3625     69.2011      0.5269      0.6920
 3      4.6235     25.7048     15.6991     84.9002      0.3689      0.8490
 4      3.8908     29.5956      7.5246     92.4248      0.2613      0.9242
 5      4.0108     33.6064      2.8661     95.2909      0.2060      0.9529

Elapsed time: 0.2991 seconds.

Все коэффициенты корреляции получились больше 0, а несколько можно даже округлить до 1. Модель можно считать хорошей, она нас устраивает. Выполняем кросс-валидацию:

cv type=loo runs=20

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

predict

Результат кросс-валидации - все коэффицинты корреляции отрицательны:

PC        SDEP          q2
--------------------------
 0      0.9658     -0.0348
 1      0.9164      0.0683
 2      0.9733     -0.0509
 3      0.9667     -0.0368
 4      0.9880     -0.0829
 5      0.9497     -0.0006

Elapsed time: 0.9765 seconds.

Результат предсказания коэффициентов корреляции - они все положительны, но близки к нулю:

PC    r2(pred)        SDEP
--------------------------
 0      0.0000      1.0362
 1      0.2655      0.8881
 2      0.3296      0.8484
 3      0.2353      0.9061
 4      0.2754      0.8821
 5      0.2536      0.8953
 
 Elapsed time: 0.0205 seconds.

3. 3DQSAR анализ c использованием выравнивания и конформаций, полученных с учетом структуры активного центра белка-мишени.

open3dalign.sh
import type=SDF file=compounds.sdf
align object_list=1
save file=aligned_2.sdf

iconv -c -f utf-8 -t ascii aligned_2.sdf > aligned_ascii_2.sdf
sed -e 's/.*HEADER.*\([0-9][0-9]\).*/\1/' -e 's/\(.*M  END.*\)/\1\n$$$$/'  aligned_ascii_2.sdf > temp
sed -n '/^[0-9a-zA-Z \$\.-]*$/ p'  temp > aligned_ok_2.sdf
rm temp

Рисунок 2. Пространственое выравнивание с учетом структуры активного центра белка-мишени.

open3dqsar.sh
import  type=sdf file=aligned_ok_2.sdf
import type=dependent file=activity.txt
box
set object_list=60-85 attribute=TEST
set object_list=86-88 attribute=EXCLUDED
calc_field type=VDW force_field=MMFF94 probe_type=CR
cutoff type=max level=5.0 field_list=1
cutoff type=min level=-5.0 field_list=1
zero type=all level=0.05
sdcut level=0.1
nlevel
remove_x_vars type=nlevel
pls
cv type=loo runs=20
predict

          Exp.   Cum. exp.        Exp.   Cum. exp.
PC    var. X %    var. X %    var. Y %    var. Y %        SDEC          r2
--------------------------------------------------------------------------
 0      0.0000      0.0000      0.0000      0.0000      0.9494      0.0000
 1     12.1342     12.1342     48.4736     48.4736      0.6815      0.4847
 2     13.2295     25.3637     14.5885     63.0621      0.5770      0.6306
 3      7.6412     33.0049     13.2040     76.2661      0.4625      0.7627
 4      8.0257     41.0305      4.3684     80.6345      0.4178      0.8063
 5      6.0521     47.0827      3.8642     84.4987      0.3738      0.8450

Elapsed time: 0.1493 seconds.

Результат кросс-валидации^ только один коэффицинт отрицателен, но все же больше, чем полученные по PLS. Хуже, чем в прошлый раз:

PC        SDEP          q2
--------------------------
 0      0.9658     -0.0348
 1      0.8027      0.2851
 2      0.7664      0.3484
 3      0.7061      0.4468
 4      0.6735      0.4968
 5      0.6401      0.5454

Elapsed time: 0.1001 seconds.

Результат предсказания коэффициентов корреляции - практически не отличаются от предыдущего варианта:

PC    r2(pred)        SDEP
--------------------------
 0      0.0000      1.0362
 1      0.3451      0.8385
 2      0.3226      0.8529
 3      0.2998      0.8671
 4      0.3012      0.8662
 5      0.2693      0.8858


Elapsed time: 0.0066 seconds.

4. Использование полученной модели для предсказания активности

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

open3dqsar.sh
import  type=sdf file=aligned_ok_2.sdf
import type=dependent file=activity.txt
box
set object_list=60-85 attribute=TRAINING
set object_list=86-88 attribute=TEST
calc_field type=VDW force_field=MMFF94 probe_type=CR
cutoff type=max level=5.0 field_list=1
cutoff type=min level=-5.0 field_list=1
zero type=all level=0.05
sdcut level=0.1
nlevel
remove_x_vars type=nlevel
pls
predict

Результат PLS:

          Exp.   Cum. exp.        Exp.   Cum. exp.
PC    var. X %    var. X %    var. Y %    var. Y %        SDEC          r2
--------------------------------------------------------------------------
 0      0.0000      0.0000      0.0000      0.0000      0.9749      0.0000
 1     12.5822     12.5822     46.4042     46.4042      0.7137      0.4640
 2     14.2226     26.8048     15.5157     61.9199      0.6016      0.6192
 3      6.7847     33.5895     11.1828     73.1027      0.5056      0.7310
 4      8.7614     42.3509      4.2898     77.3925      0.4635      0.7739
 5      4.7029     47.0537      4.5965     81.9889      0.4137      0.8199

Elapsed time: 0.2082 seconds.

Результат предсказания:

External predictions for dependent variable  1 (activity)
---------------------------------------------------------------------------------------------------
    N   ID    Name  Actual           1           2           3           4           5    Opt PC n
--------------------------------------------------------------------------------------------------
   86   86    01    0.0000      7.1119      7.5466      7.4119      7.6262      7.7234           1
   87   87    44    0.0000      6.9428      7.1202      7.0946      7.3278      7.5477           1
   88   88    72    0.0000      5.5073      5.2436      5.1697      5.4378      5.4696           3
   
   PC    r2(pred)        SDEP
--------------------------
 0      0.0000      6.6604
 1      0.0298      6.5603
 2     -0.0155      6.7118
 3      0.0082      6.6331
 4     -0.0627      6.8660
 5     -0.1011      6.9889

Elapsed time: 0.0097 seconds.

Все файлы, полученные во время выполнения практикмуа, можно найти в папке.

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

1. Построение пространственных выравниваний исследуемых веществ.

2. 3DQSAR анализ.

3. 3DQSAR анализ c использованием выравнивания и конформаций, полученных с учетом структуры активного центра белка-мишени.

4. Использование полученной модели для предсказания активности