从零开始分子对接

从零开始分子对接

  记录一下学习分子对接的流程,希望给后人提供一点帮助。

工作环境

Autodock

  点击下载链接下载,该软件为纯命令行软件,图形界面需要使用MGLTools。此外有Autodock-GPU,速度更快,但仅限于Linux环境,且需配置CUDA环境,因为用不上在此不提供演示。

MGLTools

  点击下载链接下载,由于Autodock/MGLTools均需要一定专业知识,因此更推荐使用Autodock Vina。

Autodock Vina

  点击下载链接下载,目前一般都用这个对接。建议安装1.1.2,后续版本没有log,不方便抄代码。

Quick Vina

  Vina的进一步优化,下载链接,据说性能更高。

Quick Vina-W

  适用于盲对接(blind docking),也就是不知道配体-受体结合位点,下载链接同上,可以一块下载。

Openbabel

  以前写过,不再赘叙。

PyMol

  在有编译好的whl文件,使用pip安装即可。

Chimera

  用来展示结合蛋白分子以及口袋,比PyMol更复杂一点。在下载。

受体配体处理

受体准备

数据库下载

  PDB数据库进行检索,如果有完整的X射线晶体结构则可以下载使用,遵循原则是来自人类,分辨率≤2.50 Å。

同源模建

  如果没有完整结构,但有同源序列,且相似度大于30%,则可以进行同源模建。先访问Uniprot数据库,搜索受体,获取蛋白质序列。其原理是基于相似的氨基酸序列对应着相似的蛋白质结构。也可以进行叠加证明模板可靠性,简而言之用PyMol分别打开受体与模板的pdb文件,用align A, B命令即可获得叠加图,输出作为支持材料。

Swiss model法

  打开Swiss model网页,输入直接粘贴序列运行即可。对于生成结果,网站会自动排序,需要了解的指标为:

软件法

  用MOE等软件,有同源模建模块,因为暂时用不上,先不写了。有兴趣可以自己研究。

从头计算

  如果无同源序列,那就得考虑从头计算了。

I-TASSER

  打开I-TASSER,注册账户,按流程操作。建议参考同源建模方法整理

穿线法/综合法

  可以看山东大学的mooc,有个初步介绍。

AlphaFold

  基于AlphaFold模型生成的蛋白数据库,访问AlphaFold搜索下载。

受体评估

SAVES

  SAVES是目前最常用的综合评估网页,上传pdb文件就能自动开始5种评估。

PDBsum

  PDBsum也是一个经典评估网页,上传pdb文件提交。

相关概念解释:

  • PROCHECK用于评估模型三维结构的立体化学参数,结果可用Ramachandran图表示,横坐标为多肽蛋白质立体结构中肽键内α碳原子和羰基碳原子间的键的旋转度psi(ψ),纵坐标为α碳原子和氮原子间的键的旋转度phi(φ),理论上psi和phi都可自由变化,但键的转动会带动其它原子一起转动,因此受分子各个基团的空间障碍和作用力的影响,Ramachandran图中就出现了氨基酸的允许区域和不允许构象区域。主要看拉氏图检查Cα的两面角是否合理,合格的模型超过90%的残基都应该落在红色(允许区域)和正黄色(额外允许区域)落到其他区域的残基应当被查看并修正。

  • VERIFY_3D对同源建模所得模型的质量进行了可视化分析,主要是用来判断模型与其本身氨基酸序列之间的兼容性,对氨基酸残基数大于100的蛋白质,VERIFY_3D的评估结果更准确。超过80%的残基拥有大于0.2的3D/1D值,则模型质量合格,低于0.2的部分需要进一步修正。

  • PROSA的评价结果反映的是目的蛋白三维结构氨基酸残基之间的相互作用能,由公式计算得到:ZS,C=(ES,CS,C)/σS。其中ES,C表示序列S间的残基在C构象空间中的相互作用能,而ÊS,C指序列S间的残基在C构象空间的相互作用能的平均值,下面的σS为相对应的标准偏差。将含有目的蛋白空间结构的PDB格式的文件上传到网站便可得到PROSA的评价结果,计算所得数值ZS,C即为 Z-score,该值通常为负数。在PROSA结果图中还包括所有已经实验确定的现存于RCSB PDB数据库中与目的蛋白大小相近的蛋白质的Z-score值,当目标蛋白的Z-score值分布在这些已知蛋白的Z-score值绘制的图形范围内时,则证明目标蛋白结构的能量是合理的,反之,则认为结构是不合理的。

  • whatcheck提交的蛋白结构与正常结构之间的差异,指标贼多,绿色多就当通过了。

  • errat计算0.35 nm范围之内,不同的原子类型对之间形成的非键相互作用的数目(侧链)。得分>85较好,晶体可达到95,一般来说结果在91以内,很容易过。

  • prove与预先计算好的一系列标准体积的差别,用z-score来表示,显示模版蛋白质与待测蛋白之间的匹配程度,越高越好。

