新建
从蛋白/小分子准备开始的Uni-Dock分子对接全流程介绍
yuanyn@dp.tech
推荐镜像 :docking-pipeline:v0.4.4.4
推荐机型 :c12_m46_1 * NVIDIA GPU B
赞 3
6
5
目录
从蛋白/小分子准备开始的Uni-Dock分子对接全流程介绍
代码
文本
代码
文本
©️ Copyright 2023 @ Authors
作者：
袁彦南 📨
日期：2023-07-10
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的 开始连接 按钮，选择 docking-pipeline:v0.0.4.4镜像和任意GPU机型即可开始。
代码
文本
背景
代码
文本
分子对接是通过计算模拟快速评估候选分子与目标靶点的结合能力的技术，被广泛应用于药物设计的早期阶段虚拟筛选中，用来在大规模分子库中寻找潜在的活性分子[1]。
Uni-Dock[2]是深势科技发布的基于GPU加速的分子对接引擎，可以实现在NVIDIA V100 GPU上计算速度对比单核CPU超过1600倍的加速比的极致加速。
关于Uni-Dock的更具体的介绍可见文章：https://dptechnology.feishu.cn/docx/XJZ8dE43Uo4EJUxHdZmcmHHhnAg
Uni-Dock的基本使用方法已在之前的notebook中有所介绍：https://nb.bohrium.dp.tech/detail/1288
代码
文本
本教程将从最常用的PDB格式的蛋白文件和SDF格式的小分子文件出发，通过蛋白和配体准备对蛋白和小分子进行标准化并转化为Uni-Dock所需的输入形式，利用Uni-Dock的 SDF-IN-SDF-OUT 模式来对处理过的满足Uni-Dock格式需求的SDF输入文件进行分子对接，输出SDF格式的对接结果文件，并在之后进行后处理来修复小分子坐标并获取对接结果。
代码
文本
0. 准备
代码
文本
0-1. 环境
代码
文本
Docking流程使用的环境较为复杂，推荐在notebook中直接连接使用准备好的docking-pipeline:v0.0.1-f5镜像。
代码
文本
0-2. 数据
代码
文本
本流程可以使用任意PDB格式蛋白文件和SDF格式小分子文件，本教程提供了样例文件来使用，你也可以把下面的蛋白和小分子文件URL换成自己的文件URL来完成本流程。
代码
文本
[1]
from typing import List
from pathlib import Path
import os
import traceback
import shutil
import re
import glob
import uuid
import math
from pprint import pprint
import subprocess
# 工作目录
work_dir = f'/tmp/{uuid.uuid4().hex}'
os.makedirs(work_dir, exist_ok=True)
os.environ["WORK_DIR"] = work_dir
# 蛋白文件
! cd $WORK_DIR; wget https://app.bohrium.dp.tech/download/artifacts/applications/uni-dock/artifacts/receptor.pdb
# 小分子文件
! cd $WORK_DIR; wget https://app.bohrium.dp.tech/download/artifacts/applications/uni-dock/artifacts/actives_final.sdf
receptor_file = os.path.join(work_dir, "receptor.pdb")
ligand_file = os.path.join(work_dir, "actives_final.sdf")
--2023-08-02 18:03:01-- https://app.bohrium.dp.tech/download/artifacts/applications/uni-dock/artifacts/receptor.pdb Resolving ga.dp.tech (ga.dp.tech)... 10.255.254.37, 10.255.254.7, 10.255.254.18 Connecting to ga.dp.tech (ga.dp.tech)|10.255.254.37|:8118... connected. Proxy request sent, awaiting response... 200 OK Length: 88385 (86K) [application/vnd.palm] Saving to: ‘receptor.pdb’ receptor.pdb 100%[===================>] 86.31K --.-KB/s in 0.01s 2023-08-02 18:03:01 (7.16 MB/s) - ‘receptor.pdb’ saved [88385/88385] --2023-08-02 18:03:01-- https://app.bohrium.dp.tech/download/artifacts/applications/uni-dock/artifacts/actives_final.sdf Resolving ga.dp.tech (ga.dp.tech)... 10.255.254.37, 10.255.254.18, 10.255.254.7 Connecting to ga.dp.tech (ga.dp.tech)|10.255.254.37|:8118... connected. Proxy request sent, awaiting response... 200 OK Length: 783390 (765K) [application/vnd.Kinar] Saving to: ‘actives_final.sdf’ actives_final.sdf 100%[===================>] 765.03K --.-KB/s in 0.1s 2023-08-02 18:03:02 (5.90 MB/s) - ‘actives_final.sdf’ saved [783390/783390]
代码
文本
1. 蛋白准备
代码
文本
使用 MGLTools 包中的 prepare_receptor4.py 脚本将输入的PDB格式文件转为Uni-Dock输入需要的PDBQT格式文件，由于该脚本处理金属原子存在问题，需提前对蛋白文件中的金属进行手动处理。
代码
文本
[2]
def pdb2pdbqt(protein_file:str, output_file:str):
env_path = "/opt/mamba/envs/mgltools"
script_path = os.path.join(env_path, "MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py")
pythonsh_path = f'{env_path}/bin/pythonsh'
with open(protein_file, "r") as f:
protein_lines = f.readlines()
noh_protein_file = f'/tmp/{uuid.uuid4().hex}_tmp.pdb'
with open(noh_protein_file, "w") as f:
for line in protein_lines:
if "HN" not in line[12:16].strip():
if 'ZN' in line[12:16]:
line = line[:12] + "ZN " + line[16:]
elif 'MG' in line[12:16]:
line = line[:12] + "MG " + line[16:]
f.write(line)
args = [
pythonsh_path,
script_path,
'-r', noh_protein_file,
'-o', output_file,
'-U', 'nphs_lps_nonstdres'
]
print("PDB2PDBQT args: ")
pprint(args)
response = subprocess.run(args,
stdout=subprocess.PIPE, stderr=subprocess.PIPE,
encoding="utf-8")
print(response.stdout)
if response.stderr:
print(response.stderr)
if not os.path.exists(output_file):
raise KeyError("Pdb2Pdbqt err")
Path(noh_protein_file).unlink(missing_ok=True)
processed_receptor = receptor_file+"qt"
pdb2pdbqt(receptor_file, processed_receptor)
PDB2PDBQT args: ['/opt/mamba/envs/mgltools/bin/pythonsh', '/opt/mamba/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py', '-r', '/tmp/b0246b50fb4e46dc88aaa1f09efd79cd_tmp.pdb', '-o', '/tmp/3519d8ca22ee45a98dd11c895ca1a5b0/receptor.pdbqt', '-U', 'nphs_lps_nonstdres'] Sorry, there are no Gasteiger parameters available for atom b0246b50fb4e46dc88aaa1f09efd79cd_tmp: :ZN4490:ZN
代码
文本
2. 配体准备
代码
文本
使用 Uni-LigandPrep 包对小分子进行修复，加氢，计算可旋转键/原子类型信息（用于作为Uni-Dock输入）操作，得到处理后的小分子文件。
代码
文本
[3]
def ligprep(ligand_files:List[str], output_dir:str) -> List[str]:
from uniligprep.IO import read_ligand
from uniligprep.utils.mol_utils import set_properties
from uniligprep.PropCalculation.pdbqt_info_obabel import compute_pdbqt_info_obabel
from uniligprep.Protonation.AddHWithProperty import add_hydrogen_mol
from uniligprep.format_converter.convert_fmt import save_ligand_file
os.makedirs(output_dir, exist_ok=True)
for ligand_file in ligand_files:
ligand_name = os.path.splitext(re.sub(r"[\s]", "", os.path.basename(ligand_file)))[0]
mols = read_ligand(ligand_file)
print(len(mols))
for i, mol in enumerate(mols):
try:
mol = add_hydrogen_mol(mol)
pdbqt_info = compute_pdbqt_info_obabel(mol.get_mol())
set_properties(mol.get_mol(), pdbqt_info)
save_ligand_file(mol, os.path.join(output_dir, f'{ligand_name}_{i}.sdf' if len(mols)>1 else f'{ligand_name}.sdf'))
except:
print(f'Ligand {ligand_name} {i} prepared err: {traceback.format_exc()}')
return glob.glob(os.path.join(output_dir, "*.sdf"))
prepared_ligands_dir = os.path.join(work_dir, "prepared_ligands")
os.makedirs(prepared_ligands_dir, exist_ok=True)
prepared_ligands = ligprep([ligand_file], prepared_ligands_dir)
161
代码
文本
3. 分子对接
代码
文本
3-1. 设置参数
代码
文本
这里用一个字典来定义Uni-Dock所需的对接参数，包括打分函数，模式，口袋，生成构象数量等，该字典的key为Uni-Dock参数名，value为参数值。
代码
文本
[4]
params = dict()
代码
文本
设置口袋大小：
代码
文本
[5]
params.update({
"center_x": -20,
"center_y": 30,
"center_z": 60,
"size_x": 20,
"size_y": 20,
"size_z": 20,
})
代码
文本
设置打分函数（vina/vinardo/ad4）:
代码
文本
[6]
params["scoring"] = "vina"
代码
文本
设置mode（exhaustiveness和max_step参数的结合，fast/balance/detail）:
代码
文本
[7]
params["search_mode"] = "fast"
代码
文本
设置对接输出pose数量：
代码
文本
[8]
params["num_modes"] = 1
代码
文本
3-2(可选). ad4方法前处理
代码
文本
如果对接参数中的打分函数（scoring）选择了ad4，那么需要对蛋白文件进行提前处理，生成**.map**文件，并在之后使用Uni-Dock时以这些生成的map文件作为输入。
代码
文本
[27]
params["scoring"] = "ad4"
代码
文本
[28]
def find_atom_types(ligand_files:List[str]):
AD_atm_types = ["H", "HD", "HS", "C", "A", "N", "NA", "NS", "OA", "OS", "F", "Mg", "MG", "P", "SA", "S", "Cl", "CL",
"Ca", "CA", "Mn", "MN", "Fe", "FE", "Zn", "ZN", "Br", "BR", "I", "Z", "G", "GA", "J", "Q", "O"]
atms = set()
tag = False
for ligand_file in ligand_files:
with open(ligand_file, "r") as f:
for line in f.readlines():
if line.strip():
if line.startswith("> <atomInfo>"):
tag = True
elif tag and (line.startswith("> <") or line.startswith("$$$$")):
tag = False
elif tag:
atms.add(line[13:].strip())
atms = [a for a in atms if a in AD_atm_types]
return atms
def gen_ad4_map(protein_file:str, ligand_files:List[str], mapdir:str, docking_params:dict) -> str:
os.makedirs(os.path.abspath(mapdir), exist_ok=True)
env_path = "/opt/mamba/envs/mgltools"
script_path = os.path.join(env_path, "MGLToolsPckgs/AutoDockTools/Utilities24/prepare_gpf4.py")
autogrid_dat_file = "/opt/data/AD4.1_bound.dat"
atom_types = find_atom_types(ligand_files)
name = os.path.splitext(os.path.basename(protein_file))[0]
shutil.copyfile(protein_file, os.path.join(mapdir, os.path.basename(protein_file)))
center = [docking_params["center_x"], docking_params["center_y"], docking_params["center_z"]]
box_size = [docking_params["size_x"], docking_params["size_y"], docking_params["size_z"]]
spacing = 0.375
npts = [math.ceil(s/spacing) for s in box_size]
cmd = []
cmd.append(f"{env_path}/bin/pythonsh {script_path} -r {os.path.basename(name)}.pdbqt \
-p gridcenter='{center[0]},{center[1]},{center[2]}' -p npts='{npts[0]},{npts[1]},{npts[2]}' \
-p spacing={spacing} -p ligand_types={','.join(atom_types)} -o {name}.gpf")
cmd.append(f"sed -i '1i parameter_file {autogrid_dat_file}' {name}.gpf")
cmd.append(f"{env_path}/bin/autogrid4 -p {name}.gpf -l {name}.glg")
print("AD4 map cmd:")
pprint(cmd)
response = subprocess.run("\n".join(cmd), shell=True,
stdout=subprocess.PIPE, stderr=subprocess.PIPE,
encoding="utf-8", cwd=mapdir)
print(response.stdout)
if response.stderr:
print(response.stderr)
return mapdir + "/" + name
代码
文本
3-3. 使用Uni-Dock对接
代码
文本
给定蛋白，小分子文件和对接参数，我们这里将使用Uni-Dock进行分子对接。Uni-Dock支持SDF-IN-SDF-OUT，如果小分子SDF文件满足特定的格式，即可以以这些SDF文件直接作为输入，而不用转成PDBQT格式，并且输出SDF格式文件的结果，这样能够保留小分子文件的原始键级，方便后续流程使用。
代码
文本
Uni-Dock要求的SDF格式如下：
代码
文本
> <fragInfo> (1)
<atom number 1> <atom number 2> ...
> <torsionInfo> (1)
<torsion_start_atom_number> <torsion_end_atom_number> <start_atom_frag_ind> <end_atom_frag_ind>
> <atomInfo> (1)
<atom number> <formal charge> <atom type>
代码
文本
代码
文本
在上面步骤2配体准备中，我们已完成了对小分子的标准化操作并且生成了Uni-Dock所需的可旋转键和原子信息，Uni-Dock可以直接以这些SDF文件作为输入文件并输出SDF格式的结果文件，这样可以保留小分子的完整键级信息，方便用于后续的MM PB/GBSA和FEP过程。
代码
文本
[9]
def run_docking(receptor_file:str, ligand_files:List[str], results_dir:str, params:dict):
os.makedirs(results_dir, exist_ok=True)
scoring = params["scoring"]
docking_args = ["unidock"]
if scoring == "ad4":
mapdir = os.path.join(results_dir, "map")
os.makedirs(mapdir, exist_ok=True)
mapdir = gen_ad4_map(receptor_file, ligand_files, mapdir, params)
params["maps"] = mapdir
else:
params["receptor"] = receptor_file
file_list_path = f'/tmp/{uuid.uuid4().hex}_file_list'
with open(file_list_path, "w") as f:
f.write("\n".join(ligand_files))
params["ligand_index"] = file_list_path
params["dir"] = results_dir
for k, v in params.items():
docking_args.append(f'--{k}')
docking_args.append(f'{v}')
docking_args.extend([
"--verbosity",
"0",
"--refine_step",
"3",
"--keep_nonpolar_H"
])
docking_cmd = " ".join(docking_args)
print(f"Docking cmd:{docking_cmd}")
response = subprocess.run(docking_cmd, shell=True,
stdout=subprocess.PIPE, stderr=subprocess.PIPE,
encoding="utf-8")
print(response.stdout)
if response.stderr:
print(response.stderr)
Path(file_list_path).unlink(missing_ok=True)
results_dir = os.path.join(work_dir, "results")
run_docking(processed_receptor, prepared_ligands, results_dir, params)
Docking cmd:unidock --center_x -20 --center_y 30 --center_z 60 --size_x 20 --size_y 20 --size_z 20 --scoring vina --search_mode fast --num_modes 1 --receptor /tmp/3519d8ca22ee45a98dd11c895ca1a5b0/receptor.pdbqt --ligand_index /tmp/113bfe18a3eb42848f509bdd8e02a337_file_list --dir /tmp/3519d8ca22ee45a98dd11c895ca1a5b0/results --verbosity 0 --refine_step 3 --keep_nonpolar_H Uni-Dock v1.2.3 entering done exiting done Total ligands: 161 Avaliable Memory = 10830MiB Total Memory = 11011MiB Batch 1 size: 161 Kernel running time: 15 entering done exiting done Batch 1 running time: 21249ms
代码
文本
4. 后处理/结果分析
代码
文本
4-1. 非极性氢后处理
代码
文本
由于Uni-Dock使用SDF-IN-SDF-OUT模式时会忽略小分子的非极性氢，因此对接后的结果小分子文件非极性氢位置是不正确的，如图：
代码
文本
[10]
one_result_ligand = glob.glob(os.path.join(results_dir, "*.sdf"))[0]
print(one_result_ligand)
import nglview as nv
nv.show_file(one_result_ligand)
/tmp/3519d8ca22ee45a98dd11c895ca1a5b0/results/actives_final_59_out.sdf
代码
文本
这里需要使用Uni-LigandPrep的ReProto功能对小分子非极性氢进行重新生成：
代码
文本
[13]
def regen_hydrogen(results_dir:str):
from rdkit import Chem
from uniligprep.utils.mol_utils import set_properties
from uniligprep.IO import read_ligand
from uniligprep.Protonation.reprotonation import reprotonation
from uniligprep.format_converter.convert_fmt import get_mol_contents_by_fmt
ligand_files = glob.glob(os.path.join(results_dir, "*.sdf"))
for ligand_file in ligand_files:
try:
mols = read_ligand(ligand_file)
reproto_rdkit_mols = list()
for mol in mols:
reproto_rdkit_mol = reprotonation(mol.get_mol())
set_properties(reproto_rdkit_mol, mol.get_props())
reproto_rdkit_mols.append(reproto_rdkit_mol)
with Chem.SDWriter(ligand_file) as writer:
for reproto_rdkit_mol in reproto_rdkit_mols:
writer.write(reproto_rdkit_mol)
except:
print(f'Docking result ligand {ligand_file} reproto failed: {traceback.format_exc()}')
regen_hydrogen(results_dir)
代码
文本
重生成后的结果：
代码
文本
[15]
nv.show_file(one_result_ligand)
代码
文本
4-2. 获取结果打分
代码
文本
Uni-Dock打分会写在结果小分子文件里，从文件中读出打分并排序就可挑选出本次对接的最优结果。
代码
文本
[16]
def analysis(docking_results_dir:str, table_path:str) -> object:
import pandas as pd
from rdkit import Chem
df = pd.DataFrame(columns=["name", "pose", "docking_score"])
result_files = glob.glob(os.path.join(docking_results_dir, "*.sdf"))
for result_file in result_files:
mols = Chem.SDMolSupplier(result_file, removeHs=False, sanitize=True)
for i, mol in enumerate(mols):
if mol:
score_line = mol.GetProp("Uni-Dock RESULT")
score = float(score_line.partition("LOWER_BOUND=")[0].lstrip("ENERGY="))
df.loc[df.shape[0]] = [os.path.basename(result_file), i+1, score]
df = df.sort_values(by=["docking_score"], ignore_index=True)
df.to_csv(table_path, index=None)
return df
table_path = f'/tmp/{uuid.uuid4().hex}_score.csv'
res_table = analysis(results_dir, table_path)
print(res_table.head())
name pose docking_score 0 actives_final_41_out.sdf 1 -5.407 1 actives_final_156_out.sdf 1 -5.257 2 actives_final_27_out.sdf 1 -5.166 3 actives_final_39_out.sdf 1 -5.033 4 actives_final_120_out.sdf 1 -5.010
代码
文本
[17]
view = nv.show_file(receptor_file)
view.add_component(one_result_ligand)
view
代码
文本
参考
- https://dptechnology.feishu.cn/docx/XJZ8dE43Uo4EJUxHdZmcmHHhnAg
- Yu, Y., Cai, C., Wang, J., Bo, Z., Zhu, Z., & Zheng, H. (2023). Uni-Dock: GPU-Accelerated Docking Enables Ultralarge Virtual Screening. Journal of Chemical Theory and Computation.
代码
文本
已赞3
本文被以下合集收录
App related
Charmy Niu
更新于 2024-01-17
10 篇3 人关注
分子对接
18895326035@163.com
更新于 2024-05-21
6 篇0 人关注
推荐阅读
公开
Hands-on to APEX (v1.2) on Bohrium
zhuoyli@connect.hku.hk
更新于 2024-08-08
4 赞5 转存文件
公开
OLED量化属性预测：Top4-方案
SeeFun4Sci
发布于 2023-10-31
3 赞2 转存文件