环肽

线状多肽作为药物可能天生就不稳定，在细胞内蛋白水解的可能性很高。考虑到游离酸和游离胺位于两端，它们通常也是极性的。环化促进了环结构内的分子内氢键，降低了分子的外部氢键能力，与无环前体相比，这降低了分子的极性，增加了化合物的膜渗透性。环化引起的两个β-turn局部二级结构也会减少极性表面积，这通常会增强细胞的通透性。 可根据构成环的键的类型对环肽进行分类，通常可将其分为均环肽和杂环肽，均环肽是其中环仅由正常肽键组成（即一个残基的α羧基与另一个残基的α胺之间）的环肽，如环孢菌素A。杂环肽至少含有一个非α-酰胺键，如微囊藻毒素和杆菌肽中，一个残基的侧链与另一个残基的α羧基之间的键。 除此以外还可以根据肽的环肽位置将可以环肽分为四种类型：头尾，头对侧链，侧链对尾和侧链对侧链环化。

环肽

环肽的特点

环肽的应用特点：

• 具有明确的固定构象，能够与受体很好地契合
• 环肽的代谢稳定性和生物利用度远远高于直链肽
• 由于环状结构导致的构象变化限制，环肽类化合物一般具有较大的表面积，从而使其与靶点蛋白间具有高度亲和力及识别特异性
• 相比直链多肽用量极少就能达到高生物活性。功能等同于多种线性多肽的活性总和

Rosetta中环肽设计的两个重要的算法

Kinematic Closure (KIC) 运动学闭合

运动学闭合Kinematic Closure (KIC)算法最初是为机器人领域开发的，用于解决确定将机器人的手或脚放置在所需位置所需的关节角度的问题。Rosetta已经对它们进行了调整，以便在Rosetta software suite中使用，以对具有明确起点和终点的原子链的构象进行采样。 分子运动学闭合问题可以描述如下：给定分子内原子的共价连续链，链末端的共价键固定链的起点和终点，什么可能的构象保持链内键长、键角和二面角限制的完整性？为了解决这样的问题，我们将原子链分为两段，并在链的起点（ the first pivot ）、链的终点（ the last pivot ）和两段之间的断点（ the middle pivot ）定义"pivot points" 。 做到这一点后，两个部分内的自由度可以保持固定、随机、扰动或以其他方式进行改变。这些自由度包括键长、键角和二面角。无论对这些自由度做什么，最终都会有两个段，它们仍然具有明确的刚体变换，从第一个枢轴到中间枢轴（在第一段中），以及从中间枢轴到最后一个枢轴（在第二段中）。然后随机地扰动所有的非节点氨基酸的二面角、键角、键长，从而生成新的骨架构象。最后通过计算调整三个节点氨基酸的骨架二面角来使先前的“切点”闭合。 GeneralizedKIC mover只能对现有几何结构的构象进行采样。它既不能在pose上添加氨基酸残基，也不能在残基之间建立新的键。

simple_cycpep_predict

simple_cycpep_predict应用旨在快速采样主链环化限制性多肽的闭合构象。这些肽可以由L-或D-氨基酸残基、非手性氨基酸残基，类蛋白残基，或L-或D-低聚脲残基的任何混合物组成。同时可以指定构象中必须具备一定数量的主链氢键。

simple_cycpep_predict将多肽骨架当成是一段Loop结构进行处理, 其中包括若干个处理的阶段:

1. 生成线性的多肽结构，根据给定的初始序列正确设置环状多肽的氨基酸连接关系，随机设定N端和C端的位置。比如输入的序列为ALA LYS PHE ASP DILE PRO，可能会被随机处理为PHE ASP DILE PRO ALA LYS。
2. 根据Ramachandran分布随机设定每个氨基酸的phi和psi角的值，默认设定omega角设定为180°(顺式)。如上图所示Break peptide，随机设定二面角后环形的多肽出现了断点。
3. GenKIC将随机选定一个氨基酸作为锚定点(排除最后两个氨基酸)，根据每个氨基酸的偏好不同随机选取骨架psi和phi二面角进行干扰(导致构象开环)，后续进行KIC闭环处理。闭环后生成的多个候选构象将进一步通过Filter进行过滤(排除能量较差或几何结构异常的闭环构象)，并返回能量最佳者的构象。
4. 如果设定了需要形成分子内的二硫键，那么程序将尝试使用TryDisulfPermutations Mover对骨架中所有可能存在二硫键的区域进行设计。
5. 如果设定了在特定位置形成Cross-link, 程序将使用CrosslinkerMover来进行化学共价连接两个位点，目前兼容使用的是1,3,5-tris(bromomethyl)benzene(TBMB)和trimesic acid(TMA)。
6. 进行进一步的骨架和侧链构象优化，可选FastRelax、Cartesian FastRelax，又或FastDesign对氨基酸进行设计（此处的FastDesign聊胜于无）。
7. 输出完成预测或设计的约束肽结构。
8. (非算法的步骤)完成成百上千个结构后，判断是否有能量漏斗出现，选出若干个靠谱的构象进行分子动力学模拟，观测设计的构象是否稳定。

