P450酶rosetta对接，设计时的HEM ligated CYM问题

　　P450酶是一类以HEM作为辅酶，并是在催化过程中严重依赖HEM氧化还原循环的一类酶。

　　在天然环境下，HEM辅酶的Fe原子会与蛋白中一个处于垂直位置的脱氢CYS（CYM）的S形成一个共价键，此共价键基本垂直于HEM平面，其间的距离、角度、二面角都保持在一直相对稳定的状态。

　　在实际的计算模拟过程中，尤其是用rosetta做对接或者设计的过程中，我们不可避免的会遇到CYM与HEM进行共价键连接的问题，因为rosetta不能对CYM以及HEM这类非天然残基进行识别，且在计算中还需要考虑CYM与HEM间的共价键，所以这是一个异常棘手的问题，今天我就来分享一种解决方案，此方案是我一个师弟偶然中想出来的。

　　在我以前的方案中，我会设计一个全新的非天然残基CYM，基本类似于CYS，区别在于少了一个HG，此CYM需要使用molfiletoparams命令生成相应的params，在计算中要考虑CYM与上下游氨基酸的肽键连接以及与HEM-Fe的共价连接，所以在matchCST里面至少要有3个block，这大大的增加了计算的复杂度及操作难度，因为任何一个小小的错误都有可能造成CST的失误并最终使对接和设计失败。

　　今天要提出的方案中乃是源于使用老方案对一个P450酶进行对接时，报的一个莫名的错误，且找不到出错原因，所以师弟在苦思冥想中提出了一个全新的方案，这个方案在设计matchCST时，仅需要一个block，不用额外考虑CYM与上下游氨基酸的肽键连接问题，下面对这个方案进行详解。

　　一、准备一个CYM的params文件，将其添加到rosetta fa_standard的l-caa识别库中（其他talaris库也可以以此方式添加） （~/rosetta_src/main/database/chemical/residue_type_sets/fa_standard/residue_types/l-caa） ，且这个params应该继承rosetta库里天然CYS的属性，这样我们就无需额外构建CYS与上下游氨基酸的肽键，下面的CYM.params脱胎于rosetta标准库里的CYS.params文件，我们删除了HG原子及相关信息，修改了CYS的名称，AA后改为UNK，删除了ROTAMER_AA，删除了PROTON_CHI项，保留了其他的关键条目，保证rosetta可以对其进行识别，并构建肽键：

CYM.params

#rosetta residue topology file
NAME CYM
IO_STRING CYM Z
TYPE POLYMER #residue type
AA UNK
ATOM  N   Nbb  NH1  -0.6046255 -0.350
ATOM  CA  CAbb CT1   0.0900506  0.100
ATOM  C   CObb C     0.6884871  0.550
ATOM  O   OCbb O    -0.6884871 -0.550
ATOM  CB  CH2  CT2  -0.1178426  0.000
ATOM  SG  SH1  S    -0.2463981 -0.290
ATOM  H   HNbb H     0.3987955  0.250
ATOM  HA  Hapo HB    0.1157793  0.000
ATOM 1HB  Hapo HA    0.0964167  0.000
ATOM 2HB  Hapo HA    0.0964167  0.000

LOWER_CONNECT N
UPPER_CONNECT C
BOND  N    CA 
BOND  N    H  
BOND  CA   C  
BOND  CA   CB 
BOND  CA   HA 
BOND_TYPE  C    O 2
BOND  CB   SG 
BOND  CB  1HB 
BOND  CB  2HB 
CHI 1  N    CA   CB   SG 
PROPERTIES PROTEIN ALPHA_AA L_AA SC_ORBITALS METALBINDING
METAL_BINDING_ATOMS O SG
DISULFIDE_ATOM_NAME SG
NBR_ATOM CB
# APL O to CB distance -- max O-CB observed in ubiquitin
NBR_RADIUS 3.4473
FIRST_SIDECHAIN_ATOM CB
RAMA_PREPRO_FILENAME all.ramaProb prepro.ramaProb
ACT_COORD_ATOMS SG END
ICOOR_INTERNAL    N      0.000000    0.000000    0.000000   N     CA    C  
ICOOR_INTERNAL    CA     0.000000  180.000000    1.458001   N     CA    C  
ICOOR_INTERNAL    C      0.000000   68.800003    1.523258   CA    N     C  
ICOOR_INTERNAL  UPPER  150.000000   63.800018    1.328685   C     CA    N  
ICOOR_INTERNAL    O    180.000000   59.200008    1.231015   C     CA  UPPER
ICOOR_INTERNAL    CB  -121.600000   69.400000    1.528861   CA    N     C  
ICOOR_INTERNAL    SG     0.000000   65.900000    1.808803   CB    CA    N  
ICOOR_INTERNAL   1HB   121.200000   70.500000    1.090249   CB    CA    SG 
ICOOR_INTERNAL   2HB   117.600000   70.500000    1.089821   CB    CA   1HB 
ICOOR_INTERNAL    HA  -119.000000   71.500000    1.090059   CA    N     CB 
ICOOR_INTERNAL  LOWER -150.000000   58.299995    1.328685   N     CA    C  
ICOOR_INTERNAL    H    180.000000   60.849998    1.010000   N     CA  LOWER

　　二、更改fa_standard文件夹下的residue_types.txt，在记录CYS条目的行下添加CYM的记录行：

residue_types/l-caa/CYM.params

       三、修改fa_standard文件夹下的exclude_pdb_component_list.txt，在##The 20 canonical amino acids下加入CYM；

　　四、修改prepro.ramaProb和all.ramaProb文件(/rosetta/main/database/scoring/score_functions/rama/fd/目录内)，在这里直接复制的该文件里关于CYS的记录条目，改成相应的CYM，文件太大，这里就不上传了。这样就实现了自定义拓展rosetta的params库，rosetta就可以识别我们模拟分子里的CYM残基了。

