专栏名称: 生信杂谈
生物信息学;生物信息;计算机辅助药物设计;测序分析;Python;R;机器学习;论文写作;网站制作;LOL;dota2。
目录
相关文章推荐
彬彬有理  ·  身体出现这几个症状,要提早重视 ·  20 小时前  
宛央女子  ·  穿真丝内衣的福报来了 ·  昨天  
宛央女子  ·  用永远天真无畏的杉菜送别大S ·  4 天前  
彬彬有理  ·  过年社交潜规则,没人明说,但很重要 ·  4 天前  
重庆共青团  ·  25岁徐枫灿,二级机长! ·  4 天前  
重庆共青团  ·  25岁徐枫灿,二级机长! ·  4 天前  
51好读  ›  专栏  ›  生信杂谈

SMINA教程以CDPK1为例

生信杂谈  · 公众号  ·  · 2017-11-12 11:30

正文

SMINA教程以CDPK1为例,由于微信限制所有的链接和资料可以点击阅读原文下载

原作者:David Ryan Koes
翻译:康文渊

若进行精度建议查看原文以及参考文献,本文是我根据教程快速阅读做的笔记

Q1:我们有autodock vina,为什么还要SMINA?

A1:SMINA有更加优秀的速度和更加的扩展性

Q2:用该对接软件有什么优势?

A2:首先该对接软件为免费/学术免费软件,其次其有非常良好的对接准确度,甚至在商业软件中都为上乘,同时可以很好的运用SMINA进行得分改造,进一步提升性能和速度

本教程是一个使用SMINA进行对 E.tenella 中的钙依赖性蛋白激酶1(calcium-dependent protein kinase 1 CDPK1)进行的结构为基础的虚拟筛选的教程。 E.tenella 是一种感染幼禽并且导致潜在致命球虫病的寄生虫。尽管这种酶并未解析晶体结构(译者翻译时该结构已经被解析,pdb:4YSJ),但是有相关寄生虫的该种酶的结构。 C. parvum T. gondii 的该种酶结构已经被解析,并且有生物活性数据表明他们的小分子竭抗剂。
我们将会对怎么分析已知结构,评价SMINA对于这些结构的表现,对于对接结果创建自定义的得分函数,使用SWISS-MODEL创建 E.tenella 的结构模型,并且进行 E.tenella CDPK1的虚拟筛选。(个人认为作者提到的 E.tenella 均表示的为 E.tenella CDPK1)。

假设和惯例

假设读者对Linux命令行有一定的了解,并对一些用到的软件熟悉。教程在Ubuntu Linux12.04(原作者)和Ubuntu Linux16.04(译者)测试通过。bash命令行如下:


1

echo "Hello"



python基础的PyMOL环境展示如下:


1

2

#python2.x

print "Hello"



1. 检查 C.parvum T.gondii 晶体结构

该部分内容可以参考早期博客文章同源建模系列
首先我们分析 E.tenella 序列和验证三个结构的异同(三个寄生虫中改结构)。然后我们使用 C.parvum T.gondii 创建一个虚拟筛选的测试集,并且对接这些测试集化合物。我们使用这些对接结果来进行晶体结构的筛选和模型的增加。

1.1 靶标分析

首先,我们分析 E.tenella 序列和 C.parvum T.gondii* 结构

1.1.1 E.tenella 序列分析

1.复制或者下载 E.tenella 的FASTA序列,你可以点击这里进入网址,同样也可以下载我们已经下载好的。

2.进行BLAST搜索,若你为进入的网址,可以直接点击ncbi右侧的BLAST进行搜寻(如下图)。若为下载的文件可以进入BLAST页面上传fasta文件进行搜寻。

3.进入页面。

4.点击 BLAST

5.结果如图3,可以发现得分最高的为已经解析的结构,相似度为99%,但是我们在这里假装没有。排第二的为Neospora Caninum(4M97),但是更加相近的为 T.gondii 的几个结构(4M7N,3KU2,3I79,3HX4,3I7C),相似度大于80%。 C.parvum 的首次命中(3IGO)相似度为62%。