文章解读

Accurate de novo design of membrane-traversing macrocycles

这篇文章的主要内容是关于如何通过计算机设计和实验验证来系统地研究大环肽膜渗透性和口服生物利用度的设计原则。作者们设计了184个6-12个氨基酸残基的大环肽，这些大环肽具有广泛的预测结构，包含非规范背景修饰，并确定了35个的实验结构，其中29个与计算模型非常接近。通过这种控制，作者们表明，可以通过确保所有酰胺（NH）基团参与内部氢键相互作用来系统地实现膜渗透性。在6-12个氨基酸残基大小范围内，84个设计具有表观渗透性大于1 3 10^-6 cm/s的跨膜能力。暴露NH基团的设计可以通过设计在脂质膜中偏爱的另一种等能全氢键状态来使其具有膜渗透性。这种能够稳健地设计具有高结构准确性的膜渗透性和口服生物利用度的肽类应该有助于下一代设计的大环肽治疗药物。

设计策略

策略1：尽可能减少暴露的极性基团，增加疏水性

• generalized kinematic closure (genKIC) 采样6-12个长度的多聚甘氨酸主链，把有2对以上主链氢键形成的构象取出；
• 在有内部氢键相互作用的主链集合上，引入为L-pro,D-pro和AIB为刚性结构支撑位点，强化环肽的结构刚性； （
• 与之前环肽生成算法最大的不同点）消除暴露极性基团，在Rosetta设计序列阶段将主链上存在不饱和NH的氨基酸替换为主链甲基化的氨基酸类型，并且侧链序列设计空间只允许非极性残基；
• 筛选<5个主链甲基化残基以及内部存在2个或以上氢键配对的序列进行下一步验证conformational energy landscape，把折叠态即为势能面最低点的序列筛选出来。 根据ABOXY的rama二面角进行聚类。

**接下来将对他的代码进行分析：xml文件的格式可参考
**

1. SCOREFXNS:首先定义了四个打分函数，根据不同的应用场景设置打分函数中不同打分类别的权重：

    #Define scorefunctions to be used later in the script 
<SCOREFXNS>
     <ScoreFunction name="SFXN_STD" weights="beta_nov16.wts" >
        <Reweight scoretype="atom_pair_constraint" weight="1" />
        <Reweight scoretype="dihedral_constraint" weight="1" />
        <Reweight scoretype="angle_constraint" weight="1" />
        <Reweight scoretype="aa_composition" weight="1.0" />

     </ScoreFunction>    
     <ScoreFunction name="SFXN_STD_high_hbond" weights="beta_nov16.wts" >
         <Reweight scoretype="atom_pair_constraint" weight="1" />
         <Reweight scoretype="dihedral_constraint" weight="1" />
         <Reweight scoretype="angle_constraint" weight="1" />
         <Reweight scoretype="aa_composition" weight="1.0" />
         <Reweight scoretype="hbond_sr_bb" weight="5.0" />
         <Reweight scoretype="hbond_lr_bb" weight="5.0" />

    </ScoreFunction>
     <ScoreFunction name="SFXN_CART" weights="beta_nov16_cart.wts" >
         <Reweight scoretype="atom_pair_constraint" weight="1" />
         <Reweight scoretype="dihedral_constraint" weight="1" />
         <Reweight scoretype="angle_constraint" weight="1" />
         <Reweight scoretype="aa_composition" weight="1.0" />

     </ScoreFunction>     
    <ScoreFunction name="SFXN_hbond_bb" weights="empty.wts" symmetric="0">
         <Reweight scoretype="fa_rep" weight="0.1" />
         <Reweight scoretype="fa_atr" weight="0.2" />
         <Reweight scoretype="hbond_sr_bb" weight="2.0" />
         <Reweight scoretype="hbond_lr_bb" weight="2.0" />
         <Reweight scoretype="rama_prepro" weight="0.45" />
         <Reweight scoretype="omega" weight="0.4" />
         <Reweight scoretype="p_aa_pp" weight="0.6" />
     </ScoreFunction>
</SCOREFXNS>

2. RESIDUE_SELECTORS:氨基酸残基的选择。

• 选择phi小于或大于零的区域，用L-和d-氨基酸进行设计；
• 我们需要选择距离每个残基的两个氨基酸之外的残基来计算内部氢键；
• 使用Unsat选择不满的NH基团，之后会将这些选择好的氨基酸应用到下面的动作中。

