RosettaLigand

RosettaLigand

前言

多巴胺是一种重要的神经递质，它通过五种多巴胺受体亚型发挥作用，一类g蛋白偶联受体(GPCRs)的重要成员。亚型二(D2R)和亚型二3 (D3R)通过抑制腺苷酸环化酶发挥作用，调节这两种受体有临床意义在精神分裂症治疗中的应用。但是，D2R之间的结合位点守恒程度较高D3R使得产生选择性结合其中一种的药理化合物变得困难，从而减少副作用。我们将研究eticlopride，一种D2R/D3R拮抗剂，如何与人源D3R结合。D3R和eticlopride配合物晶体结构PDB为3PBL。在此练习中，我们首先对蛋白和配体进行预对接准备，然后使用RosettaScripts进行对接，最后进行结果分析。

该教程也是一个学习RosettaScripts的材料。同时需要注意的是，rosetta版本需要在2017年之后(xml格式经过修改)

预处理

准备D3R蛋白结构
- 创建protein_prep、ligand_prep、docking文件夹
- 使用clean_pdb.py自动下载PDB文件并去除掉蛋白坐标信息之外的信息，A表示只获取A链
- clean_pdb.py输出两个文件：3PBL_A.pdb包含蛋白单链结构。3PBL_A.fasta包含对应的序列信息
- 3PBL_A.pdb为受体结构，被用于对接，将它复制到docking文件夹
  1
  2
  3
  4
  5
  6
  mkdir protein_prep
  mkdir ligand_prep
  mkdir docking
  cd protein_prep
  /path/to/rosetta/main/tools/protein_tools/scripts/clean_pdb.py 3PBL A
  cp 3PBL_A.pdb ../docking
准备配体文件，生成配体构象参数库文件

进入ligand_prep文件夹

配体文件，有两种获取方式

eticlopride.sdf：3PBL蛋白内部包含配体文件，你可以从初始蛋白中将该配体提取出来。注意：该结构配体不含氢

eticlopride_conformers.sdf: 此文件为eticlopride的构象库。下载的PDB中只包含配体的唯一构象，对于对接来说是远远不够的。因此必须扩展构象库以期达到对构象空间的充分采样效果。与此同时，我们也需要对配体加氢。即先加氢后生成构象库的生成方法

Meiler lab’s BioChemicalLibrary （目前无法访问）
OpenEye’s MOE software
Frog 2.1
DG-AMMOS

openbabel（推荐，可以本地使用）

# 使用遗传算法生成30个构象
obabel Conformer.sdf -O ga_conformers.sdf --conformer --nconf 30  --writeconformers
# 使用Confab生成分子的低能量构象分子集合
obabel Conformer.sdf -O ga_conformers.sdf --confab --conf 100

生成一个.params文件，并将PDB构象与eticlopride的Rosetta原子类型相关联。params文件是配体对接所必需的，因为Rosetta没有自定义记录数据库中的小分子. 生成的ETQ.params包含Rosetta处理配体ETQ.pdb包含生成的第一构象和剩余构象的必要信息。PDB包含其余的构象库。

cd ../ligand_prep
# 对获取的配体加完氢之后生成构象库(加氢操作可以使用DS、pymol、gaussian等软件实现)
obabel eticlopride.sdf -O eticlopride_conformers.sdf --confab --conf 50
/path/to/rosetta/main/source/scripts/python/public/molfile_to_params.py -n ETQ -p ETQ --conformers-in-one-file eticlopride_conformers.sdf
cp ETQ* ../docking

将生成的配体文件复制到docking文件夹中，检查文件夹中的文件是否齐全。
- 3PBL_A.pdb：蛋白受体结构
  - ETQ.pdb：配体的默认初始结构
  - ETQ_conformers.pdb：包含所有配体构象库的pdb文件
  - ETQ.params：Rosetta参数文件，提供配体的必要参数
检查文件夹中的XML文件(dock.xml)，输入选项文件(options.txt)，复合物结构文件(crystal_complex.pdb)
- dock.xml：RosettaScripts XML脚本文件告诉Rosetta如何采样和打分。定义打分函数、提供低解析度粗粒化采样和高精度蒙特卡洛采样参数
- options.txt：告诉Rosetta输入参数
- crystal_complex.pdb：D3-eticlopride复合物pdb。它可以作为正确的答案，让我们比较模型和实际结构