　　五、准备CYM与HEM的matchCST文件，准备方法可以参考rosetta geometric constraint file（用于match和design） - 丨o聽乄雨o丨 - 博客园 (cnblogs.com)和enzyme design 整体流程及感想 - 丨o聽乄雨o丨 - 博客园 (cnblogs.com)这两篇博客，下面也附上我们所用的CST文件内容及options和XML文件。

LIG.CST

#block 1 for covalent bond for CYM and CPD1

CST::BEGIN  
  TEMPLATE::   ATOM_MAP: 1 atom_name: FE1 N2 N1
  TEMPLATE::   ATOM_MAP: 1 residue3: HEM

  TEMPLATE::   ATOM_MAP: 2 atom_name: SG CB CA
  TEMPLATE::   ATOM_MAP: 2 residue3: CYM 
 
  CONSTRAINT:: distanceAB:    2.00   0.20 180.00  1
  CONSTRAINT::    angle_A:   98.00   3.00 100.00  360.00
  CONSTRAINT::    angle_B:  120.00  10.00  50.00  360.00
  CONSTRAINT::  torsion_A:  -98.00   5.00  50.00  360.00
  CONSTRAINT::  torsion_B:  103.00   5.00  25.00  360.00
  CONSTRAINT:: torsion_AB: -176.60   5.00   5.00  360.00
CST::END

ligand_dock_p450.flags

-in:file:s Sbai.pdb
-in:file:extra_res_fa HEM.params AGI.params
-run::preserve_header
-packing
    -ex1
    -ex2aro
    -ex2 
    -no_optH false
    -flip_HNQ true
    -ignore_ligand_chi true
-enzdes
    -cstfile LIG.cst
-parser
    -protocol ligand_dock_2018.xml

-out
#    -level 100
    -nstruct 100
    #-overwrite

ligand_dock_2018.xml

<ROSETTASCRIPTS>This protocol will simply do low-resolution followed by high-resolution docking.
It will also report the binding energy (ddg) and buried-surface area (sasa) in the score file.
    <SCOREFXNS>
        <ScoreFunction name="ligand_soft_rep" weights="ligand_soft_rep">
        </ScoreFunction>
        <ScoreFunction name="hard_rep" weights="ligand">
        </ScoreFunction>
    </SCOREFXNS>

    <LIGAND_AREAS>
        <LigandArea name="docking_sidechain_X" chain="X" cutoff="6.0" add_nbr_radius="true" all_atom_mode="true" minimize_ligand="10"/>
        <LigandArea name="final_sidechain_X" chain="X" cutoff="6.0" add_nbr_radius="true" all_atom_mode="true"/>
        <LigandArea name="final_backbone_X" chain="X" cutoff="7.0" add_nbr_radius="false" all_atom_mode="true" Calpha_restraints="0.3"/>
        
        <LigandArea name="docking_sidechain_F" chain="F" cutoff="6.0" add_nbr_radius="true" all_atom_mode="true" minimize_ligand="10"/>
        <LigandArea name="final_sidechain_F" chain="F" cutoff="6.0" add_nbr_radius="true" all_atom_mode="true"/>
        <LigandArea name="final_backbone_F" chain="F" 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="docking_sidechain_X,docking_sidechain_F"/>
        <InterfaceBuilder name="side_chain_for_final" ligand_areas="final_sidechain_X,final_sidechain_F"/>
        <InterfaceBuilder name="backbone" ligand_areas="final_backbone_X,final_backbone_F" extension_window="3"/>
    </INTERFACE_BUILDERS>
    <MOVEMAP_BUILDERS>
        <MoveMapBuilder name="docking" sc_interface="side_chain_for_docking" minimize_water="true"/>
        <MoveMapBuilder name="final" sc_interface="side_chain_for_final" bb_interface="backbone" minimize_water="true"/>
    </MOVEMAP_BUILDERS>
    <SCORINGGRIDS ligand_chain="X" width="20">
        <ClassicGrid grid_name="classic" weight="1.0"/>
    </SCORINGGRIDS>
    <MOVERS>
    single movers_X
        <AddOrRemoveMatchCsts name="cstadd" cst_instruction="add_new"/> add catalytic constraints
        <Transform name="transform_F" chain="F" box_size="7.0" move_distance="0.2" angle="20" cycles="700" repeats="1" temperature="5"/>
        <Transform name="transform_X" chain="X" box_size="8.0" move_distance="0.2" angle="20" cycles="700" repeats="1" temperature="5"/>
        <AddOrRemoveMatchCsts name="cstrem" cst_instruction="remove" keep_covalent="1"/> remove constraints
        <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"/>
        <AddOrRemoveMatchCsts name="cstfinadd" cst_instruction="add_pregenerated"/>
        <InterfaceScoreCalculator name="add_scores" chains="X,F" scorefxn="hard_rep"/>
        
    compound movers
        <ParsedProtocol name="low_res_dock">
            <Add mover_name="cstadd"/>
            <Add mover_name="transform_F"/>
            <Add mover_name="transform_X"/>
            <Add mover_name="cstrem"/>
        </ParsedProtocol>
        <ParsedProtocol name="high_res_dock">
            <Add mover_name="high_res_docker"/>
            <Add mover_name="final"/>
            <Add mover_name="cstfinadd"/>
        </ParsedProtocol>
    </MOVERS>
    <PROTOCOLS>
        <Add mover_name="low_res_dock"/>
        <Add mover_name="high_res_dock"/>
        <Add mover_name="add_scores"/>        
    </PROTOCOLS>
</ROSETTASCRIPTS>