<RESIDUE_SELECTORS>
    #Define different residue selections to be used later in the script 
    <Chain name="chainA" chains="A"/>
    <Chain name="sel_A" chains="1"/>
    
    #select the regions with phi less than or greater than zero for designing with L- and D-amino acids
    <Phi name="posPhi" select_positive_phi="true" />
    <Phi name="negPhi" select_positive_phi="false" />
    <Bin name="ProBin" bin="LPRO" bin_params_file="PRO_DPRO" />
    <Bin name="DProBin" bin="DPRO" bin_params_file="PRO_DPRO" />
    
    #we need to select residues that are outside +/-2 amino acids of each residue to calculate the internal h-bonds
    <Index name="hbonds_to_1" resnums="4,5,6,7,8"/>
    <Index name="hbonds_to_2" resnums="5,6,7,8,9"/>
    <Index name="hbonds_to_3" resnums="6,7,8,9,10"/>
    <Index name="hbonds_to_4" resnums="1,7,8,9,10"/>
    <Index name="hbonds_to_5" resnums="1,2,8,9,10"/>
    <Index name="hbonds_to_6" resnums="1,2,3,9,10"/>
    <Index name="hbonds_to_7" resnums="1,2,3,4,10"/>
    <Index name="hbonds_to_8" resnums="1,2,3,4,5"/>
    <Index name="hbonds_to_9" resnums="2,3,4,5,6"/>
    <Index name="hbonds_to_10" resnums="3,4,5,6,7"/>

    #select the residues with unsatisfied NH grups
    <Unsat name="unsat_amines" scorefxn="SFXN_STD" check_acceptors="false" consider_mainchain_only="false"/>
    <And name="ProBin_unsat" selectors="ProBin,unsat_amines" /> 
    <And name="DProBin_unsat" selectors="DProBin,unsat_amines" />   
</RESIDUE_SELECTORS>

3. SimpleMetrics:是为了返回pose的度量值

<SIMPLE_METRICS>
    <SelectedResidueCountMetric name="count_unsats" residue_selector="unsat_amines" />
</SIMPLE_METRICS>

4. TASKOPERATIONS:控制侧链的组装和设计。

• 前面的几行没什么特殊的，下面的四行是应用了限制性文件，将肽的选定区域的应用允许的残基。

<TASKOPERATIONS>
    <InitializeFromCommandline name="init"/>
    <IncludeCurrent name="current"/>
    <LimitAromaChi2 name="limitchi2" chi2max="110" chi2min="70" include_trp="True" />
    
    ### expanded chi1 and chi2 rotamer libraries
    <ExtraRotamersGeneric name="ex1_ex2" ex1="1" ex2="1"/>
    <ExtraRotamersGeneric name="ex1" ex1="1"/>
    
    #Add the allowed residues from resfile that apply to selected regions of the peptide
    <ReadResfile name="l_res" filename="l_res.txt" selector="negPhi"/>
    <ReadResfile name="d_res" filename="d_res.txt" selector="posPhi"/>
    <ReadResfile name="lpro_res" filename="pro.txt" selector="ProBin_unsat"/>
    <ReadResfile name="dpro_res" filename="dpro.txt" selector="DProBin_unsat"/>
    
</TASKOPERATIONS>

5. FILTERS: 对结果进行评估，抛弃一些不合理的结构。

• 首先计算多肽内部氨基酸的氢键数量大于2作为过滤的第一个条件；
• 删除氢键受体的氢键数量超过允许数量的设计的结果；
• 删除含有超过5个N-甲基化氨基酸的设计。

<FILTERS>
    #calculation of internal h-bonds in the peptide
    <HbondsToResidue name="hbonds1" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="1" residue_selector="hbonds_to_1"/>
    <HbondsToResidue name="hbonds2" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="2" residue_selector="hbonds_to_2"/>
    <HbondsToResidue name="hbonds3" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="3" residue_selector="hbonds_to_3"/>
    <HbondsToResidue name="hbonds4" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="4" residue_selector="hbonds_to_4"/>
    <HbondsToResidue name="hbonds5" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="5" residue_selector="hbonds_to_5"/>
    <HbondsToResidue name="hbonds6" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="6" residue_selector="hbonds_to_6"/>
    <HbondsToResidue name="hbonds7" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="7" residue_selector="hbonds_to_7"/>
    <HbondsToResidue name="hbonds8" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="8" residue_selector="hbonds_to_8"/>
    <HbondsToResidue name="hbonds9" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="9" residue_selector="hbonds_to_9"/>
    <HbondsToResidue name="hbonds10" partners="0" energy_cutoff="-0.25" backbone="1" bb_bb="1" sidechain="0" residue="10" residue_selector="hbonds_to_10"/> 
    <CombinedValue name="total_hbonds" threshold="-2.0"> #change it to the -1*number of hydrogens bonds required

        <Add filter_name="hbonds1" factor="-0.5"/>
        <Add filter_name="hbonds2" factor="-0.5"/>
        <Add filter_name="hbonds3" factor="-0.5"/>
        <Add filter_name="hbonds4" factor="-0.5"/>
        <Add filter_name="hbonds5" factor="-0.5"/>
        <Add filter_name="hbonds6" factor="-0.5"/>
        <Add filter_name="hbonds7" factor="-0.5"/>
        <Add filter_name="hbonds8" factor="-0.5"/>
        <Add filter_name="hbonds9" factor="-0.5"/>
        <Add filter_name="hbonds10" factor="-0.5"/>
    </CombinedValue>

    #Remove designs with any residue with a hydrogen bond acceptor that makes more than allowed number of h-bonds
    <OversaturatedHbondAcceptorFilter name="oversat" scorefxn="SFXN_STD" max_allowed_oversaturated="0" consider_mainchain_only="false" confidence="1" />
    
    #Remove designs with more than 5 N-methylated amino acids   
    <ResidueCount name="nme_count" max_residue_count="5" include_property="N_METHYLATED" confidence="1" />
    <ResidueCount name="nme_count2" max_residue_count="5" include_property="N_METHYLATED" confidence="1" />
    

    <SimpleMetricFilter name="C6s" metric="count_unsats" cutoff="0" comparison_type="lt_or_eq"/>