运行docking

# options.txt
#Pound signs indicate comments 
#-in:file:s option imports the protein and ligand PDB structures
#-in:file:extra_res_fa option imports the parameters for the ligand
-in
	-file
		-s '3PBL_A.pdb ETQ.pdb'
		-extra_res_fa ETQ.params
#the packing options allow Rosetta to sample additional rotamers for
#protein sidechain angles chi 1 (ex1) and chi 2 (ex2) 
#no_optH false tells Rosetta to optimize hydrogen placements
#flip_HNQ tells Rosetta to consider HIS,ASN,GLN hydrogen flips
#ignore_ligand_chi prevents Roseta from adding additional ligand rotamer
-packing
	-ex1
	-ex2
	-no_optH false
	-flip_HNQ true
	-ignore_ligand_chi true
#parser:protocol locates the XML file for RosettaScripts
-parser
	-protocol dock.xml
#overwrite allows Rosetta to write over previous structures and scores
-overwrite
#Ligand docking is not yet benchmarked with the updated scoring function
#This flag restores certain parameters to previously published values
-mistakes
	-restore_pre_talaris_2013_behavior true 
# 指定对接结构生成路径
-out 
	-path
		-all ./out

dock.xml
<ROSETTASCRIPTS>
		<SCOREFXNS>
			<ScoreFunction name="ligand_soft_rep" weights="ligand_soft_rep">
			</ScoreFunction>
			<ScoreFunction name="hard_rep" weights="ligand">
			</ScoreFunction>
		</SCOREFXNS>
		<LIGAND_AREAS>
			<LigandArea name="inhibitor_dock_sc" chain="X" cutoff="6.0" add_nbr_radius="true" all_atom_mode="false"/>
			<LigandArea name="inhibitor_final_sc" chain="X" cutoff="6.0" add_nbr_radius="true" all_atom_mode="false"/>
			<LigandArea name="inhibitor_final_bb" chain="X" cutoff="7.0" add_nbr_radius="false" all_atom_mode="true" Calpha_restraints="0.3"/>
		</LIGAND_AREAS>
		<INTERFACE_BUILDERS>
			<InterfaceBuilder name="side_chain_for_docking" ligand_areas="inhibitor_dock_sc"/>
			<InterfaceBuilder name="side_chain_for_final" ligand_areas="inhibitor_final_sc"/>
			<InterfaceBuilder name="backbone" ligand_areas="inhibitor_final_bb" extension_window="3"/>
		</INTERFACE_BUILDERS>
		<MOVEMAP_BUILDERS>
			<MoveMapBuilder name="docking" sc_interface="side_chain_for_docking" minimize_water="false"/>
			<MoveMapBuilder name="final" sc_interface="side_chain_for_final" bb_interface="backbone" minimize_water="false"/>
		</MOVEMAP_BUILDERS>
		<SCORINGGRIDS ligand_chain="X" width="15">
			<ClassicGrid grid_name="class" weight="1.0"/>
		</SCORINGGRIDS>
		<MOVERS>
			<Transform name="transform" chain="X" box_size="7.0" move_distance="0.2" angle="20" cycles="500" repeats="1" temperature="5"/>
			<HighResDocker name="high_res_docker" cycles="6" repack_every_Nth="3" scorefxn="ligand_soft_rep" movemap_builder="docking"/>
			<FinalMinimizer name="final" scorefxn="hard_rep" movemap_builder="final"/>
			<InterfaceScoreCalculator name="add_scores" chains="X" scorefxn="hard_rep" native="crystal_complex.pdb"/> 
		</MOVERS>
		<PROTOCOLS>
			<Add mover="transform"/>
			<Add mover="high_res_docker"/>
			<Add mover="final"/>
			<Add mover="add_scores"/>
		</PROTOCOLS>
</ROSETTASCRIPTS>

1
2
3

cd docking 
mkdir out
/path/to/rosetta/main/source/bin/rosetta_scripts.linuxgccrelease @options.txt -nstruct 500 -database /path/to/rosetta/main/database/

