- Предсказание активности ингибитора (подготовка файлов)
Дан набор из 88 веществ - ингибиторов тромбина (compounds.sdf). Для 85 из них активность известна, для трех - предстоит предсказать.
Для начала необходимо построить пространственное выравнивание активных конформаций исследуемых веществ. Будем считать активной конформацией (то есть конформацией, в которой вещество-ингибитор взаимодействует с белком-мишенью) наиболее энергетически выгодную конформацию (часто это вполне соответствует истине). Попробуем сгенерировать эти конформации, используя программу obconformer из пакета OpenBabel:
obconformer 100 100 compounds.sdf > compounds_best_conformer.sdf
Далее необходимо сделать выравнивание полученных конформеров. Попробуем сделать это с помощью программы Open3DALIGN (open3dalign.sourceforge.net).
open3dalign.sh
Чтобы сделать выравнивание нужно загрузить SDF файл со структурами веществ (команда import), выполнить выравнивание (чтобы использовать в качестве темплэйта, к которому выравниваются все вещества, первое вещество в списке, наберите align object_list=1), и записать выравнивание в файл (save).
Перекодировать из юникода в 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
- 3DQSAR анализ
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) находится слишком близко к атомам исследуемых содеинений, и дает слишком большую по модулю энергию. Установим ограничения на значения энергии:
cutoff type=max level=5.0 field_list=1
cutoff type=min level=-5.0 field_list=1
Слишком маленькие значения энергии приравняем к 0:
zero type=all level=0.05
Исключим из анализа ячейки, в которых вариабельность в энергии взаимодействия с зондом для разных соединений мала:
sdcut level=0.1
nlevel
remove_x_vars type=nlevel
Построим регрессионную модель:
pls
Получаем коэффициенты корреляции (r2) для пяти компонент, выделенных PLS. И они нас устраивают, для 3х компонент они близки к 1.
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
Кросс-валидация:
cv type=loo runs=20
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
Здесь коэффициенты - плохие, близки к нулю.
Предсказание активности для тестовой выборки:
predict
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
Здесь коэффициенты корреляции r2(pred) не такие плохие, как q2.
- Анализ, используя выравнивание и конформации, полученные с учетом структуры активного центра белка-мишени
Регрессионный анализ
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
Кросс-валидация
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
Предсказание активности
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
Коэффициенты улучшились при кросс-валидации и при анализе тестовой выборки (здесь незначительно).
- Предсказание активности
Для начала переделаем модель с использованием всех имеющихся данных, а вещества с неизвестной активностью обозначим как тестовую выборку:
set object_list=60-85 attribute=TRAINING
set object_list=86-88 attribute=TEST
Построение модели
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.8375 12.8375 44.4004 44.4004 0.7269 0.4440
2 14.5264 27.3638 14.3748 58.7753 0.6260 0.5878
3 6.9607 34.3245 11.2007 69.9760 0.5342 0.6998
4 8.4659 42.7904 5.4939 75.4699 0.4828 0.7547
5 4.7600 47.5503 5.7466 81.2166 0.4225 0.8122
Предсказание
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.0954 7.5090 7.3772 7.6623 7.8822 1
87 87 44 0.0000 6.9300 7.0808 6.9883 7.1990 7.4119 1
88 88 72 0.0000 5.5493 5.2836 5.1285 5.3788 5.3537 3
PC r2(pred) SDEP
--------------------------
0 0.0000 6.6604
1 0.0294 6.5616
2 -0.0102 6.6942
3 0.0265 6.5717
4 -0.0480 6.8183
5 -0.0950 6.9696
Кросс-валидация
PC SDEP q2
--------------------------
0 0.9865 -0.0240
1 0.8305 0.2743
2 0.7666 0.3816
3 0.7394 0.4247
4 0.7339 0.4332
5 0.7380 0.4269
Коэффициент q2 наилучший у компоненты 4 (кросс-валидация выдала лучше коэффициенты, чем r2(pred) в предсказании), поэтому предскажем активность, как:
-------------------
N ID 4
-------------------
86 86 7.6623
87 87 7.1990
88 88 5.3788