</FILTERS>

6. MOVERS: 对pose进行操作。

• 建立多聚甘氨酸主链；
• N端和C端之间的环化和肽键约束；
• 设置聚甘氨酸链的初始扭转；
• 基于内部氢键和过饱和的过滤条件来接受或拒绝GenKIC的结果；
• 使用GenKIC对多肽进行环化，里面的很多参数可参考；添
• 加约束条件来设置一些氨基酸的最小数量和最大数量；
• 第一轮的fastdesign：高加权氢键得分条件，以有利于内部背键到主链氢键的设计；
• 计算不满意的氨基酸残基，并进行N-甲基化；
• 第二轮fastdesign：利用上面的STD打分函数进行设计；
• 第三轮fastdesign：利用Cartesian打分函数进行设计，三种fastedesign的设计方法通过应用不同权重的打分函数来得到不同侧重的设计结果。

<MOVERS>
    #Construct the intial poly glycine chain
    <PeptideStubMover name="intial_stub" reset="true">
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />
        <Append resname="GLY" />        
    </PeptideStubMover>
    #Declare the peptide bond and cyclization constraints between the N- and C-terminal 
    <PeptideCyclizeMover name="cyclize_peptide" />

    #Set the intial phi/psi/omega torsions for the poly glycine chain   
    <SetTorsion name="torsion1">
        <Torsion residue="ALL" torsion_name="omega" angle="180.0" />
        <Torsion residue="1,2,3,4,5,6,7,8,9,10" torsion_name="rama" angle="rama_biased" custom_rama_table="flat_symm_dl_aa_ramatable"/>
    </SetTorsion>   

    #Set some filtering criteria based on internal h-bonds and oversaturation to accept or reject the solution from GenKIC 
    <ParsedProtocol name="preselection_pp">
        <Add filter="total_hbonds"/>
        <Add filter="oversat"/>
        <Add mover="cyclize_peptide" />
    </ParsedProtocol>   
    
    #GenKIC method for cyclizing the peptide
    <GeneralizedKIC name="genkic1" closure_attempts="1000" selector="lowest_energy_selector" pre_selection_mover="preselection_pp" stop_when_n_solutions_found="1" stop_if_no_solution="500" selector_scorefunction="SFXN_hbond_bb" >       
        <AddResidue res_index="7" />
        <AddResidue res_index="8" />
        <AddResidue res_index="9" />
        <AddResidue res_index="10" />
        <AddResidue res_index="1" />
        <AddResidue res_index="2" />
        <AddResidue res_index="3" />
        <AddResidue res_index="4" />
        <AddResidue res_index="5" /> #we are leaving resn 6 out of loop to be closed    
        <SetPivots atom1="CA" atom2="CA" atom3="CA" res1="7" res2="2" res3="5" />    
        <CloseBond prioratom_res="10" prioratom="CA" res1="10" atom1="C" res2="1" atom2="N" followingatom="CA" followingatom_res="1"  angle1="121.69997" angle2="116.199993" bondlength="1.328685" torsion="180.0" randomize_flanking_torsions="false" />   
        <AddPerturber effect="set_dihedral">
            <AddAtoms atom1="C" res1="10" res2="1" atom2="N" />
            <AddValue value="180.0" />
        </AddPerturber>        
        <AddPerturber effect="randomize_alpha_backbone_by_rama" custom_rama_table="flat_symm_dl_aa_ramatable">
            <AddResidue index="7" />
            <AddResidue index="8" />
            <AddResidue index="9"/>
            <AddResidue index="10"/>
            <AddResidue index="1"/>
            <AddResidue index="2"/>
            <AddResidue index="3"/>
            <AddResidue index="4"/>
            <AddResidue index="5"/>
        </AddPerturber>
        <AddFilter type="loop_bump_check" />
        <AddFilter type="rama_prepro_check" residue="7" rama_cutoff_energy="2" />
        <AddFilter type="rama_prepro_check" residue="2" rama_cutoff_energy="2" />
        <AddFilter type="rama_prepro_check" residue="5" rama_cutoff_energy="2" />        
    </GeneralizedKIC>

    #Add constraints to set min and max number of some amino acids
    <AddCompositionConstraintMover name="global" filename="comp/global_preferences.comp" selector="sel_A" />
    <AddCompositionConstraintMover name="pro_comp" filename="comp/pro.comp" selector="sel_A" />
    
    #sequence design with upweighted hydrogen bonding score terms to favor internal backbon-to-backbone hydrogen bonds
    <FastDesign name="fdesign_hbond_upweighted" scorefxn="SFXN_STD_high_hbond" repeats="1" task_operations="limitchi2,ex1_ex2,ex1,d_res,dpro_res,l_res,lpro_res" ramp_down_constraints="False"  >
        <MoveMap name="fdesign_step1_mm">
            <Chain number="1" chi="1" bb="1"/>
        </MoveMap>
    </FastDesign>

    #Calculate the number of unsats in the peptide
    <RunSimpleMetrics name="run_metrics" metrics="count_unsats" prefix="end_unsat" />

    #mutate the selected residues to their N-methyl versions
    <ModifyVariantType name="N_methylation" add_type="N_METHYLATION" residue_selector="unsat_amines" />
    
    #design amino acid sequence with STD sfxn
    <FastDesign name="fdesign" scorefxn="SFXN_STD" repeats="2" task_operations="limitchi2,ex1_ex2,ex1,d_res,l_res" ramp_down_constraints="False"  >
        <MoveMap name="fdesign_step2_mm">
            <Chain number="1" chi="1" bb="1"/>
        </MoveMap>
    </FastDesign>
    
    #design amino acid sequence with Cartesian sfxn
    <FastDesign name="fdesign_cart" scorefxn="SFXN_CART" repeats="1" task_operations="limitchi2,ex1_ex2,ex1,d_res,l_res" ramp_down_constraints="False" cartesian="True" >
        <MoveMap name="fdesign_step3_mm">
            <Chain number="1" chi="1" bb="1"/>
        </MoveMap>
    </FastDesign>