Rosetta生成的模型以3PBL_A_ETQ_前缀+四位数标识符保存。每个PDB中均包含坐标和模型相关的Rosetta分数。除此之外，模型分数也总结在score.sc中，两个主要的打分项为：
- total_score：total score为整个蛋白-配体复合物的打分，可以用来进行模型评估
- interface_delta_X：该分数为蛋白-配体复合物与蛋白-配体单体之间的差值。Interface Score对于分析配体影响和比较不同的复合物来说十分有用。
另一个需要关注的指标为Transform_accept_ratio。这是在低精度变换网格搜索期间接受的蒙特卡洛采样比率。如果这个数是0或者非常低时，搜索空间可能限制太大，无法进行适当的采样。
在拥有正确晶体结构的基准测试示例中，ligand_rms_no_super_X将在crystal_complex.pdb给出模型配体与晶体结构配体之间的RMSD差值。这是对模型相关性进行基准测试时的一个重要指标现实。当晶体结构未知时，我们也可以用最佳的方法计算模型rmsd评分结构为“true answer”。
可以使用pymol查看对接的模型结构
scripts目录中的visualize_ligand.py脚本是进行快速蛋白质-配体界面可视化的有用方式。它接受一个PDB并通过settings生成一个.psePymol会话(前提：需要安装pymol)
1
2
/path/to/rosetta/main/tools/protein_tools/scripts/visualize_ligand.py 3PBL_A_ETQ_0001.pdb
pymol 3PBL_A_ETQ_0001.pse

分析

在随机生成了一系列对接模型之后，并不是每个模型都是有用的，因此我们需要对其进行分析

进入out文件夹
1
cd out
生成score_vs_rmsd.csv、rmsds_to_best_model.data、.png文件
- score.sc：Rosetta输出的所有模型的score文件
- score_vs_rmsd.csv：以逗号分隔的文件，第一列为文件名，第二列为复合物的总分数，第三列为界面分数、第四列为配体的RMSD值(与初始结构相比)。
  该文件使用extract_score.bash脚本生成。
  1
  /path/to/rosetta/main/demos/tutorials/ligand_docking/scripts/extract_scores.bash score.sc > score_vs_rmsd.csv
rmsds_to_best_model.data：以空格分隔的文件，所有对接模型和最佳打分模型之间的RMSD比较
文件第一列为文件名、第二列为所有重原子的RMSD、第三列为没有align的配体RMSD、第四列为align之后的配体RMSD、第五列为配体周围蛋白侧链重原子的RMSD。该文件由calculate_ligand_rmsd.py计算得出。使用pymol比较具有相同残基和配体原子pdb结构，这是计算配体rmsd的简便方法。
1
/path/to/rosetta/main/demos/tutorials/ligand_docking/scripts/calculate_ligand_rmsd.py -n 3PBL_A_ETQ_0003.pdb -c X -a 7 -o rmsds_to_best_model.data *_000*.pdb
该命令将所有的模型与-n指定的模型进行比较。
-c表示配体所在链为X链。
-a告诉脚本使用7埃作为侧链RMSD的球形截断。
-o 输出文件名。
最后，提供一系列pdb。
PNG：使用上述生成数据进行绘图，使用python和matplotlib

对CSV文件进行排序，并输出前20个

1
2
3

sort -t, -nk3 score_vs_rmsd.csv | head -n 20
# 输出后20个，即不合理的结构
sort -t, -nk3 score_vs_rmsd.csv | tail -n 20

molfile_to_params.py 使用

molfile_to_params将小分子转换为Rosetta格式的.params残基文件，将配体构象原子输出为PDB HETATMs。

1	Usage: molfile_to_params.py [flags] { INPUT.mol \| INPUT.sdf \| INPUT.mol2 }

-h,--help ：帮助文件
-n NM, --name=NM：残基名 NM1，NM2...
-p FILE, --pdb=File：pdb文件前缀
-c, --centroid：以Rosetta质心模型输出文件
--chain=CHAIN：配体链名
--center=X, Y, Z：将pdb坐标平移至给定的重原子质心
-m MAX, --max-confs=MAX：如果超过总数，不扩展proton chis
--root_atom=ATOM_NUM：指定root原子（index从1开始）
--nbr_atom=ATOM_NUM：指定nbr原子
-k FILE, --kinemage=FILE：输出配体拓扑
-a ALA, --amino-acid=ALA：建立修改后的氨基酸残基参数文件，只允许mol2文件
--clobber：覆盖存在文件
--no-param：不输出参数文件
--no-pdb：不输出pdb文件
--extra_torsion_output：输出额外的旋转角文件
--keep-names：除重复之外，原子名保留
--long-names：如果指定名超过3个字符，不进行截断
--recharge=CHG：忽略点电荷，设置总电荷
--m-ctrl=FILE：从FILE中读取额外的M control
--mm-as-virt：assign mm atom types as VIRT, rather than X
--skip-bad-conformers：如果conformer中原子顺序错误，直接跳过
--conformers-in-one-file：Output 1st conformer to NAME.pdb and all others to NAME_conformers.pdb

