<jupyter_start><jupyter_text>Getting started with PyRosetta and PyMOLRosetta is a software suite for working with models of proteins (or other molecules, like small molecule ligands and RNA).We can access all of Rosetta's functions using the `pyrosetta` python package.After navigating to your proteindesign directory, activate your conda environment and start a python terminal:```bashcd /path/to/proteindesignconda activate proteindesignpython```<jupyter_code>from pyrosetta import *
pyrosetta.init()<jupyter_output>PyRosetta-4 2020 [Rosetta PyRosetta4.Release.python38.mac 2020.28+release.8ecab77aa50ac1301efe53641e07e09ac91fee3b 2020-07-07T16:41:06] retrieved from: http://www.pyrosetta.org
(C) Copyright Rosetta Commons Member Institutions. Created in JHU by <NAME> and PyRosetta Team.
[0mcore.init: [0mChecking for fconfig files in pwd and ./rosetta/flags
[0mcore.init: [0mRosetta version: PyRosetta4.Release.python38.mac r260 2020.28+release.8ecab77aa50 8ecab77aa50ac1301efe53641e07e09ac91fee3b http://www.pyrosetta.org 2020-07-07T16:41:06
[0mcore.init: [0mcommand: PyRosetta -ex1 -ex2aro -database /Users/cjmathy/.pyenv/versions/3.8.4/lib/python3.8/site-packages/pyrosetta-2020.28+release.8ecab77aa50-py3.8-macosx-10.15-x86_64.egg/pyrosetta/database
[0mbasic.random.init_random_generator: [0m'RNG device' seed mode, using '/dev/urandom', seed=340126414 seed_offset=0 real_seed=340126414
[0mbasic.random.init_random_generator: [0mRandomGenerator:init: Normal mode, seed=340126414 RG_type=mt1993[...]<jupyter_text>Storing structure information in Rosetta with the Pose classProtein structures are stored in text files that contain the sequence information and the positions of each atom in cartesian space. We'll work with the `PDB` format , but there are other useful ones that Rosetta can use, such as the newer `mmCIF` format which accommates much larger models.After reading in a `PDB` file, Rosetta stores the structure as an object of the `Pose` class. This class contains the atoms in the structure and the chemical connectivity. Additionally, it can store energies for each residue, any constraints we want to apply to the structure, and annotations from the original PDB. To manipulate a structure with Rosetta, we make operations on the `Pose` object.We read in a `PDB` file with the function `pose_from_pdb` in the `pyrosetta` package. Rosetta will print out lots of information about what it is parsing from the `PDB` file, choices it makes to pick a final representation, and compute the score of the structure using a default score function. The score is in Rosetta Energy Units/REU, which you can think of as something like a ΔG free energy.We will work with the PDB "1QYS", which is the structure of Top7, a novel fold of protein designed using Rosetta!You can download it by clicking on https://files.rcsb.org/download/1QYS.pdb. Then move it into your directory for the course. Something like:```mv ~/Downloads/1qys.pdb ~/proteindesign```<jupyter_code>pose = pyrosetta.pose_from_pdb('1qys.pdb')<jupyter_output>[0mcore.chemical.GlobalResidueTypeSet: [0mFinished initializing fa_standard residue type set.  Created 981 residue types
[0mcore.chemical.GlobalResidueTypeSet: [0mTotal time to initialize 0.890463 seconds.
[0mcore.import_pose.import_pose: [0mFile '1qys.pdb' automatically determined to be of type PDB
[0mcore.chemical.GlobalResidueTypeSet: [0mLoading (but possibly not actually using) 'MSE' from the PDB components dictionary for residue type 'pdb_MSE'
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0mReading MSE as MET!
[0mcore.io.pose_from_sfr.PoseFromSFRBuilder: [0[...]<jupyter_text>We can get lots of useful information from this `Pose` object, including:- The structure's sequence- The number of chains- The residue numbering (when residues are missing from a structure, the `Pose` numbering will not match the `PDB` numbering, as the `Pose` numbering always starts at 1.- The atomic position (X,Y,Z) coordinates of any atom, such as the alpha-carbons of two residues (and therefore the distance between them). The units are angstroms.- The energy of individual residues- The energy of the structure<jupyter_code># sequence
pose.sequence()
# chain and numbering info
print(pose.pdb_info())
# CA positions and distance
CA_1 = pose.residue(1).xyz("CA")
CA_2 = pose.residue(2).xyz("CA")

print(CA_1)
print(CA_2)
# Distance between the two residues, using the function norm()
print((CA_1 - CA_2).norm())
# the energy terms for the 5th residue
pose.energies().show(5)
# Structure energy
scorefxn = get_fa_scorefxn()  # a score function is just the equation that relates
scorefxn(pose)<jupyter_output>[0mcore.scoring.ScoreFunctionFactory: [0mSCOREFUNCTION: [32mref2015[0m<jupyter_text>A good resource for seeing everything the `Pose` object has to offer is Figure 19.2 in the following paper describing Rosetta:     - <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., ... & <NAME>. (2011). ROSETTA3: an object-oriented software suite for the simulation and design of macromolecules. Methods in enzymology, 487, 545-574.    https://www.sciencedirect.com/science/article/pii/B9780123812704000196  Visualizing your pose with PyMOLWhile Rosetta doesn't have its own GUI for viewing protein models, it is set up for easy linking to PyMOL.Open a new terminal, activate your environment, and open PyMOL:```bashconda activate proteindesignpymol```In PyMOL, change your directory to the PyRosetta directory, and then run the `PyMOLRosettaServer.py` script.```bashcd /path/to/anaconda/envs/proteindesign/lib/python3.7/site-packages/pyrosettarun PyMOLRosettaServer.py```PyMOL should display a message that says,`PyMOL  PyRosetta link started!`NOTE: You can also make a file `.pymolrc` in your home directory and put the code to run the `PyMOLRosettaServer.py` line of code there for it to run automatically when you initialize PyMOL.     Once the link is set up, we can use a new class, the `PyMOLMover` to send the pose to PyMOL. A `Mover` is a class of Rosetta objects that operate on / alter `Pose` objects.<jupyter_code>pmm = PyMOLMover()
pmm.apply(pose)<jupyter_output><empty_output><jupyter_text>You should now see Top7 in your PyMOL window. Let's color each residue according to an energy term, `fa_sol`. This term accounts for the energetic penalty or favorability of a given side chain's solvent accessibility in the structure.<jupyter_code>pmm.label_energy(pose, 'fa_sol')<jupyter_output><empty_output><jupyter_text>You should see a striped pattern along the beta sheet. This is expected, as residues whose side chains point inwards towards the core are more blue (more negative/energetically favorable), while residues with solvent-exposed side chains are more red, as they pay a solvation penalty.This favorable tertiary structure pattern of residues - a hydrophobic core - was incorporated intentionally to enable folding into a stable protein.If you'd like to examine the way the side chains fit together within the core, you can see them by typing the following command into the PyMOL console: `show sticks, not (name c+n+o) and not hydro`As you can see, PyMOL is a very powerful tool for visualizing our structures over the course of designing and modeling our proteins! The following resources are great for getting used to PyMOL- Tutorial: https://pymolwiki.org/index.php/Practical_Pymol_for_Beginners- Cheat Sheet: https://pymolwiki.org/index.php/CheatSheet- List of all commands: https://pymolwiki.org/index.php/Category:Commands  Relaxing a structureYou may have noticed that when we scored the pose above, the energy was positive (~200 REU), which means Rosetta thinks our pose is pretty unfavorable. Often, the bond lengths and angles of structures from the PDB are not ideal according to Rosetta. This is for a variety of reasons, both physical (crystal contacts introduce a bias, for example) and algorithmic (the refinement forcefields used to solve crystal structures and the score function in Rosetta each have their distinct biases). Typically when manipulating or analyzing a new PDB, we first optimize the structure for the Rosetta score function using a protocol known as `Relax`.Relaxing a structure consists of trying many slight perturbations to the backbone bonds ("minimization") as well as trying other sidechain rotamers ("repacking"). We will learn more about these protocols later, but for now let's just run the `Relax` code and see how it changes our structure.This run should take about 1 minute, and will spit out a lot of log information, as it iteratively. perturbs and re-scores the structure.<jupyter_code># make a copy of the pose
pose_relax = pose.clone()
pose_relax.pdb_info().name('1qys_relaxed')

# create a Mover that can run Relax on our pose
relax_mover = pyrosetta.rosetta.protocols.relax.FastRelax()
relax_mover.set_scorefxn(scorefxn)                   # we use the default score function, which we initialized earlier
relax_mover.constrain_relax_to_start_coords(True)    # we want to stay close to the starting coordinates, just refining a bit

# run relax
relax_mover.apply(pose_relax)
pmm.apply(pose_relax)<jupyter_output><empty_output><jupyter_text>These two poses likely look very different, but that is in fact because the relaxed pose is translated and rotated in the PyMOL coordinate system. To really compare the two structures, run the following command in the PyMOL console `align 1qys.pdb, 1qys_relaxed`. This will align the original structure to the relaxed structure. The two structures should seem much more similar. PyMOL calculates the degree of difference between the two models after alignment, reported as a Root Mean Squared Deviation (RMSD). It should be about 0.7 or 0.8.We can also compute an RMSD metric in Rosetta - the function we will use only computes distances between matched alpha-carbons between the two structures:<jupyter_code>pyrosetta.rosetta.core.scoring.CA_rmsd(pose, pose_relax)<jupyter_output><empty_output><jupyter_text>We can also see if the side chains have moved around. Type the following command into your PyMOL console: `show sticks, not (name c+n+o) and not hydro`. You should see that some of the residues have different rotamers for their side chains, such a Lys31 or Asp38.Finally, let's see if the `Relax` protocol did improve the score (energy) of the structure, as intended:<jupyter_code>initial_score = scorefxn(pose)
relaxed_score = scorefxn(pose_relax)

print(initial_score)
print(relaxed_score)
print(relaxed_score - initial_score)<jupyter_output>206.2147741464653
-304.9643401505464
-511.17911429701167<jupyter_text>Are you surprised by the large difference in score? Despite the overall change in structure being less than an Angstrom (using the RMSD metrics), there was a change of about ~500 REU! This was all from slight optimizations in protein backbone angles and side chain rotamer interactions. The score function is parameterized for 1 REU to be on the order of 0.1-1 kcal/mol, so this is a sizeable difference.  Mutating a poseProtein design involves iteratively changing a protein sequence with new mutations and seeing if the resulting structure scores better in a desired metric (usually energy). Let's try this out by mutating a residue in the relaxed Top7 pose and seeing how this changes the overall energy. We will use a convenient function, `mutate_residue` to mutate a hydrophobic residue (Val 8) in the core of protein to a less favorable polar residue (Thr). We will need to get the pose number for residue 8 first.<jupyter_code>pose_V8T = pose_relax.clone()
pose_V8T.pdb_info().name('1qys_V8T')

res8 = pose_V8T.pdb_info().pdb2pose(chain='A',res=8)
pyrosetta.toolbox.mutate_residue(pack_or_pose=pose_V8T,
                                 mutant_position=res8,
                                 mutant_aa='T',
                                 pack_radius=0)<jupyter_output>[0mcore.scoring.ScoreFunctionFactory: [0mSCOREFUNCTION: [32mref2015[0m
[0mcore.pack.task: [0mPacker task: initialize from command line()
[0mcore.pack.pack_rotamers: [0mbuilt 108 rotamers at 1 positions.
[0mcore.pack.interaction_graph.interaction_graph_factory: [0mInstantiating PDInteractionGraph<jupyter_text>We can confirm that we mutated the residue correctly by querying the name of the residue at that position in both our original and mutated pose<jupyter_code># confirm we mutated correctly
print(pose_relax.residue(res8).name())
print(pose_V8T.residue(res8).name())<jupyter_output>VAL
THR<jupyter_text>Now let's score the two poses<jupyter_code>mutated_score = scorefxn(pose_V8T)

print(mutated_score - relaxed_score)<jupyter_output>4.795649293731344<jupyter_text>We increased the score by ~5 REU, so Rosetta agrees with our intuition that the hydrophobic sidechain Val was preferred over Thr at that core position.We can look at the individual score terms to see which physical force accounts for this change in score<jupyter_code>print('Val8 terms, relaxed structure:\n------------')
pose_relax.energies().show(res8)
print('\nThr8 terms, mutated structure:\n------------')
pose_V8T.energies().show(res8)<jupyter_output>Val8 terms, relaxed structure:
------------
[0mcore.scoring.Energies: [0mE               fa_atr        fa_rep        fa_sol  fa_intra_repfa_intra_sol_x   lk_ball_wtd       fa_elec     pro_close   hbond_sr_bb   hbond_lr_bb   hbond_bb_sc      hbond_sc     dslf_fa13         omega        fa_dun       p_aa_pp yhh_planarity           ref   rama_prepro
[0mcore.scoring.Energies: [0mE(i)   6         -6.97          1.36          1.16          3.17          0.05         -0.12         -1.63          0.00          0.00          0.00          0.00          0.00          0.00          0.29          0.01         -1.13          0.00          2.64         -0.12

Thr8 terms, mutated structure:
------------
[0mcore.scoring.Energies: [0mE               fa_atr        fa_rep        fa_sol  fa_intra_repfa_intra_sol_x   lk_ball_wtd       fa_elec     pro_close   hbond_sr_bb   hbond_lr_bb   hbond_bb_sc      hbond_sc     dslf_fa13         omega        fa_dun       p_aa_pp yhh_planarity           ref   rama_[...]<jupyter_text>We see that there was an increase in the `fa_sol` term, which penalizes polar residues like Threonine being buried in the core:<jupyter_code>fa_sol = pyrosetta.rosetta.core.scoring.ScoreType.fa_sol

V8_fa_sol = pose_relax.energies().residue_total_energies(res8)[fa_sol]
T8_fa_sol = pose_V8T.energies().residue_total_energies(res8)[fa_sol]

print(T8_fa_sol - V8_fa_sol)<jupyter_output>2.0033733586316864
{'fraction_non_alphanumeric': 0.07695396145610278, 'fraction_numerical': 0.03305674518201285, 'mean_word_length': 4.673880030372057, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 52, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '120381', 'license': 'MIT', 'n_tokens_mistral': 5041, 'n_tokens_neox': 4527, 'n_words': 1686, 'path': 'course_website/00_getting_started_with_rosetta.ipynb', 'repo_name': 'Kortemme-Lab/protein_design_course_2021', 'size': '14980'}