</MOVERS>

7. PROTOCOLS:将上述方法组合对pose进行操作；

<PROTOCOLS>
    <Add mover="intial_stub" />
    <Add mover="torsion1" />
    <Add mover="cyclize_peptide" />
    <Add mover="genkic1" />
    <Add mover="cyclize_peptide" />
    <Add mover="global" />
    <Add mover="fdesign_hbond_upweighted" />
    <Add filter="total_hbonds" />
    <Add filter="oversat" />
    <Add mover="N_methylation" />
    <Add filter="nme_count" />
    <Add mover="pro_comp" />
    <Add mover="fdesign" />
    <Add mover="N_methylation" />
    <Add mover="fdesign_cart" />
    <Add mover="cyclize_peptide" />
    <Add mover="run_metrics" />

</PROTOCOLS>

PROTOCOLS

开始 -> 选择宏环大小 -> 构建线性聚甘氨酸骨架 -> 设置N-to-C端环化的距离、角度和二面角约束 -> 随机选择一个残基作为“锚定残基”，另外三个残基作为“支点残基” -> 将锚定残基的二面角随机选择一个值，将支点残基的二面角从对称拉曼表中随机选择一个值 -> 使用GeneralizedKIC从线性多肽中识别环肽 -> 对于非支点和非锚定残基，从对称拉曼表中随机选择二面角 -> 对于支点残基，使用运动学闭合算法计算二面角，找到给出N-to-C环肽骨架的二面角组合 -> 定义闭合的条件 ->选择最低能量的GeneralizedKIC返回的环肽 ->使用Rosetta FastDesign 对环肽进行设计 -> 将“不满足”骨架NH基团的氨基酸突变为其N-甲基化 -> 对生成的模型进行过滤 -> 使用simple_cycpep_predict进行下一步的结构预测 -> 选择较低能量的构象作为结果的输出。

8.OUTPUT: 输出

<OUTPUT scorefxn="SFXN_CART"/>

策略2：基于变构透膜的方式设计通过设计存在两个种构象状态的多肽

使其在透膜过程中可以通过变构状态，将原本不饱和的极性基团自发形成氢键。从而形成可以透膜的“疏水态”。
3种设计方法：

1. 将策略1种生成的大量多肽进行大批量验证conformational energy landscape， 将存在2个或以上低能折叠态的序列选出。*（两个状态的能量差值在5 kcal/mol以内）
2. 对于折叠态大于5 kcal/mol的case，使用Rosetta multi-state的设计方法同时优化两个状态的序列，使得这些序列结构的能量相似。
3. 从晶体结构出发，引入不稳定的突变，并重复方法1的过程。

from rosetta import *
import glob, sys, math
import operator
import numpy as np
import numpy, os
from numpy import *
from pyrosetta import *
from rosetta.protocols.cyclic_peptide import DeclareBond
db = DeclareBond()


pyrosetta.init('-ex1 -ex2 -ex3 -ex4 -mute all -score:symmetric_gly_tables true -bbdep_omega false')#