xml解析

<SCOREFXNS>用于指定打分函数，可参考之前推文

<LIGAND_AREAS>

1
2
3

<LIGAND_AREAS>
<LigandArea name="[name_of_this_ligand_area]" chain="[string]" cutoff="[float]" add_nbr_radius="[true|false]" all_atom_mode="[true|false]" minimize_ligand="[float]" Calpha_restraints="[float]" high_res_angstroms="[float]" high_res_degrees="[float]" tether_ligand="[float]" />
</LIGAND_AREAS>

LIGAND_AREAS描述了特定于每个配体的参数，用于多个配体对接的研究。cutoff为氨基酸的beta-C原子与配体之间的阶段距离，该残基也被视作interface的一部分。all_atom_mode可以设为true或false。为真时，只要蛋白残基上的beta-C原子与配体任意原子的距离在截断距离之内时，残基即可视为界面的一部分。为假时，蛋白残基上的beta-C原子必须与配体相邻原子的距离在截断距离之内。如果add_nbr_radius为真，截断值会随着蛋白残基邻域半径的增大而增大。邻域半径为残基repack时移动范围的估计值，添加它可以用于补偿repack时蛋白侧链的移动。通过指定大于0的minimize_ligand来开启配体最小化，该值表示配体扭转角旋转的标准差大小(degree)。设置Calpha_restraints大于0，使骨架flexiable，该值表示Calpha移动的标准差大小(埃)。

在高精度对接过程中，少量的配体平移和旋转与rotamer trials或者repack周期相结合。这些值可以分别由'high_res_angstrom'和'high_res_degrees'值控制。tether_ligand值(以埃为单位)将限制配体，以便多个小平移循环叠加起来不会产生大平移。

<INTERFACE_BUILDERS>
1
2
3
<INTERFACE_BUILDERS>
<InterfaceBuilder name="[name_of_this_interface_builder]" ligand_areas="(comma separated list of predefined ligand_areas)" extension_window="(int)"/>
</INTERFACE_BUILDERS>
Interface builder描述了如何选择蛋白-配体界面的残基，这些残基被用于repacking、rotamers trials和骨架最小化。ligand_areas是一个用逗号分隔的字符串列表，匹配之前描述的ligand_areas。extension_window用于标记为'near interface'的残基包围界面残基。这对主链最小化很重要，因为残基的主链不能真正移动，除非它是柔性残基延伸的一部分
<MOVEMAP_BUILDERS>
1
2
3
<MOVEMAP_BUILDERS>
<MoveMapBuilder name="[name_of_this_movemap_builder]" sc_interface="(string)" bb_interface="(string)" minimize_water="[true|false]"/>
</MOVEMAP_BUILDERS>
movemap builder构造一个movemap，movemap是一个2xN的布尔值列表，N表示蛋白/配体复合物的残基数目。两列分别用于骨架和侧链移动。MoveMapBuilder将骨架和侧链界面结合起来，如果不想最小化骨架结构，就不要指定bb_interface。minimize_water选项为全局选项，如果对接过程中将水视为额外的残基(多配体对接)，则应通过LIGAND_AREAS和INTERFACE_BUILDERS指定。

1
2
3

<SCORINGGRIDS ligand_chain="(string)" width="(real)" name="(string, optional)" resolution="(real, optional)" normalize_mode="(string, optional)" >
<[scoring grid type name] grid_name="(string)" weight="(real)"/>
</SCORINGGRIDS>

SCORINGGRID被用于定义ligand打分网格(目前只被用于Tranform mover)。网格起初在指定配体链中心创建，且是以埃为单位的指定宽度立方体。Rosetta自动决定网格是否需要重新计算。如果任何非配体原子改变位置，网格将被重新计算。为每个网格指定的权重乘以该网格的配体分数。当前可以指定以下grid_types:

ClassicGrid:最初在Davis, 2009年实现的RosettaLigand中使用的评分网格。
AtrGrid:一个有引力的势能网格。
ChargeGrid: 电场的近似。
HbdGrid:一种基于知识的势能衍生网格，近似于蛋白质氢键供体相互作用。
HbaGrid:一种基于知识的势能衍生网格，近似于蛋白质氢键受体相互作用。
LigandPropertyScore: ?(目前无文献)。
RepGrid:排斥势网格。
ShapeGrid:一种基于知识的势能导出网格，近似形状互补。
SiteGrid: ?(目前无文献)。
SolvationGrid:溶剂化势网格。

顶层的“name”参数用于指定GridSet名称，该名称可以与Transform mover一起使用，以选择正在使用的网格集。(如果没有给出名称，网格集将以“default”名称加载。)

<MOVERS> <StartFrom>

<StartFrom name="(&string)" chain="(&string)" use_nbr="(0 &bool)" >
   <Coordinates x="(&float)" y="(&float)" z="(&float)" pdb_tag="('' &string)"/>
   <File filename="(&string)" />
   <PDB filename="(&string)" atom_name="('' &string)" pdb_tag="('' &string)" />
</StartFrom>

StartFrommover将链(通常是配体)移动到指定的坐标。这用于将对接配体放置在与输入PDB不同的位置。如果指定了多个坐标，则开始坐标将随机选择。默认情况下，指定链的质心(所有原子/残基的平均位置)以给定坐标为中心。如果设置了"use_nbr"标志，则将相邻原子(在异构体重新包装过程中叠加的原子)放置在给定位置。设置use_nbr意味着给定的链由单个残基组成。

坐标可以由一个或多个子标签指定。可以同时指定多个不同类型的子标签——可能的坐标集是所有子标记的坐标的组合。

每个子标签都可以接受一个可选的结构说明符(pdb_tag、file_name或hash)。如果结构与给定的结构说明符匹配，则只考虑该结构说明符的坐标。如果一个结构不匹配任何提供的结构说明符，那么将使用默认坐标(那些没有指定结构说明符的坐标)。

Coordinates subtag

每个标记包含一个XYZ坐标位置，该位置将被添加到有效坐标集中。

File subtag

提供一个包含起始位置的JSON格式文件。如果要将配体连接到大量的蛋白质结构中，这是非常有用的。文件格式为:

[
    {
        "file_name" : "infile.pdb",
        "x" : 0.0020,
        "y" : -0.004,
        "z" : 0.0020,
        "hash" : 14518543732039167129
    }
]

file_name和hash均为可选参数。如果两者都存在，则忽略“file_name”条目。

“哈希”是无配体结构的哈希表示。目前，浮点数的Boost散列非常依赖于构建和平台。您不应该认为这些哈希值是可移植的。这个问题将在将来解决。

PDB subtag

PDB文件中的重原子位置为起始坐标。默认情况下使用所有重原子位置。如果给出atom_name，则只使用与该原子名匹配的原子位置。注意，标准的Rosetta PDB加载机制被绕过了，所以对于所提供的文件中的所有残基/原子，您不需要params文件等。(尽管所有剩余数/原子名组合应该是唯一的。)

生成该文件的一种简单方法是图形化。如果你将初始结构和包含水的PDB一起加载到PyMol中(可以是相同的结构)，你可以使用“3键编辑”或“2键编辑”鼠标模式。按住Ctrl键(Mac上的命令)，你可以左键点击/拖动移动单个水氧到你感兴趣的起始位置。切换回查看模式，选择您移动的单个原子，然后File->Save Molecule并将该选择(只是单个水分子)保存为一个新的PDB文件。其他结构查看程序可能有类似的功能。详情请参阅手册。

注意:PDB应该只包含描述起始点的原子。它不应该包含任何上下文残基(例如，任何受体蛋白残基，或没有定义结合袋假定中心的配体原子)。

小结

此为一个英文教程的复现，具体操作时需要自己去调参数
尤为重要的是你的配体和蛋白的初始坐标对于对接至关重要。在原教程中直接使用的解析结构，且配体分子分子量较小，因此不容易出错，但是实际案例时，需要考虑多方面的因素
比如说，对于配体分子对接来说，首先就要设定其初始坐标，即在MOVER模块中加入StartFrom，在原教程中并没有添加。
Rosetta的功能太过强大，参数也十分的多。本文也是根据教程改写而来，笔者能力精力有限，读者可以自行开发。