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

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

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

obconformer 100 100 compounds.sdf > compounds_best_conformer.sdf

Далее необходимо сделать выравнивание полученных конформеров с помощью программы Open3DALIGN.

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
In [1]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import time

cmd = ServerProxy(uri="http://localhost:9123/RPC2")

cmd.bg_color('white')
Dir = '/Users/aleksandrabezmenova/Documents/fbb/term8/Golovin/Pr7/'

ImPath = '/tmp/pymolpng.png'
def MakeImage():
    cmd.ray(500,400)
    cmd.png(ImPath)
    time.sleep(2)
    
def focus(sele):
    cmd.center(sele)
    cmd.zoom(sele)
In [2]: 
MakeImage()
Image(ImPath)
Out[2]:

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

          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

При построении модели коэффициенты достигают почти 1, в предсказании коэффициенты меньше, при кросс-валидации - примерно 0.

Выполним тот же анализ, но используя выравнивание и конформации, полученные с учетом структуры активного центра белка-мишени.

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

iconv -c -f utf-8 -t ascii aligned1.sdf > aligned_ascii1.sdf
sed -e 's/.*HEADER.*\([0-9][0-9]\).*/\1/' -e 's/\(.*M  END.*\)/\1\n$$$$/'  aligned_ascii1.sdf > temp
sed -n '/^[0-9a-zA-Z \$\.-]*$/ p'  temp > aligned1_ok.sdf
rm temp
In [3]: 
MakeImage()
Image(ImPath)
Out[3]:
open3dqsar.sh
import  type=sdf file=aligned1_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
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

          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
 
cv type=loo runs=20

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
 
predict

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

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

pls
predict

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
 

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

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

Соединение Активность
86 7.0954
87 6.9300
88 5.5493