class MSD:
    def __init__(self):
        self.chm = rosetta.core.chemical.ChemicalManager.get_instance()
        self.rts = self.chm.residue_type_set( 'fa_standard' )
        
        #State 1
        self.A = pose_from_pdb( sys.argv[1].strip())
        #State 2
        self.B = pose_from_pdb( sys.argv[2].strip())
                                            
        self.sfxn = create_score_function('ref2015')
        
        self.movemap = MoveMap()
        self.movemap.set_bb(False)
        self.movemap.set_chi(True)

        self.minmover = pyrosetta.rosetta.protocols.minimization_packing.MinMover(self.movemap, self.sfxn, "dfpmin_armijo", 0.01, True)#"dfpmin_armijo

        #Which residues to design. Replace res_to_design_string with a comma deliminated list of residue positions you wish to be mutated during the search. For example for a nine residue peptide wher you wanted to explore mutations at all positions put [1,2,3,4,5,6,7,8,9]
        self.res_to_design = res_to_design_string

        #Packer pallete at each position
        #replace allowed_mut_string with a comma deliminated list of integers ranging that signify which types of residues should be sampled at a position. 1 = no-N-methyl L-residue, 2 = no-N-methyl D-residue, 3 = N-methyl L residue, 4 = N-methyl D residue. For example for the peptide JC_9_99 a list of [3,4,2,1,2,1,2,1,1] is used
        self.allowed_mut = allowed_mut_string

        self.n_seq = 100

        self.d_all = {}
        self.d_all[1] = ['LEU', 'ALA', 'VAL', 'PHE', 'ILE']
        self.d_all[2] = ['DLEU', 'DALA', 'DVAL', 'DPHE', 'DILE']
        self.d_all[3] = ['LEU:N_Methylation', 'ALA:N_Methylation', 'PRO', 'VAL:N_Methylation', 'PHE:N_Methylation']
        self.d_all[4] = ['DLEU:N_Methylation', 'DALA:N_Methylation', 'DPRO', 'DVAL:N_Methylation', 'DPHE:N_Methylation']

        self.sequences = []

        while len(self.sequences) < self.n_seq:
            t_seq = []
            pose_A = Pose()
            pose_B = Pose()
            for i, j in enumerate(self.res_to_design):
                ri = random.randint(0, len(self.d_all[self.allowed_mut[i]]))
                t_seq.append(self.d_all[self.allowed_mut[i]][ri])
            self.sequences.append([pose_A,pose_B,t_seq, 0.00, 0.00, 0.00])
                
        #if you want to seed a starting sequence you can add it here
        #To do this replace seed_sequence_string with a comma deliminated list of the three/four character codes for the residues in a sequence. For example, for JC_9_99 the list ['PRO', 'DLEU', 'DLEU', 'PHE', 'DPRO', 'PHE', 'DLEU', 'LEU', 'VAL'] is provided
        self.sequences[0][2] = seed_sequence_string 

    def scoring(self):
        for i_s, seq in enumerate(self.sequences):
            self.tA = Pose(self.A)
            self.tB = Pose(self.B)
            for i, res in enumerate(self.res_to_design):
                self.tA.replace_residue(res, rosetta.core.conformation.ResidueFactory.create_residue( self.rts.name_map(seq[2][i])), True)
                self.tB.replace_residue(res, rosetta.core.conformation.ResidueFactory.create_residue( self.rts.name_map(seq[2][i])), True)


            db.set(seq_len, 'C', 1, 'N', True)
            db.apply(self.tA)
            db.set(seq_len, 'C', 1, 'N', True)
            db.apply(self.tB)

            task = rosetta.core.pack.task.TaskFactory.create_packer_task(self.tA)
            task.initialize_extra_rotamer_flags_from_command_line()
            task.restrict_to_repacking()
            packer = pyrosetta.rosetta.protocols.minimization_packing.PackRotamersMover(self.sfxn, task)
            
            packer.apply(self.tA)
            self.minmover.apply(self.tA)


            task = rosetta.core.pack.task.TaskFactory.create_packer_task(self.tB)
            task.initialize_extra_rotamer_flags_from_command_line()
            task.restrict_to_repacking()
            packer = pyrosetta.rosetta.protocols.minimization_packing.PackRotamersMover(self.sfxn, task)
            packer.apply(self.tB)
            self.minmover.apply(self.tB)
            

            eA = self.sfxn(self.tA)
            eB = self.sfxn(self.tB)
            seq[0]=self.tA
            seq[1]=self.tB

            if eA < 10.0 and eB < 10.0 and math.fabs(eA-eB) < 6.0:
                self.sequences[i_s][-1] = (math.fabs(eA - eB)*-5)-(eA+eB)  #-(0.5*(eA + eB)))
            else:
                self.sequences[i_s][-1] = -1000.0

            self.sequences[i_s][-2] = eB
            self.sequences[i_s][-3] = eA

            if i_s == 0:
                self.tA.center()
                self.tB.center()
                self.tA.dump_pdb('state_A.pdb')
                self.tB.dump_pdb('state_B.pdb')

        self.sequences.sort(key=operator.itemgetter(-1), reverse=True)

    def propogate(self):
        seq_l = []
        for s in self.sequences:
            seq_l.append(s[2])
        for i in range(0, int(self.n_seq/2)-1):
            r1 = random.randint(0, len(self.res_to_design)-1)
            r2 = random.randint(0, len(self.d_all[self.allowed_mut[r1]]))

            t_seq = list(self.sequences[i][2])
            t_seq[r1] = self.d_all[self.allowed_mut[r1]][r2]
            while t_seq in seq_l:
                r1 = random.randint(0, len(self.res_to_design)-1)
                r2 = random.randint(0, len(self.d_all[self.allowed_mut[r1]]))
                t_seq = list(self.sequences[i+int(self.n_seq/2)][2])
                t_seq[r1] = self.d_all[self.allowed_mut[r1]][r2]

            pose_A = Pose()
            pose_B = Pose()
            self.sequences[i+int(self.n_seq/2)] = [pose_A, pose_B, t_seq, 0.00, 0.00, 0.00]

            seq_l = []
            for s in self.sequences:
                seq_l.append(s[2])
                
    def print_seqs(self, i, name):
        f = open('msd'+name+'_' + str(i) + '_out.txt', 'w')
        for s in self.sequences:
            f.write(str(s[2]).replace('[', '|').replace(']', '|')+str(s[3])+'|'+str(s[4])+'|'+str(s[5])+'\n')
        f.close()     

    def dump_best(self):
        self.sequences[0][0].dump_pdb("best_sequence_A.pdb")
        self.sequences[0][1].dump_pdb("best_sequence_B.pdb")

