摘要
VirtualFlow 是一款开源的、基于结构的、超高通量虚拟筛选平台，Gorulla 等人成功地用8000个核心花两周时间完成对10多亿化合物进行虚拟筛选，发现了 nM 水平的 KEAP1 抑制剂。本文演示了如何在 Cloudam 上使用 VirtualFlow 进行此类超高通量虚拟筛选，以 VirtualFlow 自带的算例 VFVS_GK 演示了：1）如何准备输入文件；2）编辑配置文件；3）规划计算任务；4）执行虚拟筛选。
关于 VirtualFlow
基于结构的虚拟筛选是发现先导化合物的重要途经，其发现的苗头化合物的质量随着筛选的化合物数量增加而提高①。虽然已经有很大的化合物库，但是进行这么大规模基于结构虚拟筛选的计算能力一直是个问题，比如进行此类计算集群的可获得性、高效性与灵活性一直是个挑战。
VirtulFlow 是哈佛医学院 Gorulla 等人②③开发的开源超高通量虚拟筛选平台，它可以让你以巨量并行的方式在计算集群上完成超大规模数据库的虚拟筛选任务。据报道 NAF2-KEAP1 通路在氧化应激与炎症中起到关键的作用与多种疾病相关，目前已有10个靶向 KEAP1 的药物处于临床研究阶段，9个药物处于临床前研究阶段。在 Gorulla 等人4的研究中，用 VirutalFlow 将13亿的 Real 化合物库与3.3亿的 ZINC 化合物库对接到 KEAP1-NAF2 相互作用界面的 KEAP1 结合位点，识别出一系列结构多样与 KEAP1 结合、阻断 NAF2-KEAP1 相互作用的潜在化合物。生物实验测试结果表明，其中一个化合物与 KEAP1 的结合亲和力达到 nM 水平（Kd=144nM）并且阻断 NAF2-KEAP1 相互作用。
Figure 1. VirtualFlow的虚拟筛选流程
Gorulla 的虚拟筛选分两步完成：1）对常规蛋白刚性、配体柔性的策略将13亿 Real 数据库以及3.3亿 ZINC 数据库对接到 Keap1，取打分最高的300万化合物。本步骤在8000核心的机器上花了4周时间完成。2）用蛋白柔性、配体柔性策略，将上一步打分最高的300万化合物进一步打分，将 KEAP1 与 NAF2 相互界面的13残基考虑为构象柔性，用 Smina Vinardo 与 AutoDock Vina 进行计算，分别重复两次。对第一步打分最高的98个化合物以及第二步打分最高500个化合物进行了实验验证。
Figure 2. KEAP1的虚拟筛选策略
上述虚拟筛选采用 QVINA2 分子对接软件进行初筛、再用 SMINA 与 VINA 进行两次的蛋白柔性对接打分两次可以通过 VirtualFLow 的参数文件简单的设置而实现。本文针对这些参数的设置进行演示。
关于云E算力平台
Cloudam 的云E算力平台整合了全球主流公有云近50个地域的高性能计算资源，开箱即用，无需硬件投入、无需运营维护、无需任务排队，支持一键式提交作业，具备自主学习与深度学习能力，能够根据用户的实例类型与计算要求智能推荐匹配合适的计算机型，将计算的虚拟损耗降至最低，提升计算效率，降低计算成本。
虽然 VirtualFlow 是开源的，但是将该代码迁移到自己的云服务器上还需要仔细配置与调试、需要一定的学习成本。本着降低学习成本、“开箱即用”的客户服务原则，Cloudam 的工程师们预先在云端配置好了软件，仅需对少数与项目相关的参数进行设置即可对13亿的 Real 数据库与3.3亿的 ZINC 数据库进行基于结构的虚拟筛选。本文将重点演示如何在 Cloudam上 用 VirtualFlow 进行基于结构的虚拟筛选。
详细操作步骤
VirtualFlow的文档与教程可以从它的主页获得：https://virtual-flow.org。本文的主要目的是演示如何在Cloudam上使用它。
1. 准备工作
我们假设：（1）你计划用VirtualFlow自带的13亿Real数据库与3.3亿ZINC数据库进行虚拟筛选；（2）你已经准备好了PDBQT格式的对接用的蛋白结构文件；（3）你计划用QVINA2进行第一轮的蛋白刚性，配体柔性分子对接；第二轮用AutoDock VINA进行蛋白柔性的分子对接。
2. 登录管理节点
你首先要在cloudam上创建账号，并开启管理节点，然后webssh登录管理节点。
3. 准备输入文件
以VirtualFlow自带的算例VFVS_GS为例，包含输入文件都包含在该文件夹里：
\---VFVS_GK| todo.all| vf.conf| +---input-files| \---ligand-library| | README.md| | | +---GA| | GACAAD.tar| | GACAAE.tar| | GACAAF.tar| | | \---HA| HACABE.tar| HACABF.tar| HACACE.tar| +---receptors| 4no7_prot.pdbqt| \---scenarios+---qvina02_rigid_receptor1| config.txt| \---smina_rigid_receptor1| config.txt
input-files
input-files/ligand-library目录存放事先准备好的配体数据库文件。VirtualFlow采用AutoDock Vina系列软件做为分子对接计算程序，此处为保存已经准备好对接用的PDBQT格式的化合物数据库。VirtualFlow提供了两个已经准备好的数据库可以直接使用：1）含14亿化合物的Real数据库；2）含12亿化合物的ZINC15数据库。
对于自己尚未准备过的化合物库，可以用VFLP来准备，或者按照
VFVS_GK/input-files/ligand-library的格式自行准备。
receptors
准备好的、可直接用于对接计算的蛋白结构文件，比如4no7_prot.pdbqt。
scenarios
提供给不同对接软件用的配置文件，比如
qvina02_rigid_rceptor01/config.txt,其内容如下:
receptor = ../input-files/receptors/4no7_prot.pdbqtcenter_x = -8.654center_y = 2.229center_z = 19.715size_x = 24.0size_y = 26.25size_z = 22.5exhaustiveness = 8cpu = 1
todo.all
是待对接计算的化合物列表。在
input-files/ligand-library里存储的化合物有26多亿，但是这些化合物只是准备好可以用来计算，而不必全部都需要被对接计算。只有被列在todo.all里的那些化合物才被计算。以前6行为例：
GACAAD_00000 9GACAAF_00000 22GACABD_00000 27GACABF_00000 48GACACC_00000 1GACACD_00000 30
VirtualFlow 按化学空间位置管理化合物，空间中相邻的几个化合物被打包压缩在一个文件，第一列对应该压缩文件的名称，第二列对应其中的化合物数。我们需要第二列的数值的总和，来规划后面每节点需处理的化合物数。
vf.conf
vf.conf是Virtualflow的配置文件，是最重要的文件，解释如下：
# Values have to be separated by colons ":" and without spaces, e.g: docking_scenario_names=vina-rigid:smina-rigid# Used for instance for the folder names in which the output files are stored# Should not be changed during runtime, and be the same for all joblines# Settable via range control files: Yes# 对接场景名称,可以任意命名。出于管理方便，建议与对接场景scenario名称一致# 计算结果输出目录将以此命名，可以进行一个或多个对接计算场景，使用":"隔离docking_scenario_names=qvina02_rigid_receptor1:smina_rigid_receptor1
# Relative path wrt to the tools folders# In each input folder must be the file config.txt which is used by the docking program to specify its options# If other input files are required by the docking type, usually specified in the config.txt file, they have to be in the same folder# Should not be changed during runtime, and be the same for all joblines# Settable via range control files: Yes# 配置文件config.txt所在目录，每个对接计算需要一个, ../input-files/为固定位置，不可修改，仅需将目录scenarios下的子目录填入即可docking_scenario_inputfolders=../input-files/qvina02_rigid_receptor1:../input-files/smina_rigid_receptor1
# Possible values: qvina02, qvina_w, vina, smina_rigid, smina_flexible, adfr# Values have to be separated by colons ":" and without spaces, e.g: docking_scenario_programs=vina:smina# The same programs can be used multiple times# Should not be changed during runtime, and be the same for all joblines# Settable via range control files: Yes# 对接计算所用的软件，合法的值为:qvina02, qvina_w, vina, smina_rigid, smina_flexible, adfrdocking_scenario_programs=qvina02:smina_rigid
# Series of integers separated by colons ":"# The number of values has to equal the number of docking programs specified in the variable "docking_programs"# The values are in the same order as the docking programs specified in the variable "docking_scenario_programs# e.g.: docking_scenario_replicas=1:1# possible range: 1-99999 per field/docking program# The docking scenario is comprised of all the docking types and their replicas# Should not be changed during runtime, and be the same for all joblines# Settable via range control files: Yes# 每个ligand需要被重复计算的次数docking_scenario_replicas=1:1
# When the folders are initially prepared the first time, the central todo list will be split into pieces of size . One task corresponds to one collection. # Recommended value: < 100000, e.g. 10000# Possible values: Positive integer# The smaller the value, the faster the ligand collections can be distributed.# For many types of clusters it is recommended if the total number of splitted todo lists stays below 10000.# Settable via range control files: Yes# todo.all文件中按找多少行数进行文件分割处理，建议按照分割后的文件数量在2000一下来设置此参数central_todo_list_splitting_size=1000
# Used as a limit of ligands for the to-do lists# This value should be divisible by the next setting "ligands_todo_per_refilling_step"# Settable via range control files: Yes# 每个处理进程需要处理的ligand的数量，注：总queque数=申请的总核心数， todo.all中累加为总的待处理ligand数量# 可以按照每个ligand处理需要的时长 x(乘) ligands_todo_per_queue得出每个轮次需要的时长# 可以理解为每个核心需要计算的化合物数，全部计算完保存起来，如果计算失败，那么这些化合物没有计算输出ligands_todo_per_queue=100
# The to-do files of the queues are filled with ligands per refill step # A number roughly equal to the average of number of ligands per collection is recommended# Settable via range control files: Yes# 单次计算时填充数量# ligands_todo_per_queue数量的化合物，分几批送进去计算，每次送的化合物数ligands_per_refilling_step=10
# Not (yet) available for LSF and SGE (is always set to 1)# Should not be changed during runtime, and be the same for all joblines# Settable via range control files: Yes# 应为申请的总节点数steps_per_job=2500
# Sets the slurm cpus-per-task variable (task = step) in SLURM# In LSF this corresponds to the number of slots per node# Should not be changed during runtime, and be the same for all joblines# Not yet available for SGE (always set to 1)# Settable via range control files: Yes# 应为每节点的核心数cpus_per_step=4
#steps_per_job * cpus_per_step = 使用的总节点数#steps_per_job=625， cpus_per_step=16，表示启动16核的机器共625台，总共10000核4. 虚拟筛选
为了方便计算，我们准备了 creatworkflow.sh 脚本，将该脚本复制到 home 目录，并按需要修改前面两行：
DOCKERING_CONFIG_DIR=/home/cloudam/VFVS_GKDOCKERING_WORKING_DIR=/home/cloudam/myfirsttest
其中第1行是我们刚刚准备的输入文件目录；第2行是进行虚拟筛选的工作目录，你可以任意命名，这里我命名为 myfirsttest 现在可以运行该脚本开始虚拟筛选：
./creatworkflow.sh
文献
Lyu, J.; Wang, S.; Balius, T. E.; Singh, I.; Levit, A.; Moroz, Y. S.; O’Meara, M. J.; Che, T.; Algaa, E.; Tolmachova, K.; Tolmachev, A. A.; Shoichet, B. K.; Roth, B. L.; Irwin, J. J. Ultra-Large Library Docking for Discovering New Chemotypes. Nature 2019, 566 (7743), 224–229. https://doi.org/10.1038/s41586-019-0917-9.
Gorgulla, C. VirtualFLow https://github.com/VirtualFlow (accessed Jul 19, 2020).
Gorgulla, C. VirtualFlow. https://www.virtual-flow.org (accessed Jul 19, 2020).
Gorgulla, C.; Boeszoermenyi, A.; Wang, Z.-F.; Fischer, P. D.; Coote, P. W.; Padmanabha Das, K. M.; Malets, Y. S.; Radchenko, D. S.; Moroz, Y. S.; Scott, D. A.; Fackeldey, K.; Hoffmann, M.; Iavniuk, I.; Wagner, G.; Arthanari, H. An Open-Source Drug Discovery Platform Enables Ultra-Large Virtual Screens. Nature 2020, 580 (7805), 663–668. https://doi.org/10.1038/s41586-020-2117-z.
本文转载自广州市墨灵格信息科技有限公司，由云端及墨灵格共同撰写
