Пример использования трехмерного 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
Построение пространственных выравниваний исследуемых веществ
Для начала необходимо построить пространственное выравнивание
активных конформаций исследуемых веществ. Будем считать активной
конформацией (то есть конформацией, в которой вещество-ингибитор взаимодействует
с белком-мишенью) наиболее энергетически выгодную конформацию
(часто это вполне соответствует истине).
Попробуем сгенерировать эти конформации, используя программу 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
Перекодировать из юникода в 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. Пространственное выравнивание исследуемых веществ |
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.3070 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: 1.0015 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.0185 seconds.
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.1526 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.0985 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.0065 seconds.
Использование полученной модели для предсказания активности
Для начала, переделаем модель с использованием всех имеющихся данных, а вещества с неизвестной активностью обозначим как тестовую выборку, затем построим модель и предскажем активность трех веществ:
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.2171 seconds.
Результат предсказания:
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.0070 seconds. External predictions for dependent variable 1 (activity) -------------------------------------------------------------------------------------- N Name Actual 1 2 3 4 5 OptPC n -------------------------------------------------------------------------------------- 86 01 0.0000 7.1119 7.5466 7.4119 7.6262 7.7234 1 87 44 0.0000 6.9428 7.1202 7.0946 7.3278 7.5477 1 88 72 0.0000 5.5073 5.2436 5.1697 5.4378 5.4696 3
Сравнивая активность с коэффициентами корреляции, можно заключить с использованием первой компоненты:
N Activity ------------- 86 7.1119 87 6.9428 88 5.5073
Скрипт для получени картинок в PyMol - script.pml.