从中我们得出的结论为 T.gondii 的CDPK1相比于 C.parvum 具有更高的相似度。建议 T.gondii 结构选为同源建模的模板更为合理。

1.1.2 结构分析

为了鉴定 T.gondii E.tenella 之间CDPK1结合位点的差别,我们创建了一个 E.tenella 同源模型。并且使用color_by_mutation和PyMOL进行结构可视化序列的不同。
1.在SWISS-MODEL进行自动建模,你也可以参考之前的教程

译者注:我认为该方法并不是建模的最佳方法,最佳方法为找到具有共同抑制剂的晶体结构进行基于配体的多模板建模,这样可以保证重要的活性腔的相似,而不是使用没有结合配体的晶体结构进行建模。

2.输入 E.tenella 的FASTA序列和3I79进行建模(图4)

3.下载结果文件,你可以在这里下载我已经建模好的。

4.在PyMOL中加载模型,3I79和3NCG


1

2

3

load swissmodel.pdb

fetch 3I79

fetch 3NCG


译者注:在这里 3NCG 来的有点突然,搜索一下可以发现3NCG为 C.parvum 的激活状态下的晶体结构

5.比对结构


1

alignto


6.从3NCG中提取BK1配体,并且隐藏该pdb的其余结构,因为3I79为非结合状态,我们使用这个配体来占有结合位点


1

2

extract bk1,resn BK1

hide (/3NCG/)


7.探索模型和3I79结合位点周围的残基,值得注意的是SWISS-MODEL忠实的复制了侧链构像。

8.着色两者之间的突变位点


1

2

run color\_by\_mutation.py

color\_by\_mutation swissmodel,3i79


9.探索结合位点,注意的是3I79中GLY128在 E.tenella 模型中被突变为THR 108,其余在结合位点是大体不变的。

10.重复分析3NCG,可以发现 C.parvum 在结合位点附近具有更多的突变

1.2 测试集的创建

PubChem上有 T.gondii C.parvum 的CDPK1诘抗击。我们使用这些数据作为虚拟筛选的测试集对这两个靶标进行筛选.

1.2.1组装活性化合物

(a)在PubChem搜索CDPK1然后排序活性化合物,每个靶点选择三个最大的化合物数据集,为了保证和教程一致,我们选择的为和原教程相同的数据集.
T.gondii
aid:677062,600378,672956
C.parvum
aid:600379,672955,66097

(b)选择active,然后点击Structures下载,之后选择没有压缩的SMILES格式,可以点此下载我们下载好的

(c)修正SMILES结构。下载的SMILES文件字符串和标题顺序相反。另外,我们需要去除salt,并将关键字“active”添加到所有活性化合物的标题中。我们每个文献执行如下命令(例如aid677062.txt为例):


1

awk '{print $2,$1"_active"}' aid677062.txt | sed 's/\.Cl//' > aid677062.smi



awk主要是进行了前后的颠倒,并且在标题后增加active标签,sed主要是将Cl进行替换,替换成空,中间用管道链接,输出到smi文件格式。

(d)将每个靶点的内容合并为单个文件


1

cat aid*.smi > activaes.msi


T.gondii C.pavum 分别包含156个和89个活性化合物。你可以在这里下载合并好的文件。

1.2.2 创建诱饵数据集

该在线网站说只需要10min,实际上需要1-2天左右时间
由于这些靶点的无活性配体数量较少,我们使用诱饵数据库Database of Useful Decoys:Enhanced(DUDE)方法来从ZINC数据库中进行采样。这些诱饵被选择的方法为化学性质不相同(即和有活性的配体相比较不相似,认为是不结合的)。但分子质量,logP,旋转键和氢键受体供体这些简单的分子参数是相同的。

