{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rosetta Scripts in PyRosetta (Advanced)\n",
    "\n",
    "@Author: Jian Huang\n",
    "\n",
    "@E-mail: jian.huang@xtalpi.com"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过第一节的学习,理论而言,只要我们掌握了RS的各个tag的含义和定义,就可以进行非常复杂且用户定制化的protocols。本节的内容提供关于Minimization和packing的RS高级设置,让我们可以与之前学习过的TaskOperation和MoveMap知识进行联动。并进一步了解RS的编写和用法。\n",
    "\n",
    "</br>\n",
    "\n",
    "**写在开头:一点回顾**\n",
    "\n",
    "*关于packing*\n",
    "\n",
    "Rosetta中的packing过程:在一个蛋白pose构象中,由于每个位置的氨基酸都具有其本身的rotamer分布,想要去找到哪种组合的rotamer会让体系能量最低的过程。rotamer其实就是氨基酸侧链的构象分布,在rosetta中使用的rotamer library是在PDB库中采取那些精度特别高(1.8 A)的晶体结构进行统计的结果,虽然rotamer在空间中分布是连续的,但是概率出现最大的几个位置被选出来作为能量优势构象,然后连续的rotamer分布面就被几个分离的概率较大的rotamer作为代表而简化处理了。\n",
    "\n",
    "即便如此,在N个氨基酸组成的蛋白pose中,假设每个氨基酸由M种rotamer,那么搜索的空间就是M的N次方,可以看到随着N的增加,搜索空间将无比巨大。Rosetta的Packer就是使用一种蒙特卡洛的随机搜索算法,以接受能够使构象能量降低的构象变化或以一定概率接受能量升高的构象变化,目的使寻找这些N个氨基酸的哪些rotamer的组合,可以给出能量最优点的构象。从上述描述中可以看到,该搜索算法使随机的,它只能够保证运行时间结束后,得到的构象能量会比之前的小,但却并不一定使全局能量最优构象。(“尽量去搜索最优点,但不保证能够得到最优点”)\n",
    "\n",
    "*关于minimizer*\n",
    "\n",
    "第一节使用的MinMover是一个使用最朴实的Minimization的例子,其在底层调用了rosetta中的底层算法——minimizer。\n",
    "\n",
    "之前在学习minimizer的时候,提到所有的minimizer的操作使用了rosetta的**Movemap**来确定构象的自由度范围。\n",
    "\n",
    "MoveMap在minimizer中是可以允许用户精准控制能量最小化限制的操作。MoveMap可以允许用户精确定义是否移动backbone(BB)的扭转角度变化(也就是phi角和psi角度)以及侧链扭转角变化(chi角)。其实从这里的描述也可以看到,minimizer与packer十分相似,相比于packer是去寻找哪种组合的rotamer会让能量最低,minimizer也是在空间中搜索哪一些组合的扭转角会让构象的能量最低,只是他们用的搜索算法和优化能量采取的研究对象不同(一个是rotamer,一个是扭转角)罢了。同时也可以看到,如果我们先用不做design的packer优化一下能量,但由于我们采用的rotamer库使用的是分离的侧链构象,近似为理想化的情况,此时做minimizer就很有必要了,因为它可以对构象的扭转角进行一些微调,从而采样到那些并不是理想化rotamer的情况的构象。\n",
    "\n",
    "值得注意的是,尽管我们可以在MoveMap中禁止某些侧链扭转角的变化,也并不代表该残基在空间中不会发生位移。其实这也很好理解,蛋白质链都是依次序连接在一起,如果上游的某个扭转角改变了,势必会影响它接下来的原子在空间中的位置,尽管后续的扭转角都没有变化。在rosetta中,其实是采取了一个叫FoldTree的概念来描述这种连接的关系。(关于FoldTree,请参考本教程相关章节)\n",
    "\n",
    "在使用MoveMap进行minimization的时候,我们必须注意FoldTree。rosetta中默认的FoldTree定义第一条链A的第一个残基作为根残基(root),依次后面称第二个残基是子代残基1,第三个残基是子代残基2,......等等;如果存在多条链,jump会转接到另一条链的第一个残基,而后重复上一步操作,形成完成的FoldTree。任何对于指定残基的自由度的变化都会影响到它在Foldtree中下游子代。\n",
    "\n",
    "在MoveMap中可以自由定义各种自由度。与MinMover相似,packing操作在RS中其实也是通过定义mover来进行,其在底层会调用packer以及控制packer的**TaskOperation**进行精确的packer操作。\n",
    "\n",
    "***\n",
    "\n",
    "我们应该注意,**Movemap**的定义具有不可交换性,某一个残基的Movemap定义以其最后一次出现的定义为准;而**TaskOperation**具有可交换性,交换一个packing过程的各个TaskOperation定义等同。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. RS in Pyrosetta:精准定义MoveMap的Minimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在第一节中,我们给出了官网上MinMover的定义范式,在```<MinMover> </MinMover>```标签内含有允许用户定义的亚标签——```<MoveMap></MoveMap>```。\n",
    "\n",
    "我们可以利用该亚标签对MinMover进行细致、精确地定义,以满足我们对某一pose的rosetta采样需求。\n",
    "\n",
    "Example:\n",
    "\n",
    "```\n",
    "<SCOREFXNS>\n",
    "        <ScoreFunction name=\"ref_2015\" weights=\"ref2015\"/>\n",
    "</SCOREFXNS>\n",
    "<MOVERS>\n",
    "    <MinMover name=\"min_torsion\" scorefxn=\"ref_2015\" chi=\"true\" bb=\"1\" cartesian=\"F\" >\n",
    "        <MoveMap name=\"min_torsion_mm\">                         \n",
    "            <Span begin=\"1\" end=\"999\" chi=\"false\" bb=\"false\" /> \n",
    "            <Span begin=\"1\" end=\"50\" chi=\"true\" bb=\"true\" />    \n",
    "            <Span begin=\"5\" end=\"10\" chi=\"true\" bb=\"false\" />   \n",
    "        </MoveMap>                                              \n",
    "    </MinMover>\n",
    "</MOVERS>\n",
    "<PROTOCOLS>\n",
    "    <Add mover=\"min_torsion\" />\n",
    "</PROTOCOLS>\n",
    "```\n",
    "\n",
    "过程:\n",
    "\n",
    "1. 先限制1~999号残基的BB和CHI移动;\n",
    "2. 然后允许1~55号残基的CHI和BB运动;\n",
    "3. 最后限制5~10号残基的BB不运动、CHI可以运动。\n",
    "\n",
    "NOTE: \n",
    "这里先定义了一个很宽的范围1 - 999,注意在Rosetta中即使我们现在正在模拟的蛋白长度没有这么多,也是允许的!只是其他的限制仍然还是在输入构象基础上进行限制的(比如1 - 50,5 - 10)。(假定我们模拟的蛋白长度是76,我们知道MoveMap遵循依次定义限制的规则,一个残基的自由度由最后一次出现它定义的地方指定,比如5 - 10号残基先是禁止所有运动,而后开放了CHI和BB的运动,最后又限制了BB运动,开放CHI运动——最终的结果就是5 - 10号允许CHI运动、不允许BB运动;而1 - 4和11 - 50具有BB和CHI自由度;51 - 76不允许任何移动)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyRosetta-4 2021 [Rosetta PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release 2021.31+release.c7009b3115c22daa9efe2805d9d1ebba08426a54 2021-08-07T10:04:12] retrieved from: http://www.pyrosetta.org\n",
      "(C) Copyright Rosetta Commons Member Institutions. Created in JHU by Sergey Lyskov and PyRosetta Team.\n",
      "\u001b[0mcore.init: {0} \u001b[0mChecking for fconfig files in pwd and ./rosetta/flags\n",
      "\u001b[0mcore.init: {0} \u001b[0mRosetta version: PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release r292 2021.31+release.c7009b3115c c7009b3115c22daa9efe2805d9d1ebba08426a54 http://www.pyrosetta.org 2021-08-07T10:04:12\n",
      "\u001b[0mcore.init: {0} \u001b[0mcommand: PyRosetta -ex1 -ex2aro -database /opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database\n",
      "\u001b[0mbasic.random.init_random_generator: {0} \u001b[0m'RNG device' seed mode, using '/dev/urandom', seed=-1031901739 seed_offset=0 real_seed=-1031901739 thread_index=0\n",
      "\u001b[0mbasic.random.init_random_generator: {0} \u001b[0mRandomGenerator:init: Normal mode, seed=-1031901739 RG_type=mt19937\n"
     ]
    }
   ],
   "source": [
    "from pyrosetta.rosetta.protocols.rosetta_scripts import *\n",
    "from pyrosetta import *\n",
    "init()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0mcore.chemical.GlobalResidueTypeSet: {0} \u001b[0mFinished initializing fa_standard residue type set.  Created 983 residue types\n",
      "\u001b[0mcore.chemical.GlobalResidueTypeSet: {0} \u001b[0mTotal time to initialize 0.67358 seconds.\n",
      "\u001b[0mcore.import_pose.import_pose: {0} \u001b[0mFile './data/1ubq_clean.pdb' automatically determined to be of type PDB\n"
     ]
    }
   ],
   "source": [
    "# 读入初始pose,并且拷贝一份\n",
    "original_pose = pose_from_pdb(\"./data/1ubq_clean.pdb\")\n",
    "pose = original_pose.clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 也可以将rosetta scripts直接在python中用字符串表示\n",
    "# 省去读入文件一步\n",
    "\n",
    "RS_mini = \\\n",
    "\"\"\"\n",
    "<ROSETTASCRIPTS>\n",
    "    <SCOREFXNS>    \n",
    "        <ScoreFunction name=\"ref_2015\" weights=\"ref2015\" />\n",
    "    </SCOREFXNS>\n",
    "    <RESIDUE_SELECTORS>\n",
    "    </RESIDUE_SELECTORS>\n",
    "    <TASKOPERATIONS>\n",
    "    </TASKOPERATIONS>\n",
    "    <FILTERS>\n",
    "    </FILTERS>\n",
    "    <MOVERS>\n",
    "        <MinMover name=\"min_torsion\" scorefxn=\"ref_2015\" chi=\"true\" bb=\"1\" cartesian=\"F\" >\n",
    "            <MoveMap name=\"min_torsion_mm\">                         \n",
    "                <Span begin=\"1\" end=\"999\" chi=\"false\" bb=\"false\" /> \n",
    "                <Span begin=\"1\" end=\"50\" chi=\"true\" bb=\"true\" />    \n",
    "                <Span begin=\"5\" end=\"10\" chi=\"true\" bb=\"false\" />   \n",
    "            </MoveMap>                                              \n",
    "        </MinMover>\n",
    "    </MOVERS>\n",
    "    <APPLY_TO_POSE>\n",
    "    </APPLY_TO_POSE>\n",
    "    <PROTOCOLS>\n",
    "        <Add mover=\"min_torsion\" />    \n",
    "    </PROTOCOLS>\n",
    "    <OUTPUT scorefxn=\"ref_2015\" />\n",
    "</ROSETTASCRIPTS>\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mGenerating XML Schema for rosetta_scripts...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mInitializing schema validator...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mValidating input script...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mParsed script:\n",
      "<ROSETTASCRIPTS>\n",
      "\t<SCOREFXNS>\n",
      "\t\t<ScoreFunction name=\"ref_2015\" weights=\"ref2015\"/>\n",
      "\t</SCOREFXNS>\n",
      "\t<RESIDUE_SELECTORS/>\n",
      "\t<TASKOPERATIONS/>\n",
      "\t<FILTERS/>\n",
      "\t<MOVERS>\n",
      "\t\t<MinMover bb=\"1\" cartesian=\"F\" chi=\"true\" name=\"min_torsion\" scorefxn=\"ref_2015\">\n",
      "\t\t\t<MoveMap name=\"min_torsion_mm\">\n",
      "\t\t\t\t<Span bb=\"false\" begin=\"1\" chi=\"false\" end=\"999\"/>\n",
      "\t\t\t\t<Span bb=\"true\" begin=\"1\" chi=\"true\" end=\"50\"/>\n",
      "\t\t\t\t<Span bb=\"false\" begin=\"5\" chi=\"true\" end=\"10\"/>\n",
      "\t\t\t</MoveMap>\n",
      "\t\t</MinMover>\n",
      "\t</MOVERS>\n",
      "\t<APPLY_TO_POSE/>\n",
      "\t<PROTOCOLS>\n",
      "\t\t<Add mover=\"min_torsion\"/>\n",
      "\t</PROTOCOLS>\n",
      "\t<OUTPUT scorefxn=\"ref_2015\"/>\n",
      "</ROSETTASCRIPTS>\n",
      "\u001b[0mcore.scoring.ScoreFunctionFactory: {0} \u001b[0mSCOREFUNCTION: \u001b[32mref2015\u001b[0m\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0mStarting energy table calculation\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0msmooth_etable: changing atr/rep split to bottom of energy well\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0msmooth_etable: spline smoothing lj etables (maxdis = 6)\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0msmooth_etable: spline smoothing solvation etables (max_dis = 6)\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0mFinished calculating energy tables.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/hbonds/ref2015_params/HBPoly1D.csv\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/hbonds/ref2015_params/HBFadeIntervals.csv\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/hbonds/ref2015_params/HBEval.csv\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/hbonds/ref2015_params/DonStrength.csv\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/hbonds/ref2015_params/AccStrength.csv\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/fd/all.ramaProb\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/fd/prepro.ramaProb\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/omega/omega_ppdep.all.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/omega/omega_ppdep.gly.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/omega/omega_ppdep.pro.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/omega/omega_ppdep.valile.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/P_AA_pp/P_AA\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/P_AA_pp/P_AA_n\n",
      "\u001b[0mcore.scoring.P_AA: {0} \u001b[0mshapovalov_lib::shap_p_aa_pp_smooth_level of 1( aka low_smooth ) got activated.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/P_AA_pp/shapovalov/10deg/kappa131/a20.prop\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0mStarting energy table calculation\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0msmooth_etable: changing atr/rep split to bottom of energy well\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0msmooth_etable: spline smoothing lj etables (maxdis = 6)\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0msmooth_etable: spline smoothing solvation etables (max_dis = 6)\n",
      "\u001b[0mcore.scoring.etable: {0} \u001b[0mFinished calculating energy tables.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/PairEPotential/pdb_pair_stats_fine\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/InterchainPotential/interchain_env_log.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/InterchainPotential/interchain_pair_log.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/EnvPairPotential/env_log.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/EnvPairPotential/cbeta_den.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/EnvPairPotential/pair_log.txt\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/EnvPairPotential/cenpack_log.txt\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mshapovalov_lib::shap_rama_smooth_level of 4( aka highest_smooth ) got activated.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/shapovalov/kappa25/all.ramaProb\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/avg_L_rama.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/avg_L_rama.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/sym_all_rama.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/sym_all_rama.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/sym_G_rama.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/sym_G_rama.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/sym_P_rama.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/sym_P_rama.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/avg_L_rama_str.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/avg_L_rama_str.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/sym_all_rama_str.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/sym_all_rama_str.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/sym_G_rama_str.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/sym_G_rama_str.dat.\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/rama/flat/sym_P_rama_str.dat\n",
      "\u001b[0mcore.scoring.ramachandran: {0} \u001b[0mReading custom Ramachandran table from scoring/score_functions/rama/flat/sym_P_rama_str.dat.\n",
      "\u001b[0mprotocols.jd2.parser.ScoreFunctionLoader: {0} \u001b[0mdefined score function \"ref_2015\" with weights \"ref2015\"\n",
      "\u001b[0mprotocols.minimization_packing.MinMover: {0} \u001b[0mOptions chi, bb: 1, 1 omega: 1\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mDefined mover named \"min_torsion\" of type MinMover\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mParsedProtocol mover with the following settings\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mAdded mover \"min_torsion\"\n"
     ]
    }
   ],
   "source": [
    "# it's gonna take a while for Rosetta to parse the RS\n",
    "xml = XmlObjects.create_from_string(RS_mini)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0m=======================BEGIN MOVER MinMover - min_torsion=======================\n",
      "\u001b[0mcore.select.residue_selector.ResidueSpanSelector: {0} \u001b[0m\u001b[1m[ WARNING ]\u001b[0m Residue span end designation 'Residue 999' is outside the Pose!\n",
      "\u001b[0mbasic.io.database: {0} \u001b[0mDatabase file opened: scoring/score_functions/elec_cp_reps.dat\n",
      "\u001b[0mcore.scoring.elec.util: {0} \u001b[0mRead 40 countpair representative atoms\n",
      "\u001b[0mcore.pack.dunbrack.RotamerLibrary: {0} \u001b[0mshapovalov_lib_fixes_enable option is true.\n",
      "\u001b[0mcore.pack.dunbrack.RotamerLibrary: {0} \u001b[0mshapovalov_lib::shap_dun10_smooth_level of 1( aka lowest_smooth ) got activated.\n",
      "\u001b[0mcore.pack.dunbrack.RotamerLibrary: {0} \u001b[0mBinary rotamer library selected: /opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database/rotamer/shapovalov/StpDwn_0-0-0/Dunbrack10.lib.bin\n",
      "\u001b[0mcore.pack.dunbrack.RotamerLibrary: {0} \u001b[0mUsing Dunbrack library binary file '/opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database/rotamer/shapovalov/StpDwn_0-0-0/Dunbrack10.lib.bin'.\n",
      "\u001b[0mcore.pack.dunbrack.RotamerLibrary: {0} \u001b[0mDunbrack 2010 library took 0.168718 seconds to load from binary\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0msetting status to success\n"
     ]
    }
   ],
   "source": [
    "protocol = xml.get_mover(\"ParsedProtocol\")\n",
    "protocol.apply(pose)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. RS in Pyrosetta:精准定义TaskOperation的packing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "</br>\n",
    "\n",
    "与Minimizer控制二面角自由度类似,packing过程通过TaskOperation控制packer的氨基酸design和packing的自由度。\n",
    "\n",
    "RS中对应经典的packer的Mover名为:PackRotamersMover"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "</br>\n",
    "\n",
    "同样地,我们可以从官方网页上找到定义PackRotamersMover的定义范式:\n",
    "(https://new.rosettacommons.org/docs/latest/scripting_documentation/RosettaScripts/Movers/movers_pages/PackRotamersMover)\n",
    "\n",
    "```\n",
    "<TASKOPERATIONS>\n",
    "    # Some TaskOperation definition here\n",
    "</TASKOPERATIONS>\n",
    "<PackRotamersMover name=\"(&string;)\" nloop=\"(1 &non_negative_integer;)\"\n",
    "        scorefxn=\"(&string;)\"\n",
    "        task_operations=\"(&task_operation_comma_separated_list;)\"\n",
    "        packer_palette=\"(&named_packer_palette;)\" />\n",
    "        \n",
    "```\n",
    ">nloop: Equivalent to \"-ndruns\".Number of complete packing runs before an output (best score) is produced.\n",
    ">\n",
    ">scorefxn: Name of score function to use\n",
    ">\n",
    ">task_operations: A comma-separated list of TaskOperations to use.\n",
    ">\n",
    ">packer_palette: A previously-defined PackerPalette to use, which specifies the set of residue types with which to design (to be pruned with TaskOperations).\n",
    "\n",
    "\n",
    "**注意**\n",
    "\n",
    "1. RS中有单独的```<TASKOPERATIONS></TASKOPERATIONS>```标签定义全局的TaskOperation内容。然后通过在PackRotamersMover的选项:task_operation控制需要使用的TaskOperations。\n",
    "\n",
    "2. PackRotamersMover中的task_operation选项可以直接删掉(定义为默认值),但是默认的PackRotamersMover会允许每一个位置使用20种的氨基酸的所有rotamer进行packing!\n",
    "\n",
    "***\n",
    "\n",
    "故此,我们还需要了解如何定义```<TASKOPERATIONS></TASKOPERATIONS>```。官方例子如下:\n",
    "\n",
    "```\n",
    "...\n",
    "<TASKOPERATIONS>\n",
    "  <ReadResfile name=\"rrf\"/>\n",
    "  <ReadResfile name=\"rrf2\" filename=\"resfile2\"/>\n",
    "  <PreventRepacking name=\"NotBeingUsedHereButPresenceOkay\"/>\n",
    "  <RestrictResidueToRepacking name=\"restrict_Y100\" resnum=\"100\"/>\n",
    "  <RestrictToRepacking name=\"rtrp\"/>\n",
    "  <OperateOnResidueSubset name=\"NoPackHelix\">\n",
    "      <SecondaryStructure ss=\"H\" />\n",
    "      <PreventRepackingRLT/>\n",
    "  </OperateOnResidueSubset>\n",
    "</TASKOPERATIONS>\n",
    "...\n",
    "<MOVERS>\n",
    "  <PackRotamersMover name=\"packrot\" scorefxn=\"sf\" task_operations=\"rrf,NoPackHelix,rtrp,restrict_Y100\"/>\n",
    "</MOVERS>\n",
    "...\n",
    "```\n",
    "\n",
    "\n",
    "Rosetta为了设定TaskOperation的便利,内置定义了多个独特的TASKOPERATIONS,例如PreventRepacking、RestrictResidueToRepacking等等。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "除此之外,TaskOperations还会控制侧链如何被采样的具体过程,rosetta中默认是基于rotamer的采样,但常常可能希望加入一些额外的rotamer的采样。例如,添加一些偏移标准rotamer库的因子,ExtraRotamersGeneric的TaskOperations允许我们控制rotamer采样的级别,一般来说在CHI1和CHI2中加入一些额外的采样比较有效,虽然可能会花费更长的时间进行packing。此外InitializeFromCommandline也允许使用-ex1 -ex2的选项进行rotamer的采样控制。\n",
    "\n",
    "```\n",
    "\n",
    "<ExtraRotamersGeneric name=\"(&string)\"\n",
    "ex1=\"(0 &boolean)\" ex2=\"(0 &boolean)\" ex3=\"(0 &boolean)\" ex4=\"(0 &boolean)\"\n",
    "ex1aro=\"(0 &boolean)\" ex2aro=\"(0 &boolean)\" ex1aro_exposed=\"(0 &boolean)\" ex2aro_exposed=\"(0 &boolean)\"\n",
    "ex1_sample_level=\"(7 &Size)\" ex2_sample_level=\"(7 &Size)\" ex3_sample_level=\"(7 &Size)\" ex4_sample_level=\"(7 &Size)\"\n",
    "ex1aro_sample_level=\"(7 &Size)\" ex2aro_sample_level=\"(7 &Size)\" ex1aro_exposed_sample_level=\"(7 &Size)\" ex2aro_exposed_sample_level=\"(7 &Size)\"\n",
    "exdna_sample_level=\"(7 &Size)\"\n",
    "extrachi_cutoff=\"(18 &Size)\"/>\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义本例中的RS\n",
    "RS_pack = \\\n",
    "\"\"\"\n",
    "<ROSETTASCRIPTS>\n",
    "        <SCOREFXNS>\n",
    "                <ScoreFunction name=\"r15\" weights=\"ref2015\" />\n",
    "        </SCOREFXNS>\n",
    "        <RESIDUE_SELECTORS>\n",
    "        </RESIDUE_SELECTORS>\n",
    "        <TASKOPERATIONS>\n",
    "                <RestrictToRepacking name=\"no_design\" />    \n",
    "                <ExtraRotamersGeneric name=\"extrachi\" ex1=\"1\" ex2=\"1\" ex1_sample_level=\"1\" ex2_sample_level=\"1\" />\n",
    "        </TASKOPERATIONS>\n",
    "        <FILTERS>\n",
    "        </FILTERS>\n",
    "        <MOVERS>\n",
    "                <PackRotamersMover name=\"pack1\" scorefxn=\"r15\" task_operations=\"no_design,extrachi\" />\n",
    "        </MOVERS>\n",
    "        <APPLY_TO_POSE>\n",
    "        </APPLY_TO_POSE>\n",
    "        <PROTOCOLS>\n",
    "                <Add mover=\"pack1\" />\n",
    "        </PROTOCOLS>\n",
    "        <OUTPUT scorefxn=\"r15\" />\n",
    "</ROSETTASCRIPTS>\n",
    "\"\"\"\n",
    "\n",
    "# 第一个TASKOPERATION:RestrictToRepacking 会禁止design;\n",
    "# 第二个TASKOPERATION:ExtraRotamersGeneric会允许额外rotamer的采样"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mGenerating XML Schema for rosetta_scripts...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mInitializing schema validator...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mValidating input script...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mParsed script:\n",
      "<ROSETTASCRIPTS>\n",
      "\t<SCOREFXNS>\n",
      "\t\t<ScoreFunction name=\"r15\" weights=\"ref2015\"/>\n",
      "\t</SCOREFXNS>\n",
      "\t<RESIDUE_SELECTORS/>\n",
      "\t<TASKOPERATIONS>\n",
      "\t\t<RestrictToRepacking name=\"no_design\"/>\n",
      "\t\t<ExtraRotamersGeneric ex1=\"1\" ex1_sample_level=\"1\" ex2=\"1\" ex2_sample_level=\"1\" name=\"extrachi\"/>\n",
      "\t</TASKOPERATIONS>\n",
      "\t<FILTERS/>\n",
      "\t<MOVERS>\n",
      "\t\t<PackRotamersMover name=\"pack1\" scorefxn=\"r15\" task_operations=\"no_design,extrachi\"/>\n",
      "\t</MOVERS>\n",
      "\t<APPLY_TO_POSE/>\n",
      "\t<PROTOCOLS>\n",
      "\t\t<Add mover=\"pack1\"/>\n",
      "\t</PROTOCOLS>\n",
      "\t<OUTPUT scorefxn=\"r15\"/>\n",
      "</ROSETTASCRIPTS>\n",
      "\u001b[0mcore.scoring.ScoreFunctionFactory: {0} \u001b[0mSCOREFUNCTION: \u001b[32mref2015\u001b[0m\n",
      "\u001b[0mprotocols.jd2.parser.ScoreFunctionLoader: {0} \u001b[0mdefined score function \"r15\" with weights \"ref2015\"\n",
      "\u001b[0mprotocols.jd2.parser.TaskOperationLoader: {0} \u001b[0mDefined TaskOperation named \"no_design\" of type RestrictToRepacking\n",
      "\u001b[0mprotocols.jd2.parser.TaskOperationLoader: {0} \u001b[0mDefined TaskOperation named \"extrachi\" of type ExtraRotamersGeneric\n",
      "\u001b[0mcore.pack.task.xml_util: {0} \u001b[0mObject pack1 reading the following task_operations: Adding the following task operations\n",
      "no_design extrachi\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mDefined mover named \"pack1\" of type PackRotamersMover\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mParsedProtocol mover with the following settings\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mAdded mover \"pack1\"\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0m=======================BEGIN MOVER PackRotamersMover - pack1=======================\n",
      "\u001b[0mcore.pack.pack_rotamers: {0} \u001b[0mbuilt 2246 rotamers at 76 positions.\n",
      "\u001b[0mcore.pack.pack_rotamers: {0} \u001b[0mRequesting all available threads for interaction graph computation.\n",
      "\u001b[0mcore.pack.interaction_graph.interaction_graph_factory: {0} \u001b[0mInstantiating DensePDInteractionGraph\n",
      "\u001b[0mbasic.thread_manager.RosettaThreadManager: {?} \u001b[0mCreating a thread pool of 1 threads.\n",
      "\u001b[0mbasic.thread_manager.RosettaThreadPool: {?} \u001b[0mLaunched 0 new threads.\n",
      "\u001b[0mcore.pack.rotamer_set.RotamerSets: {0} \u001b[0mCompleted interaction graph pre-calculation in 1 available threads (1 had been requested).\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0msetting status to success\n"
     ]
    }
   ],
   "source": [
    "pose_2 = original_pose.clone()\n",
    "\n",
    "xml_2 = XmlObjects.create_from_string(RS_pack)\n",
    "\n",
    "protocol_2 = xml_2.get_mover(\"ParsedProtocol\")\n",
    "\n",
    "protocol_2.apply(pose_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. RS in Pyrosetta VS python scripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "</br>\n",
    "\n",
    "下面我们以pyrosetta的python脚本和RS分别实现同样的过程,过程中感受一下两种方式的优缺点。\n",
    "\n",
    "这里我们结合以上两个小节:精确定义MoveMap的Minimization(MinMover)和精确定义TaskOperation的Packing(PackRotamersMover)结合起来构成一个新的**RS protocol**,要求:\n",
    "\n",
    "1. 全部过程使用ref2015能量函数;\n",
    "2. 先使用定义好的pack1(与第二小节保持一致)进行packing;\n",
    "3. 再使用定义好的min_torsion(与第一小节保持一致)进行minimization;\n",
    "4. 计算初始pose和结束pose的能量值\n",
    "\n",
    "其次,我们再使用**pyrostta的python脚本版本**进行这样的操作。比较两者。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**使用RS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mGenerating XML Schema for rosetta_scripts...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mInitializing schema validator...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mValidating input script...\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0m...done\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mParsed script:\n",
      "<ROSETTASCRIPTS>\n",
      "\t<SCOREFXNS>\n",
      "\t\t<ScoreFunction name=\"ref_2015\" weights=\"ref2015\"/>\n",
      "\t</SCOREFXNS>\n",
      "\t<RESIDUE_SELECTORS/>\n",
      "\t<TASKOPERATIONS>\n",
      "\t\t<RestrictToRepacking name=\"no_design\"/>\n",
      "\t\t<ExtraRotamersGeneric ex1=\"1\" ex1_sample_level=\"1\" ex2=\"1\" ex2_sample_level=\"1\" name=\"extrachi\"/>\n",
      "\t</TASKOPERATIONS>\n",
      "\t<FILTERS/>\n",
      "\t<MOVERS>\n",
      "\t\t<MinMover bb=\"1\" cartesian=\"F\" chi=\"true\" name=\"min_torsion\" scorefxn=\"ref_2015\">\n",
      "\t\t\t<MoveMap name=\"min_torsion_mm\">\n",
      "\t\t\t\t<Span bb=\"false\" begin=\"1\" chi=\"false\" end=\"999\"/>\n",
      "\t\t\t\t<Span bb=\"true\" begin=\"1\" chi=\"true\" end=\"50\"/>\n",
      "\t\t\t\t<Span bb=\"false\" begin=\"5\" chi=\"true\" end=\"10\"/>\n",
      "\t\t\t</MoveMap>\n",
      "\t\t</MinMover>\n",
      "\t\t<PackRotamersMover name=\"pack1\" scorefxn=\"ref_2015\" task_operations=\"no_design,extrachi\"/>\n",
      "\t</MOVERS>\n",
      "\t<APPLY_TO_POSE/>\n",
      "\t<PROTOCOLS>\n",
      "\t\t<Add mover=\"pack1\"/>\n",
      "\t\t<Add mover=\"min_torsion\"/>\n",
      "\t</PROTOCOLS>\n",
      "\t<OUTPUT scorefxn=\"ref_2015\"/>\n",
      "</ROSETTASCRIPTS>\n",
      "\u001b[0mcore.scoring.ScoreFunctionFactory: {0} \u001b[0mSCOREFUNCTION: \u001b[32mref2015\u001b[0m\n",
      "\u001b[0mprotocols.jd2.parser.ScoreFunctionLoader: {0} \u001b[0mdefined score function \"ref_2015\" with weights \"ref2015\"\n",
      "\u001b[0mprotocols.jd2.parser.TaskOperationLoader: {0} \u001b[0mDefined TaskOperation named \"no_design\" of type RestrictToRepacking\n",
      "\u001b[0mprotocols.jd2.parser.TaskOperationLoader: {0} \u001b[0mDefined TaskOperation named \"extrachi\" of type ExtraRotamersGeneric\n",
      "\u001b[0mprotocols.minimization_packing.MinMover: {0} \u001b[0mOptions chi, bb: 1, 1 omega: 1\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mDefined mover named \"min_torsion\" of type MinMover\n",
      "\u001b[0mcore.pack.task.xml_util: {0} \u001b[0mObject pack1 reading the following task_operations: Adding the following task operations\n",
      "no_design extrachi\n",
      "\u001b[0mprotocols.rosetta_scripts.RosettaScriptsParser: {0} \u001b[0mDefined mover named \"pack1\" of type PackRotamersMover\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mParsedProtocol mover with the following settings\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mAdded mover \"pack1\"\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0mAdded mover \"min_torsion\"\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0m=======================BEGIN MOVER PackRotamersMover - pack1=======================\n",
      "\u001b[0mcore.pack.pack_rotamers: {0} \u001b[0mbuilt 2246 rotamers at 76 positions.\n",
      "\u001b[0mcore.pack.pack_rotamers: {0} \u001b[0mRequesting all available threads for interaction graph computation.\n",
      "\u001b[0mcore.pack.interaction_graph.interaction_graph_factory: {0} \u001b[0mInstantiating DensePDInteractionGraph\n",
      "\u001b[0mcore.pack.rotamer_set.RotamerSets: {0} \u001b[0mCompleted interaction graph pre-calculation in 1 available threads (1 had been requested).\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0m=======================BEGIN MOVER MinMover - min_torsion=======================\n",
      "\u001b[0mcore.select.residue_selector.ResidueSpanSelector: {0} \u001b[0m\u001b[1m[ WARNING ]\u001b[0m Residue span end designation 'Residue 999' is outside the Pose!\n",
      "\u001b[0mprotocols.rosetta_scripts.ParsedProtocol: {0} \u001b[0msetting status to success\n"
     ]
    }
   ],
   "source": [
    "# RS\n",
    "RS_pack_min = \\\n",
    "\"\"\"\n",
    "<ROSETTASCRIPTS>\n",
    "    <SCOREFXNS>    \n",
    "    \n",
    "        <ScoreFunction name=\"ref_2015\" weights=\"ref2015\" />\n",
    "        \n",
    "    </SCOREFXNS>\n",
    "    <RESIDUE_SELECTORS>\n",
    "    </RESIDUE_SELECTORS>\n",
    "    <TASKOPERATIONS>\n",
    "    \n",
    "        <RestrictToRepacking name=\"no_design\" />  \n",
    "        \n",
    "        <ExtraRotamersGeneric name=\"extrachi\" ex1=\"1\" ex2=\"1\" ex1_sample_level=\"1\" ex2_sample_level=\"1\" />\n",
    "        \n",
    "    </TASKOPERATIONS>\n",
    "    <FILTERS>\n",
    "    </FILTERS>\n",
    "    <MOVERS>\n",
    "    \n",
    "        <MinMover name=\"min_torsion\" scorefxn=\"ref_2015\" chi=\"true\" bb=\"1\" cartesian=\"F\" >\n",
    "            <MoveMap name=\"min_torsion_mm\">                         \n",
    "                <Span begin=\"1\" end=\"999\" chi=\"false\" bb=\"false\" /> \n",
    "                <Span begin=\"1\" end=\"50\" chi=\"true\" bb=\"true\" />    \n",
    "                <Span begin=\"5\" end=\"10\" chi=\"true\" bb=\"false\" />   \n",
    "            </MoveMap>                                              \n",
    "        </MinMover>\n",
    "        \n",
    "        <PackRotamersMover name=\"pack1\" scorefxn=\"ref_2015\" task_operations=\"no_design,extrachi\" />\n",
    "        \n",
    "    </MOVERS>\n",
    "    <APPLY_TO_POSE>\n",
    "    </APPLY_TO_POSE>\n",
    "    <PROTOCOLS>\n",
    "        <Add mover=\"pack1\" /> \n",
    "        <Add mover=\"min_torsion\" />\n",
    "    </PROTOCOLS>\n",
    "    <OUTPUT scorefxn=\"ref_2015\" />\n",
    "</ROSETTASCRIPTS>\n",
    "\"\"\"\n",
    "\n",
    "pose_3 = original_pose.clone()\n",
    "\n",
    "xml_3 = XmlObjects.create_from_string(RS_pack_min)\n",
    "\n",
    "protocol_3 = xml_3.get_mover(\"ParsedProtocol\")\n",
    "\n",
    "protocol_3.apply(pose_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32.66393846892796 -204.98731243884194\n"
     ]
    }
   ],
   "source": [
    "scorefxn_RS = xml_3.get_score_function(\"ref_2015\")\n",
    "before = scorefxn_RS.score(original_pose)\n",
    "after = scorefxn_RS.score(pose_3)\n",
    "print(before, after)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**使用Python Script**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: UserWarning: Import of 'rosetta' as a top-level module is deprecated and may be removed in 2018, import via 'pyrosetta.rosetta'.\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "#Python\n",
    "from pyrosetta import *\n",
    "from pyrosetta.rosetta import *\n",
    "\n",
    "#Core Includes\n",
    "from rosetta.core.kinematics import MoveMap\n",
    "from rosetta.core.kinematics import FoldTree\n",
    "from rosetta.core.pack.task import TaskFactory\n",
    "from rosetta.core.pack.task import operation\n",
    "from rosetta.core.simple_metrics import metrics\n",
    "from rosetta.core.select import residue_selector as selections\n",
    "from rosetta.core import select\n",
    "from rosetta.core.select.movemap import *\n",
    "from rosetta.protocols import minimization_packing as pack_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyRosetta-4 2021 [Rosetta PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release 2021.31+release.c7009b3115c22daa9efe2805d9d1ebba08426a54 2021-08-07T10:04:12] retrieved from: http://www.pyrosetta.org\n",
      "(C) Copyright Rosetta Commons Member Institutions. Created in JHU by Sergey Lyskov and PyRosetta Team.\n",
      "\u001b[0mcore.init: {0} \u001b[0mChecking for fconfig files in pwd and ./rosetta/flags\n",
      "\u001b[0mcore.init: {0} \u001b[0mRosetta version: PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release r292 2021.31+release.c7009b3115c c7009b3115c22daa9efe2805d9d1ebba08426a54 http://www.pyrosetta.org 2021-08-07T10:04:12\n",
      "\u001b[0mcore.init: {0} \u001b[0mcommand: PyRosetta -ex1 -ex2aro -database /opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database\n",
      "\u001b[0mbasic.random.init_random_generator: {0} \u001b[0m'RNG device' seed mode, using '/dev/urandom', seed=-208345754 seed_offset=0 real_seed=-208345754 thread_index=0\n",
      "\u001b[0mbasic.random.init_random_generator: {0} \u001b[0mRandomGenerator:init: Normal mode, seed=-208345754 RG_type=mt19937\n"
     ]
    }
   ],
   "source": [
    "init()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PDB file name: ./data/1ubq_clean.pdb\n",
      "Total residues: 76\n",
      "Sequence: MQIFVKTLTGKTITLEVEPSDTIENVKAKIQDKEGIPPDQQRLIFAGKQLEDGRTLSDYNIQKESTLHLVLRLRGG\n",
      "Fold tree:\n",
      "FOLD_TREE  EDGE 1 76 -1 \n"
     ]
    }
   ],
   "source": [
    "pose_4 = original_pose.clone()\n",
    "# 查看pose基本信息\n",
    "print(pose)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use TaskFactory to control TaskOperations -- packing\n",
    "tf = TaskFactory()\n",
    "# InitializeFromCommandline会调用初始init()中的选项\n",
    "tf.push_back(operation.InitializeFromCommandline())\n",
    "# RestrictToRepacking可以限制不允许进行design\n",
    "tf.push_back(operation.RestrictToRepacking())\n",
    "# ExtraRotamersGeneric允许额外rotamer采样\n",
    "extra_rotamer = operation.ExtraRotamersGeneric()\n",
    "extra_rotamer.ex1(True)\n",
    "extra_rotamer.ex2(True)\n",
    "extra_rotamer.ex1_sample_level(pyrosetta.rosetta.core.pack.task.ExtraRotSample.EX_ONE_STDDEV)\n",
    "extra_rotamer.ex2_sample_level(pyrosetta.rosetta.core.pack.task.ExtraRotSample.EX_ONE_STDDEV)\n",
    "tf.push_back(extra_rotamer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "selection_1_4 = selections.ResidueIndexSelector(\"1-4\")\n",
    "selection_5_10 = selections.ResidueIndexSelector(\"5-10\")\n",
    "selection_11_50 = selections.ResidueIndexSelector(\"11-50\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "residue 1 - 4 vector1_bool[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
      "residue 5 - 10 vector1_bool[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
      "residue 11 - 50 vector1_bool[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n"
     ]
    }
   ],
   "source": [
    "print(\"residue 1 - 4\", selection_1_4.apply(pose))\n",
    "print(\"residue 5 - 10\", selection_5_10.apply(pose))\n",
    "print(\"residue 11 - 50\", selection_11_50.apply(pose))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "-------------------------------\n",
      "  resnum     Type  TRUE/FALSE \n",
      "-------------------------------\n",
      " DEFAULT      BB     FALSE\n",
      " DEFAULT      SC     FALSE\n",
      " DEFAULT      NU     FALSE\n",
      " DEFAULT  BRANCH     FALSE\n",
      "     001      BB      TRUE\n",
      "              SC      TRUE\n",
      "     002      BB      TRUE\n",
      "              SC      TRUE\n",
      "     003      BB      TRUE\n",
      "              SC      TRUE\n",
      "     004      BB      TRUE\n",
      "              SC      TRUE\n",
      "     005      SC      TRUE\n",
      "     006      SC      TRUE\n",
      "     007      SC      TRUE\n",
      "     008      SC      TRUE\n",
      "     009      SC      TRUE\n",
      "     010      SC      TRUE\n",
      "     011      BB      TRUE\n",
      "              SC      TRUE\n",
      "     012      BB      TRUE\n",
      "              SC      TRUE\n",
      "     013      BB      TRUE\n",
      "              SC      TRUE\n",
      "     014      BB      TRUE\n",
      "              SC      TRUE\n",
      "     015      BB      TRUE\n",
      "              SC      TRUE\n",
      "     016      BB      TRUE\n",
      "              SC      TRUE\n",
      "     017      BB      TRUE\n",
      "              SC      TRUE\n",
      "     018      BB      TRUE\n",
      "              SC      TRUE\n",
      "     019      BB      TRUE\n",
      "              SC      TRUE\n",
      "     020      BB      TRUE\n",
      "              SC      TRUE\n",
      "     021      BB      TRUE\n",
      "              SC      TRUE\n",
      "     022      BB      TRUE\n",
      "              SC      TRUE\n",
      "     023      BB      TRUE\n",
      "              SC      TRUE\n",
      "     024      BB      TRUE\n",
      "              SC      TRUE\n",
      "     025      BB      TRUE\n",
      "              SC      TRUE\n",
      "     026      BB      TRUE\n",
      "              SC      TRUE\n",
      "     027      BB      TRUE\n",
      "              SC      TRUE\n",
      "     028      BB      TRUE\n",
      "              SC      TRUE\n",
      "     029      BB      TRUE\n",
      "              SC      TRUE\n",
      "     030      BB      TRUE\n",
      "              SC      TRUE\n",
      "     031      BB      TRUE\n",
      "              SC      TRUE\n",
      "     032      BB      TRUE\n",
      "              SC      TRUE\n",
      "     033      BB      TRUE\n",
      "              SC      TRUE\n",
      "     034      BB      TRUE\n",
      "              SC      TRUE\n",
      "     035      BB      TRUE\n",
      "              SC      TRUE\n",
      "     036      BB      TRUE\n",
      "              SC      TRUE\n",
      "     037      BB      TRUE\n",
      "              SC      TRUE\n",
      "     038      BB      TRUE\n",
      "              SC      TRUE\n",
      "     039      BB      TRUE\n",
      "              SC      TRUE\n",
      "     040      BB      TRUE\n",
      "              SC      TRUE\n",
      "     041      BB      TRUE\n",
      "              SC      TRUE\n",
      "     042      BB      TRUE\n",
      "              SC      TRUE\n",
      "     043      BB      TRUE\n",
      "              SC      TRUE\n",
      "     044      BB      TRUE\n",
      "              SC      TRUE\n",
      "     045      BB      TRUE\n",
      "              SC      TRUE\n",
      "     046      BB      TRUE\n",
      "              SC      TRUE\n",
      "     047      BB      TRUE\n",
      "              SC      TRUE\n",
      "     048      BB      TRUE\n",
      "              SC      TRUE\n",
      "     049      BB      TRUE\n",
      "              SC      TRUE\n",
      "     050      BB      TRUE\n",
      "              SC      TRUE\n",
      "-------------------------------\n",
      " jumpnum     Type  TRUE/FALSE \n",
      "-------------------------------\n",
      " DEFAULT     JUMP    FALSE\n",
      "-------------------------------\n",
      "  resnum  atomnum     Type  TRUE/FALSE \n",
      "-------------------------------\n",
      " DEFAULT               PHI    FALSE\n",
      " DEFAULT             THETA    FALSE\n",
      " DEFAULT                 D    FALSE\n",
      " DEFAULT               RB1    FALSE\n",
      " DEFAULT               RB2    FALSE\n",
      " DEFAULT               RB3    FALSE\n",
      " DEFAULT               RB4    FALSE\n",
      " DEFAULT               RB5    FALSE\n",
      " DEFAULT               RB6    FALSE\n",
      "-------------------------------\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Use MoveMapFactory to control MoveMap -- Minimization\n",
    "# 默认为自由度全部关闭状态!按照需求打开自由度\n",
    "mmf = MoveMapFactory()\n",
    "# 开放1 - 4 残基骨架和侧链的自由度\n",
    "mmf.add_chi_action(mm_enable, selection_1_4)\n",
    "mmf.add_bb_action(mm_enable, selection_1_4)\n",
    "# 仅开放5 - 10号允许CHI运动\n",
    "mmf.add_chi_action(mm_enable, selection_5_10)\n",
    "# 开放11 - 50 残基骨架和侧链的自由度\n",
    "mmf.add_chi_action(mm_enable, selection_11_50)\n",
    "mmf.add_bb_action(mm_enable, selection_11_50)\n",
    "# 关闭51 - 76任何自由度:无需操作,默认关闭\n",
    "mm  = mmf.create_movemap_from_pose(pose)\n",
    "# 查看MoveMap,观察是否与设想的一致\n",
    "print(mm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0mcore.pack.task: {0} \u001b[0mPacker task: initialize from command line()\n",
      "\u001b[0mcore.pack.pack_rotamers: {0} \u001b[0mbuilt 2246 rotamers at 76 positions.\n",
      "\u001b[0mcore.pack.pack_rotamers: {0} \u001b[0mRequesting all available threads for interaction graph computation.\n",
      "\u001b[0mcore.pack.interaction_graph.interaction_graph_factory: {0} \u001b[0mInstantiating DensePDInteractionGraph\n",
      "\u001b[0mcore.pack.rotamer_set.RotamerSets: {0} \u001b[0mCompleted interaction graph pre-calculation in 1 available threads (1 had been requested).\n",
      "32.66393846892796 -208.75941216807576\n"
     ]
    }
   ],
   "source": [
    "# run protocol as follows:\n",
    "scorefxn = create_score_function( \"ref2015\" )\n",
    "\n",
    "# 首先进行packing -- 对应于pack1\n",
    "packer = pack_min.PackRotamersMover()\n",
    "packer.score_function(scorefxn)\n",
    "packer.task_factory(tf)\n",
    "packer.apply(pose_4)\n",
    "\n",
    "# 再进行minimization -- 对应于min_torsion\n",
    "minimizer = pack_min.MinMover()\n",
    "minimizer.score_function(scorefxn)\n",
    "minimizer.cartesian(False)\n",
    "minimizer.set_movemap(mm)\n",
    "minimizer.apply(pose_4)\n",
    "\n",
    "# calculate scores\n",
    "before = scorefxn.score(original_pose)\n",
    "after = scorefxn.score(pose_4)\n",
    "print(before, after)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 小结与练习\n",
    "\n",
    "在RS中进行的定义可以与Pyrosetta的python版本一一对应。\n",
    "\n",
    "1. 对比两种方式,各有什么优缺点?你会选择什么方式定制自己的protocol呢?\n",
    "\n",
    "2. 请将./data/Example2-MinMover.xml改为python script,并尝试在两个MinMover中设置不同的MoveMap。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}