APP_MSD = MSD()

#number of generations to run
generations = 1000

for i in range(0, generations):
    print( i )
    APP_MSD.scoring()
    APP_MSD.print_seqs(i, 'test')
    APP_MSD.propogate()
    APP_MSD.dump_best()

代码逻辑非常简单，对环肽的部分位置极性氨基酸的替换，替换成
['LEU', 'ALA', 'VAL', 'PHE', 'ILE']
['DLEU', 'DALA', 'DVAL', 'DPHE', 'DILE']
['LEU:N_Methylation', 'ALA:N_Methylation', 'PRO', 'VAL:N_Methylation', 'PHE:N_Methylation']
['DLEU:N_Methylation', 'DALA:N_Methylation', 'DPRO', 'DVAL:N_Methylation', 'DPHE:N_Methylation']，并计算替换前后的得分，最后得到结果。

Comment

多肽、蛋白质和基因类生物大分子药物在许多疾病的治疗中发挥着重要的作用，但由于其尺寸较大 和低亲脂特性，导致其不能有效穿过细胞膜，使得这些具有治疗价值的生物大分子在医学等领域的 应用受到了极大限制。这篇文章为解决以上问题的解决提供了新的思路。但环肽采样方法上并没有特别重大的改进，与以前的Anchor extension的方法大同小异。

蛋白的环肽binder设计

Anchor extension: a structure-guided approach to design cyclic peptides targeting enzyme active sites

这篇文章的主要内容是关于循环肽的计算设计，旨在通过扩展锚点来针对蛋白质界面，从而设计出高亲和力的蛋白质-肽界面。文章中介绍了一种计算方法，称为“锚点扩展”，该方法通过将肽链围绕非规范氨基酸残基锚点进行扩展，以便更好地固定在蛋白质功能位点上。作者使用组合动力学环路闭合方法在Rosetta软件中构建了一个环状肽，以将锚点放置在具有结合能力的方向上，并通过提供计算设计期间引入的额外相互作用来增强其与目标的结合。文章还介绍了作者使用该方法设计出的循环肽，这些循环肽可以抑制组蛋白去乙酰化酶2和6（HDAC2和HDAC6），并且比原始锚点具有更高的亲和力。这些结果突显了高亲和力蛋白-肽界面的全新设计方法。

Computationally designed peptide macrocycle inhibitors of New Delhi metallo-β-lactamase 1

这篇文章的主要内容是关于计算设计的肽环抑制剂，用于抑制新德里金属-β-内酰胺酶1（NDM-1）的。NDM-1是一种细菌酶，可以降解β-内酰胺类抗生素，而抗生素耐药性的崛起需要新的治疗方法。这篇文章介绍了一种计算方法，用于设计由L-和D-氨基酸混合构建的肽环，这些肽环能够结合和抑制治疗感兴趣的目标。该方法明确考虑了肽是否倾向于有一个结合竞争构象，我们发现这种倾向性可以预测七个设计的NDM-1抑制肽的实验观察到的IC50值的排名顺序。我们能够确定三个设计的抑制剂与NDM-1形成复合物的X射线晶体结构，在其中三个中，肽的构象非常接近计算设计模型。在三个结构中，与NDM-1的结合模式也非常类似于设计模型，在第三个结构中，我们观察到一种可能由于设计形状内部对称性与靶标的灵活性相结合而产生的替代结合模式。虽然在强大地预测靶向脊椎变化、结合模式以及突变对结合亲和力的影响方面仍然存在挑战，但我们为设计有序、具有结合竞争力的肽环的方法应该具有广泛适用性提供了方法。