受体预处理

  主要就是去除原有小分子(ligand)、加氢(polar)、导出为pdbqt文件。大概操作如下:

  • PyMOL打开pdb文件,右下角点S显示序列,删掉配体和多余的水分子,Export molecular保存成为pdb文件;
  • 用MGLTools打开上面的pdb文件,在左侧的方框中勾选上protein,点击Edit—Hydrogens—Add,为蛋白质结构加氢原子;
  • 选择Grid—Macromolecule—Choose,点击选择protein作为AutoDock4的对接分子(Select Molecule),点击确定,保存为pdbqt文件。

口袋预测

在线预测

  使用DoGSiteScorer网站预测,上传pdb文件,选择go,再选择右边的DoGSiteScorer—Calculate,然后download文件。该文件为压缩包,其中pockets为预测的口袋,解压后用Chimera可以加载,通过叠加受体蛋白pdb,可以获得口袋与受体相对关系的图片。

设置口袋

  residues则是预测的残基位点,使用MGLTools的ADT打开选定口袋的pdb文件以及受体pdb文件,将残基的C点上,显示为球状氨基酸。基于此来设置口袋区域,Grid—Grid box,在Options中先把Spacing调到1,然后通过调整box大小和坐标让box正好包住球状区域。然后将参数复制到config.txt文件。可以看这个

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# box数据只是示例
# 需要自行修改受体配体文件名
receptor = model_01.pdbqt

center_x = 19.369
center_y = 11.437
center_z = -32.072
size_x = 32
size_y = 26
size_z = 26