(a)使用DUDE网站创建诱饵数据集,将会生成dude-decoys.tar.gz文件,里面包含50个参数匹配的诱饵。

(b)将诱饵结合进入单个文件


1

cat decoys/* | grep -v ligand > decoys.msi


点此下载decoys文件,( T.gondii 诱饵和活性化合物。 C.pavum 诱饵和活性化合物)。

1.2.3 生成三维结构

我们将使用开源软件RDKit利用2D SMILES生成3D构像

(a)将RDKit写入你的PYTHONPATH,并且执行rdconf.py脚本


1

2

wget http://bits.csb.pitt.edu/tdtCDPK1/rdconf.py

chmod + x rdconf.py


(b)每个active/decoy生成单一构像


1

2

rdconf.py --maxconfs 1 decoys.smi decoys.sdf

rdconf.py --maxconfs 1 actives.smi actives.sdf


(c)将文件合并


1

cat actives.sdf decoys.sdf > combined.sdf


由于有两个文件,我们将 T.gondii C.pavum 分别命令为gondii.sdf和parvum.sdf

1.3 对接

我们使用SMINA通过对接化合物进入受体结构进行虚拟筛选,采用AutoDock Vina得分功能。由于我们对一个固定的受体对接,所以我们选择一个好的受体结构是非常重要的。我们将会使用CDPK1的测试集来模拟评价我们的对接工具和选择受体结构。

1.鉴定PDB中所有 C.parvum T.gondii 文件,搜索’calcium-dependent protein kinase 1’,并且选择恰当的生物(organism.)

C. parvum : 2QG5 2WEI 3DFA 3F3Z 3HKO 3IGO 3L19 3LIJ 3MWU 3NCG

T. gondii : 3I79 3I7B 3I7C 3KU2 3N51 3NYV 3SX9 3SXF 3T3U 3T3V 3UPX 3UPZ 3V51 3V5P 3V5T 4M84

由于教程较老,故和自己搜索有差异

2.下载这些pdb文件,你可以用如下的简单方法,或者直接下载我下载好的。

将以下代码保存为脚本,例如 down.sh


1

2

3

4

5

6

7

8

9

#!/bin/bash  

str="2QG5 2WEI 3DFA 3F3Z 3HKO 3IGO 3L19 3LIJ 3MWU 3NCG"  

arr=(${str// / })

for i in ${arr[@]}  

do  

wget https://files.rcsb.org/download/$i.pdb

done



执行如下命令:


1

2

chmod +x down.sh

./down.sh


3.比对和提取结构。在pymol中打开靶标pdb文件。比对他们。除去水和原子。提取每个配体到自己的对象。


1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

#pymol中执行

from glob import glob

for fil in glob("*.pdb"): cmd.load(fil)

alignto

remove solvent

remove metal

remove GOL

#选择配体BK1周围5埃外的小分子,并且删除,因为并未在结合位点内,不是配体

select unligand,byres resname BK1 gap 5 and hetatm

remove unligand

#选择配体

select ligand,hetatm

extract parvum,ligand

#保存配体

save parvumlig.pdb,parvum

remove ligand


4.保存受体和配体文件,注意的是 3LI9 没有酶结构域, 3DFA 2QG5 为未结合态结构,综合考虑我们将其删除。


1

remove 3LI9 and 3DFA and 2QG5


保存(注意一下需要保存为python文件后 run 进行操作)


1

for name in cmd.get_names():cmd.save(name + ".pdb" , name )


5.将配体导入进单一文件,由于我们之前做了,该步骤无需做

6.对接测试集化合物到每个受体结构,例如:


1

2

smina --seed 0 --autobox_ligand parvumlig.pdb -r 2WEI.pdb \

-l combined.sdf -o 2WEI_docked.sdf.gz


在这里我们制定了一个随机种子。结合盒子采用 autobox_ligand 设置,创建的盒子大小为对于提供的配体8埃缓冲。配体,受体和输出文件可以使用任意OpenBabel支持的格式。

7.获取每个化合物输出得分最高的对接姿势。


1

2

sdsorter -sort minimizedAffinity -reduceconfs 1 2WEI_docked.sdf.gz \

-print -c > 2WEI_docked.txt


8.分析对接结果,原文是使用的R进行的分析,大家可以以下做参考,当然也可以依葫芦画瓢使用Python进行处理。

(a)R版本方法如下:


1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

install.packages("pROC")

library("pROC")

rarr=list()

for (file in Sys.glob("*docked.txt")) {

d = read.table(file,header=T)

r = roc(grepl("active",d $Title) , d$minimizedAffinity,direction=">")

rarr[[length(rarr)+1]]=r

}

names(rarr) = Sys.glob("*docked.txt")

for (n in names(rarr)) {

r = rarr[[n]]

a = auc(r,partial.auc=c(1,.9),partial.auc.correct= T )

s = sprintf("%s %.4 f%.4f\n" , sub("_docked.txt" ,"",n),r$auc[1],a[1])

cat(s)

}


图大致如下:

AUC和Partial AUC的值分别为 0.8520 , 0.6485 。此数据运用的为官方数据,python版本为我自己写的方法,供参考

(b) ython版本如下:

首先我们使用 awk 处理一下原始文本


1

2

awk 'BEGIN{FS='_';OFS='\t'}{print $1,$2}' 2WEI_docked.txt rocdata.txt

#打开文件,在Title后加type项


得到的文件可以从这里下载,python脚本如下:


1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

#-*-coding:utf-8-*-

#python=3.x

import pandas as pd

from sklearn.metrics import roc_curve, roc_auc_score

import matplotlib.pyplot as plt

# 读取数据

df=pd.read_csv('rocdata.txt',sep='\t')

df.loc[lambda df: df.type=='active',['type']]='1'

df.loc[lambda df: df.type=='unactive',['type']]='0'

df['minimizedAffinity']=df['minimizedAffinity']/-10.76520

y=df['type'].astype('int64')

scores=df['minimizedAffinity']

# ROC

fpr,tpr,thresholds=roc_curve(y,scores,pos_label=1)

#auc值

auc=roc_auc_score(df['type'].astype('int64'),scores)

#绘图

plt.figure()

plt.plot(fpr,tpr,color='darkorange',lw=2,label='pavum ROC curve,auc %0.2f' % auc)

plt.xlim([0.0,1.0])

plt.ylim([0.0,1.05])

plt.plot([0,1],[0,1],color='navy',lw=2,linestyle='--')

plt.xlabel('False Positive Rate')

plt.ylabel('True Positive Rate')

plt.title('2WEI ROC curve')

plt.legend(loc="lower right")

plt.show()

图9


2. E.tenella CDPK1 模型和 测试集预测

现在我们创建一个 E.tenella 模型来进行虚拟筛选和发展基于先前得分结果的用户自定义的得分功能。

2.1 建模

基于之前章节的对接结果, T.gondii 在全局和结合位点的序列中是更加接近 E.tenella 的.3T3U是 T.gondii 最优得分性能的一个结构(AUC和早期经验性AUC评价),因此我们使用SWISS-Model构建swissmodel3T3U.pdb

2.2 对接测试集

这里有提供22个化合物和详细实验数据作为模型的测试集

1.转化为SMILES格式


1

2

awk ' NR>1{print $2, $1} ' tdt2-challenge3-cdpk1_externaltestset.txt \

> testset.smi


2.生成构像


1

rdconf.py --maxconfs 1 testset.smi testset.sdf


3.对接测试集


1

2

smina --autobox_ligand allligs.pdb -r swissmodel3T3U.pdb  \

-l ../testset.sdf -o testset_docked.sdf --seed 0


4.选择和排序得分


1

2

sdsorter -sort minimizedAffinity testset_docked.sdf -reduceconfs 1 \

testset_docked_best.sdf -print -c > testset_default_ranking.txt


最佳得分构像如下图

![图7](/images/2017/smina/SMINA8.png)

2.3 用户自定义得分

我们发展了用户自定义得分功能用以 T.gondii 测试集参数化。增加到默认的得分功能上。

1.删除atom-types项目(如果考虑所有可能的原子类型组合,这将是一个非常大的项目),同样包括默认的第二Gaussian项目,他同样存在于默认得分之中。所有的权重均设为1.0


1

2

3

4

smina --print_terms | \

grep -v -E 'atom_type|constant|num_|ligand_length' > allterms

echo "gauss(o=3,_w=2,_c=8)" >> allterms

sed -i 's/^/1.0 /' allterms


grep -v -E 的用法比较少, -v 表示不包含匹配文本的所有行, -E 为启用正则表达式

2.得分最高的构像3T3U测试集拿来进行比较,由于我这里只做了2WEI,所以实际上此步我用的2WEI


1

2

3

4

5

6

7

8

# 提取最高得分构像

sdsorter -sort minimizedAffinity -reduceconfs 1 \

3T3U_docked.sdf.gz 3T3U_docked_single.sdf.gz

# 重新打分,不用再继续跑

smina --custom_soring allterms --score_only

-l 3T3U_docked_single.sdf.gz -r 3T3U.pdb | \

grep "##" | sed "s/##//" | \

awk '{print $1 ~ /_active/, $0}' > allscores


原文档中最后一行为/active/,但是个人觉得并没有达到区分active和unactive的用途,故本人对其进行了修改。可以从此下载文件。

3.作者使用R中的 rms 包和 logistical 回归来进行后向变量选择来契合活性数据的得分。为了防止过拟合,只选择p-值较低的功能(<.0001>


1

2

3

4

5

6

7

8

9

10

install.packages("rms")

library(rms)

scores=read.table("allscores",header=T)[c(-2)]   #remove Name  colum

colnames(scores)[1]="activity"

formula=reformulate(name(scores)[c(-1)],respinse="activity")

fit=lrm(formula,data=scores,x=T,y=T)

fit2=fastbw(fit,rule="p",sls=0.0001)

for(n in names(fit2$coefficients)){

cat(sprintf("%f %s\n",fit2$coefficients[[n]],n))

}


4.保存测定的系数coefficients以创建新的评分函数。因为项目名称中包含特征符,其在R中是无效命名,所以需要编辑输入以重新储存适当的terms,得分函数的结果如图下:




5.通过将系数乘以每一项的平均值来创建该评分函数中每个项的贡献的平均估计


1

2

3

4

aves=apply(scores[name(coef(fit2))[-1]],2,mean)*coef(fit2)[-1]

for(n in names(aves)){

cat(sprintf("%f %s\n",aves[[n]],n))

}


结果如下图

疏水性,VDW,氢键和溶剂化为有利相互作用,非疏水性,排斥性和受体-受体为负相关作用。VDW和gauss具有非常强的权重,线性氢键和Lennard Jones氢键展示了一些权重。

4*.Python包实现,Python的scikit-learn中对于logistical并没有P-值,译者对与R也不是特别理解,不知道能否领会作者的含义,以下仅供参考。 我们采用scikit-learn包中的logistical回归,同时将数据为训练集(75%)与测试集,检测是否过拟合,若没有测试集效果也不错,就使用所有的作为训练集再次训练。







请到「今天看啥」查看全文


推荐文章
彬彬有理  ·  身体出现这几个症状,要提早重视
20 小时前
宛央女子  ·  穿真丝内衣的福报来了
昨天
宛央女子  ·  用永远天真无畏的杉菜送别大S
4 天前
重庆共青团  ·  25岁徐枫灿,二级机长!
4 天前
重庆共青团  ·  25岁徐枫灿,二级机长!
4 天前
世纪风云潮  ·  七月间的国际大事
7 年前
华中科技大学  ·  小科要参加全运会了,你来不来?
7 年前