CASEA的奇思妙想

RosettaLigand

发表于 分类于
本文字数： 13k 阅读时长 ≈ 12 分钟

RosettaLigand

RosettaLigand

前言

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

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

预处理

  1. 准备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

  2. 准备配体文件，生成配体构象参数库文件

  • 进入ligand_prep文件夹

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

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

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

      • Meiler lab’s BioChemicalLibrary （目前无法访问）

      • OpenEye’s MOE software

      • Frog 2.1

      • DG-AMMOS

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

        1
        2
        3
        4
        5
        # 使用遗传算法生成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包含其余的构象库。

    1
    2
    3
    4
    5
    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

  1. 将生成的配体文件复制到docking文件夹中，检查文件夹中的文件是否齐全。

    • 3PBL_A.pdb：蛋白受体结构
      • ETQ.pdb：配体的默认初始结构
      • ETQ_conformers.pdb：包含所有配体构象库的pdb文件
      • ETQ.params：Rosetta参数文件，提供配体的必要参数

  2. 检查文件夹中的XML文件(dock.xml)，输入选项文件(options.txt)，复合物结构文件(crystal_complex.pdb)

    • dock.xml：RosettaScripts XML脚本文件告诉Rosetta如何采样和打分。定义打分函数、提供低解析度粗粒化采样和高精度蒙特卡洛采样参数
    • options.txt：告诉Rosetta输入参数
    • crystal_complex.pdb：D3-eticlopride复合物pdb。它可以作为正确的答案，让我们比较模型和实际结构

  3. 运行docking

    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
    32
    # 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
    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
    32
    33
    34
    35
    36
    37
    38
    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/

  4. Rosetta生成的模型以3PBL_A_ETQ_前缀+四位数标识符保存。每个PDB中均包含坐标和模型相关的Rosetta分数。除此之外，模型分数也总结在score.sc中，两个主要的打分项为：

    • total_score：total score为整个蛋白-配体复合物的打分，可以用来进行模型评估
    • interface_delta_X：该分数为蛋白-配体复合物与蛋白-配体单体之间的差值。Interface Score对于分析配体影响和比较不同的复合物来说十分有用。

  5. 另一个需要关注的指标为Transform_accept_ratio。这是在低精度变换网格搜索期间接受的蒙特卡洛采样比率。如果这个数是0或者非常低时，搜索空间可能限制太大，无法进行适当的采样。

  6. 在拥有正确晶体结构的基准测试示例中，ligand_rms_no_super_X将在crystal_complex.pdb给出模型配体与晶体结构配体之间的RMSD差值。这是对模型相关性进行基准测试时的一个重要指标现实。当晶体结构未知时，我们也可以用最佳的方法计算模型rmsd评分结构为“true answer”。

  7. 可以使用pymol查看对接的模型结构

  8. scripts目录中的visualize_ligand.py脚本是进行快速蛋白质-配体界面可视化的有用方式。它接受一个PDB并通过settings生成一个.pse Pymol会话(前提：需要安装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

分析

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

  1. 进入out文件夹

    1
    cd out

  2. 生成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

  3. 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。

  4. PNG：使用上述生成数据进行绘图，使用python和matplotlib

  5. 对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解析

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

  2. <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值(以埃为单位)将限制配体，以便多个小平移循环叠加起来不会产生大平移。

  3. <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’的残基包围界面残基。这对主链最小化很重要，因为残基的主链不能真正移动，除非它是柔性残基延伸的一部分

  4. <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指定。

  5. <SCORINGGRIDS ligand_chain=“X” width=“15”>

    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”名称加载。)

  6. <MOVERS>
    <StartFrom>

1
2
3
4
5
<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格式文件。如果要将配体连接到大量的蛋白质结构中，这是非常有用的。文件格式为:

1
2
3
4
5
6
7
8
9
[
    {
        "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的功能太过强大，参数也十分的多。本文也是根据教程改写而来，笔者能力精力有限，读者可以自行开发。