exhaustiveness = 20
energy_range = 4
num_modes = 20

  实在不知道怎么设置,就把蛋白全部罩住吧,不清楚后果(逃。

配体准备

手动绘制

  高手可以用ChemDraw绘制,不赘述。

下载

  去Zinc等数据库下载小分子文件。

对接

  由于smile直接转mol2会遇到丢失结构的bug,血泪教训,说实话也不知道怎么解决,还好不是所有分子都这样。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
cd ~   # 进入home目录
mkdir docking # 新建工作目录

# 建立对接脚本
cd docking # 进入docking目录
mkdir log_file
chmod +x /home/catkin/mgltools_x86_64Linux2_1.5.6/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_ligand4.py # 让prepare_ligand4.py可以直接被执行
vim vina.sh # 建立shell脚本

#!bin/bash
#vina批量对接脚本
#本脚本及配体受体文件在~/docking 目录下面,prepare_ligand4.py和vina使用时根据情况自行调整。

for m in mol_*.sdf; do
b=`basename $m .sdf`
echo Processing ligand $b
mkdir -p $b
obabel -i sdf $m -o mol2 -O ${b}.mol2 -h --gen3d
/home/catkin/mgltools_x86_64Linux2_1.5.6/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_ligand4.py -l ${b}.mol2 -o ${b}.pdbqt
/home/catkin/autodock_vina_1_1_2_linux_x86/bin/vina --config config.txt --ligand ${b}.pdbqt --out ${b}/${b}_out.pdbqt --log log_file/${b}.log
done

数据处理

提取亲和力数据

  基于Becky的脚本修改,需要自己调整一下路径和正则。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
"""
# @Time : Wed Sep 11 11:05:47 2024
# @Author : Catkin
# @Website : blog.catkin.moe
"""
import re
import glob

library = 'FDA'

log_file = glob.glob('..\\Docking\\Docking_results\\' + '*.log')
file1 = open('..\\Docking\\Docking_results\\' + library + '.csv',"w")
for i in log_file:
file2 = open(i,'r')
data = file2.readlines()
file2.close()
patterns_x = re.compile('mol_(.*?).log')
x = re.findall(patterns_x,i)[0]
for j in data:
if len(j.split())>=2 and j.split()[0] == "1":
affinity = str(x) + "," + j.split()[1]
file1.writelines(affinity + '\n')
file1.close()
构建蛋白/配体复合物

  运行ADT,选择Analyze—Docking—Open autoDock vina result—OK,左右切换不同构象,选择最左边的构象,File—Save—Write PDB—OK。随后运行PyMol,分别载入刚才的分子和蛋白。然后File—Export Molecule—PDB options—Save,选择PDB格式保存。再重新打开复合物文件,选择配体,然后A—find—halogen-bond interactions。

检查是否正常形成氢键

  运行ADT,选择Analyze—Docking—Open autoDock vina result—OK,选择Analyze—Macromolecule打开受体文件,再Analyze—Docking—Show interactions,选择受体的C隐藏球状文件以防干扰,取消display msms,close contact,勾选hydrogen bonds。

绘制蛋白/小分子互作图
2D

  参考视频,主要可以用的工具如下:

  以Proteins Plus的2D作图为例,访问页面,上传复合结构,看到分子,点create pocket,右边PoseView,输入Ligand id(在左边点Ligands可以查看),然后会出现2D相互作用图,可以下载svg格式保存。

  如果是DS Visualizer,File—Open选择配体复合结构,然后选择Receptor-Ligand Interactions,再选择Show 2D Diagram,然后可以保存为PNG格式。

3D

  参考视频

  1. 三维结构需要先把复合结构用Pymol打开,在右下角切换成Chains,可以分别选出蛋白质和小分子并命名,为了方便区分给蛋白质换个颜色;
  2. 点击小分子,A—find—polar contacts—to other atoms in object显示氢键;
  3. 点蛋白质,S—licorice—sticks,显示球棍模式,这时候就能方便选中残基,右下角切换回Residues,选择氢键连接的残基,右键命名,然后选回蛋白质H—licorice—sticks,把球棍模式隐藏;
  4. 选择刚才的残基,显示回球棍模式,为了方便区分再换个颜色;
  5. 主屏幕右键还原最初位置,然后右上角Draw输出全视角图片;
  6. 选择蛋白质,S—surface显示表面,Setting—Transparency—Surface设置一个透明度,H—cartoon蛋白质结构,右上角Draw输出第二张图片;
  7. 还原至第5步,放大局部区域,使氢键和残基可以很好展示,再选择氢键,S—labels,选择残基,L—residues,再Setting—Label—Size调一下字体大小,右下角把Viewing切成Editing,按住Ctrl可以拖动Label,调整好输出图片。

参考链接

[新版]AutoDock Vina分子对接小白1小时入门教程

使用SWISS-MODEL的同源建模及模型检验

同源建模方法整理

Linux系统下,使用open babel、MGLTools、AutoDock Vina批量筛选ZINC数据库的方法记录

【GROMACS】蛋白-配体动力学模拟之一:准备拓扑文件

分子对接 - 如何生成蛋白和小分子的2D相互作用图片?

……

评论