Target-templated de novo design of macrocyclic d-/l-peptides: discovery of drug-like inhibitors of PD-1

主要的设计方案

总体来说，首先要选择的PPI中对结合力贡献较大hotspot位点，并以此hotspot作为构建环肽的锚点。在这个锚点中，同时保持目标蛋白作为输入，通过在锚点的N端和C端分别添加多个个GLY残基来生成骨干，并保持氨基酸的数量在8-10个左右。甘氨酸是唯一自然存在的氨基酸，探索拉马钱德拉图的两侧，聚甘氨酸骨架是设计异手性拓扑结构的合适模型。然后多聚甘氨酸在GenKIC算法的加持下进行采样，生成成百上千个构象，之后利用ddG、interface sasa、氢键数量等指标对构象进行过滤筛选。具体的实现方法与上面的穿膜环肽设计的方法类似。

其中不同的是，因为hotspot位点的存在，我们可以在设计的过程中让这些hotspot不发生改变，具体的实现方法为：

<TASKOPERATIONS>
<OperateOnResidueSubset name="no_design" selector="nochange_res">
        <RestrictToRepackingRLT/>
</OperateOnResidueSubset>
</TASKOPERATIONS>

其次是聚甘氨酸环的构建，只需要在锚点的前面及后面添加氨基酸即可：

<MOVERS>
<PeptideStubMover name="intial_stub" reset="0">
    <Prepend resname="GLY" anchor_rsd="351" repeat="4"/>
    <Append resname="GLY" repeat="3"/>
</PeptideStubMover>
</MOVERS>

然后就是添加能量、形状互补等过滤条件:

<FILTERS>
    #total energy score
    <ScoreType name="score" scorefxn="SFXN_STD" score_type="total_score" threshold="-650" confidence="0" />
    #filter to discard very bad structures after genKIC
    <ScoreType name="score_low" scorefxn="SFXN_STD" score_type="total_score" threshold="0"/>
    #interface shape complementarity
    <ShapeComplementarity name="sc_filter" verbose="0" min_sc="0.6" write_int_area="1" jump="1" confidence="1"/>
    #computes ddG with and without minimizing
    <Ddg name="ddg" threshold="-1" jump="1" repeats="5" repack="1" confidence="1" scorefxn="SFXN_STD"/>
    <Ddg name="ddg_norepack" threshold="-1" jump="1" repeats="1" repack="0" confidence="1" scorefxn="SFXN_STD"/>
    #computes the interface sasa (total, hydrophobic, and polar)
    <Sasa name="interface_buried_sasa" confidence="0"/>
    <Sasa name="interface_hydrophobic_sasa" confidence="0" hydrophobic="True"/>
    <Sasa name="interface_polar_sasa" confidence="0" polar="True"/> 
</FILTERS>

Comment

前两篇文章都是baker实验室的成果，但是都是关于酶的，其中设计的时候要配有金属离子作为催化。后一篇是利用PPI的结果来进行设计，但是结果中的结合力仅在um级别，效果一般。

其他

Cyclic peptide structure prediction and design using AlphaFold | bioRxiv

基于AlphaFold的环肽结构预测和设计方法：这是一种利用深度学习模型AlphaFold进行环肽结构预测和设计的方法，通过修改相对位置编码的输入来扩展AlphaFold，使其能够用于环肽的结构预测，并测试了这些变化在预测PDB中可用的环肽结构方面的准确性。然后，作者报告了一种使用AlphaFold重新设计大环主链序列的方法，以提高它们折叠成设计结构的倾向。最后，他们描述了一种从头幻觉设计新环肽的方法，并使用它来列举7-13元环肽丰富多样性。

但是我是用之后发现这种环肽预测受限于训练集中无非天然氨基酸，所以无论是在序列设计还是从头幻觉上，都无法得到包含D型或者其他非天然氨基酸。这对后面的使用及优化产生了限制。

其他文章

• Accurate de novo design of hyperstable constrained peptides
• The emerging role of computational design in peptide macrocycle drug discovery
• Cyclic peptide drugs approved in the last two decades (2001–2021)
• Sub-angstrom accuracy in protein loop reconstruction by robotics-inspired conformational sampling
• Comprehensive computational design of ordered peptide macrocycles

参考

Anchor extension：基于热点残基拓展的环肽设计方法 - 知乎 (zhihu.com)
Rosetta simple_cycpep_predict: 预测环肽结构 - 知乎 (zhihu.com)
高精度从头设计透膜环肽分子 - 知乎 (zhihu.com)