I will not lie: I often struggle to find a snippet of code that did something in PyRosetta or I spend hours facing a problem caused by something not working as I expect it to. I recently did a tricky project involving RFdiffusion and I kept slipping on the PyRosetta side. So to make future me, others, and ChatGTP5 happy, here are some common operations to make working with PyRosetta for RFdiffusion easier.
Import
For easy copypasting, let’s start PyRosetta and import things. As said in previous post I am not keen on star imports, because they are not meant to be used indiscriminately and can cause issues —a common pythonic example is when someone mixes collections.Counter
and typing.Counter
.
import pyrosetta import pyrosetta_help as ph from types import ModuleType # Better than star imports: prc: ModuleType = pyrosetta.rosetta.core prp: ModuleType = pyrosetta.rosetta.protocols pru: ModuleType = pyrosetta.rosetta.utility prn: ModuleType = pyrosetta.rosetta.numeric pr_conf: ModuleType = pyrosetta.rosetta.core.conformation pr_scoring: ModuleType = pyrosetta.rosetta.core.scoring pr_res: ModuleType = pyrosetta.rosetta.core.select.residue_selector pr_options: ModuleType = pyrosetta.rosetta.basic.options logger = ph.configure_logger() pyrosetta.distributed.maybe_init(extra_options=ph.make_option_string(no_optH=False, ex1=None, ex2=None, ignore_unrecognized_res=True, load_PDB_components=False, ignore_waters=True, ) )
Load from a PDB file:
pyrosetta.pose_from_file('👾👾👾.pdb')
or from string
pose = pyrosetta.Pose() prc.import_pose.pose_from_pdbstring(pose, 'ATOM 👾👾👾...')
or from a mmCIF file. Note this will not work straight from PyMol because an extra entry at end is needed as discussed previously. I should say working with mmCIF is uncommon, but does have some nice metadata handling advantages.
pose: pyrosetta.Pose = prc.import_pose.pose_from_file('👾👾👾.cif', read_fold_tree=True, type=prc.import_pose.FileType.CIF_file)
There are several other pose importers. Another useful one, for threading especially, is pyrosetta.pose_from_sequence
is rather helpful —just watch out for the eye of Sauron when threading.
Extra complication happens when dealing with residue types that are not in the database. These need to be added to the pose before loading the file.
As always with PyRosetta, if you pass an illegal value, say prc.import_pose.pose_from_file(Exception)
you’ll get a TypeError
that will tell you the accepted arguments of the overloaded function, which sometimes differs from the help(fun
) info.
A new residue type can be added from a params file or from memory. The official way to create a params file is from the Python 2.7 script from the Rosetta download, but that was not good for my purposes so I wrote a new converter, the rdkit_to_params module, which I also made into a web app due to popular demand.
A params file can be passed via the initialisation options (-extra_res_fa
) or via the following:
pru: ModuleType = pyrosetta.rosetta.utility params_filenames: List[str] = ... pose: pyrosetta.Pose = ... params_vector = pru.vector1_string() params_vector.extend([f for f in params_filenames if f]) pyrosetta.generate_nonstandard_residue_set(pose, params_vector)
Having added a residue type does not mean a residue was added to the pose. For that you’d need to do:
prc: ModuleType = pyrosetta.rosetta.core resiset: prc.chemical.ResidueTypeSet = pyrosetta.generate_nonstandard_residue_set(pose, params_vector) new_res = prc.conformation.ResidueFactory.create_residue( resiset.name_map( name ) ) pose.append_residue_by_jump(new_res, pose.num_jump() + 1)
In the rdkit-to-params module, Params.add_residuetype
adds a params but from a string basically.
def add_residuetype(self, pose: pyrosetta.Pose, params_block: str, name3: str, reset:bool=True) \ -> pyrosetta.rosetta.core.chemical.ResidueTypeSet: """ Adds the params to a copy of a residue type set of a pose. If reset is True it will also save it as the default RTS —it keeps other custom residue types. """ rts = pose.conformation().modifiable_residue_type_set_for_conf(prc.chemical.FULL_ATOM_t) buffer = pyrosetta.rosetta.std.stringbuf(params_block) stream = pyrosetta.rosetta.std.istream(buffer) new = prc.chemical.read_topology_file(stream, name3, rts) rts.add_base_residue_type(new) if reset: pose.conformation().reset_residue_type_set_for_conf(rts) return rts
Connections, either polymeric (LOWER
and UPPER
) or otherwise (CONN1
etc.) are a topic that would require its own blog post, so I won’t cover it here.
PDBInfo
PDB information (residue numbering, chain, occupancy and more) is stored in a pyrosetta.rosetta.core.pose.PDBInfo
instance attached to the pose (accessible via .pdb_info()
method, which returns the instance attached and not a copy).
This very handy. Two of its methods are pose2pdb(res=👾)
and pdb2pose(chain=👾, res=👾)
which allows one to jump between the two standards. A nice thing is that it does not segfault when a pose residue index that does not exist is called (RuntimeError
). If the a PDB residue is request that does not exist it will return 0.
A parenthetical warning — one of its methods is .remarks()
, which is a no go. In the PDB format, REMARK is basically a comment that has really strict tools-specific grammar and is a total mess. The .remarks()
does not work in PyRosetta and segfaults even when the file loaded did indeed have REMARK lines. And passing a pyrosetta.rosetta.core.io.Remarks
instance will not work nor will appending to a Remarks
or even making one:
remarks = pyrosetta.rosetta.core.io.Remarks() remark = pyrosetta.rosetta.core.io.RemarkInfo() remark.value = 'Hello world' remarks.append(remark) # ValueError: vector
So were metadata needed to be added to a file, editing the PDB block might be a better bet —editing the mmCIF dictionary data is even betterer. In a pinch appending comment lines to a PDB helps too (but don’t tell anyone).
Certain operations (RemodelMover, grafting etc.) will make the PDB information obsolete or make it get lost.
The Pose class has a method split_by_chain
which returns the chains (as defined in the fold tree) preserving the PDBInfo, but the function (external to pyrosetta.Pose
) pyrosetta.rosetta.core.pose.append_pose_to_pose
does not, so it needs correcting:
def add_chain(built: pyrosetta.Pose, new: pyrosetta.Pose, reset: bool = False) -> None: """ Add a chain ``new`` to a pose ``built`` preserving the residue numbering. :param built: this is the pyrosetta.Pose that will be built into... :param new: the addendum :param reset: resets the PDBInfo for the chain present to A :return: """ built_pi = built.pdb_info() if built_pi is None or reset: built_pi = prc.pose.PDBInfo(built) built.pdb_info(built_pi) for r in range(1, built.total_residue() + 1): built_pi.set_resinfo(res=r, chain_id='A', pdb_res=r) for chain in new.split_by_chain(): offset: int = built.total_residue() pyrosetta.rosetta.core.pose.append_pose_to_pose(built, chain, new_chain=True) # the new `built` residues will not have PDBinfo chain_pi = chain.pdb_info() for r in range(1, chain.total_residue() + 1): built_pi.set_resinfo(res=r + offset, chain_id=chain_pi.chain(r), pdb_res=chain_pi.number(r)) built_pi.obsolete(False)
The functions pyrosetta.rosetta.protocols.grafting.delete_region
or pyrosetta.rosetta.protocols.grafting.return_region
do not impact the PDBinfo.
Alignment and superposition
First, a silly pedantic lexical side note. When one has a sequence of elements and one spaces out the elements so they match, that is called to align
. When one places one object over another but letting one still be visible, say with a stagger like the traces in a joy plot or translating without rotating (e.g. ▽+△=✡︎) that is called to superimpose
, when one rototranslates fully (e.g. ▽+△↺=▽), that is called to superpose
. Here I will align sequences, to superposed structures.
In PyRosetta one can superpose two poses with the following:
mobile: pyrosetta.Pose = ... ref: pyrosetta.Pose = ... atom_map: prs.map_core_id_AtomID_core_id_AtomID = ... # this is a hash-mapping of mobile `pyrosetta.AtomID` to reference `pyrosetta.AtomID` rmsd: float = pr_scoring.superimpose_pose(mod_pose=mobile, ref_pose=ref, atom_map=atom_map)
There are actually several functions that allow superpositions, including presets for CA atoms or all atoms and so forth. The above is customisable as you have to tell it what atom goes to which.
So where we to align two poses by the CA atoms of a specified PDB chain that ought to be common, we can do:
prc: ModuleType = pyrosetta.rosetta.core prp: ModuleType = pyrosetta.rosetta.protocols pru: ModuleType = pyrosetta.rosetta.utility prn: ModuleType = pyrosetta.rosetta.numeric prs: ModuleType = pyrosetta.rosetta.std pr_conf: ModuleType = pyrosetta.rosetta.core.conformation pr_scoring: ModuleType = pyrosetta.rosetta.core.scoring pr_res: ModuleType = pyrosetta.rosetta.core.select.residue_selector pr_options: ModuleType = pyrosetta.rosetta.basic.options def superpose_pose_by_chain(pose, ref, chain: str, strict: bool=True) -> float: """ superpose by PDB chain letter :param pose: :param ref: :param chain: :return: """ atom_map = prs.map_core_id_AtomID_core_id_AtomID() chain_sele: pr_res.ResidueSelector = pr_res.ChainSelector(chain) r: int # reference pose residue number (Fortran) m: int # mobile pose residue number (Fortran) for r, m in zip(pr_res.selection_positions(chain_sele.apply(ref)), pr_res.selection_positions(chain_sele.apply(pose)) ): if strict: assert pose.residue(m).name3() == ref.residue(r).name3(), 'Mismatching residue positions!' ref_atom = pyrosetta.AtomID(ref.residue(r).atom_index("CA"), r) mobile_atom = pyrosetta.AtomID(pose.residue(m).atom_index("CA"), m) atom_map[mobile_atom] = ref_atom return pr_scoring.superimpose_pose(mod_pose=pose, ref_pose=ref, atom_map=atom_map)
If pyrosetta.AtomID
seems familiar, that is because for constraints one has to play around with them a lot.
Sequence alignment comes into play when superposing protein by the conserved regions. For that Bio.Align.PairwiseAligner
can be used (since pairwise2
was deprecated).
In the default settings, the character ‘-‘ in the input sequence is treated like a regular character and not a free gap. Were one hellbent to do so a custom substitution matrix would be needed.
Were one to want align to homologues, this would be fine:
from Bio.Align import PairwiseAligner, Alignment, substitution_matrices def superpose_by_alignment(mobile: pyrosetta.Pose, ref: pyrosetta.Pose) -> float: """ Pairwise alignment of the sequences of the poses. return (ref_index, mobile_index) :param mobile: :param ref: :return: """ # ## align aligner = PairwiseAligner() aligner.substitution_matrix = substitution_matrices.load("BLOSUM62") ref_seq: str = ref.sequence() pose_seq: str = mobile.sequence() aln: Alignment = aligner.align(ref_seq, pose_seq)[0] # `aln.indices` has the mapping # an index of -1 is a map to a gap aln_map: Dict[int, int] = {t: q for t, q in zip(aln.indices[0], aln.indices[1]) if q != -1 and t != -1} # ## make pyrosetta atom map # for the purpose of explanation the following block is not in its own function, but will be repeated in a minute. # where these all part of the same code, the repeated part would need to be its own function as it is very modular anyway! atom_map = prs.map_core_id_AtomID_core_id_AtomID() for r, m in aln_map.items(): ref_atom = pyrosetta.AtomID(ref.residue(r + offset).atom_index("CA"), r + offset) mobile_atom = pyrosetta.AtomID(mobile.residue(m + offset).atom_index("CA"), m + offset) atom_map[mobile_atom] = ref_atom # ## superpose and return RMSD return pr_scoring.superimpose_pose(mod_pose=mobile, ref_pose=ref, atom_map=atom_map)
However, if the poses to superpose differed by a span that any identity is coincidental, such as a redesigned part of a domain, then we need to strongly penalise mismatches and not penalise gap extensions.
from Bio.Align import PairwiseAligner, Alignment, substitution_matrices def superpose_by_common_alignment(mobile: pyrosetta.Pose, ref: pyrosetta.Pose) -> float: """ Pairwise alignment of the sequences of the poses. This time only spans that are common. return (ref_index, mobile_index) :param mobile: :param ref: :return: """ # ## align aligner = PairwiseAligner() # We don't want to get mismatches aligned, so no to: # aligner.substitution_matrix = substitution_matrices.load("BLOSUM62") aligner.internal_gap_score = -10 aligner.extend_gap_score = -0.01 aligner.end_gap_score = -0.01 ref_seq: str = ref.sequence() pose_seq: str = mobile.sequence() aln: Alignment = aligner.align(ref_seq, pose_seq)[0] # like before but with an extra condition `ref_seq[t] == pose_seq[q]` aln_map: Dict[int, int] = {t: q for t, q in zip(aln.indices[0], aln.indices[1]) if q != -1 and t != -1 and ref_seq[t] == pose_seq[q]} # ## make pyrosetta atom map ... # return RMSD return pr_scoring.superimpose_pose(mod_pose=mobile, ref_pose=ref, atom_map=atom_map)
In the case that one end does not matter in the sequence, there are attributes that control the left (N-terminal) an right (C-terminal) side gaps. For example, the attributes target_right_gap_score
and target_right_extend_gap_score
could be set to zero, which make the C-terminal difference an inkshed of gaps.
So what operations move the pose? A lot. The main two to keep an eye out for are RFdiffusion itself and Relax and its pesky drift.
Thread
In PyRosetta one can mutate a residue via
prp: ModuleType = pyrosetta.rosetta.protocols prp.simple_moves.MutateResidue(1, 'ALA').apply(pose)
This allows mutations to not only base residues, but also patched residues, such as NtermProteinFull
(with extra N-terminal proton) or acetylated
etc. To see if in the database folder there is what you are after it is honestly easier simply looking in the folder:
import pyrosetta, itertools from pathlib import Path from itertools import chain fa = 'database/chemical/residue_type_sets/fa_standard' # residue types print([p.name for p in (Path(pyrosetta.__file__).parent / fa / 'residue_types' ).glob('*/*.params')]) # params types patches_folder = (Path(pyrosetta.__file__).parent / fa / 'patches' ) print([p.name for p in itertools.chain(patches_folder.glob('*.txt'), patches_folder.glob('*/*.txt'))])
Just be aware, that MutateResidue
segfaults on error.
To remove the terminus patches one can do the following:
import pyrosetta prc: ModuleType = pyrosetta.rosetta.core pru: ModuleType = pyrosetta.rosetta.utility def remove_terminus_patches(pose): clean_template_pose = template_pose.clone() prc.pose.remove_nonprotein_residues(pose) ### find lowers = pru.vector1_std_pair_unsigned_long_protocols_sic_dock_Vec3_t() uppers = pru.vector1_std_pair_unsigned_long_protocols_sic_dock_Vec3_t() prc.sic_dock.get_termini_from_pose(pose, lowers, uppers) ### remove for upper, _ in uppers: prc.conformation.remove_upper_terminus_type_from_conformation_residue(clean_template_pose.conformation(), upper) for lower, _ in lowers: prc.conformation.remove_lower_terminus_type_from_conformation_residue(clean_template_pose.conformation(), lower)
However, mutating a whole pose based on a sequence is a pillar of homology modelling, namely threading. I have discussed in an other post threading and hybridisation in PyRosetta and most examples here are simply ported from a helper module of mine, so I’ll be brief. Threading in Rosetta requires a “Grishin file”, which is a specific form of pairwise sequence alignment file. The mover is ThreadingMover
with the aid of StealSideChainsMover
, namely:
prc: ModuleType = pyrosetta.rosetta.core prp: ModuleType = pyrosetta.rosetta.protocols clean_template_pose = remove_terminus_patches(pose) # no termini from earlier! target_pose = pyrosetta.Pose() prc.pose.make_pose_from_sequence(target_pose, target_sequence, 'fa_standard') ## Thread align = prc.sequence.read_aln(format='grishin', filename=aln_file)[1] threader = prp.comparative_modeling.ThreadingMover(align=align, template_pose=clean_template_pose) threader.apply(target_pose) ## Steal sidechains qt = threader.get_qt_mapping(target_pose) steal = prp.comparative_modeling.StealSideChainsMover(clean_template_pose, qt) steal.apply(target_pose)
The main point of the aforementioned post about threading was loop modelling by cannibalising AF2, which can be done with the following settings before using the ThreadingMover
mover:
# optional set for loop modelling from reference poses poses: List[pyrosetta.Pose] = ... # poses to cannibilise for loops lengths: List[int] = [3,] # 3,6,9 are traditional choices fragment_sets = pru.vector1_std_shared_ptr_core_fragment_FragSet_t(len(lengths)) for i, l in enumerate(lengths): fragment_sets[i+1] = prc.fragment.ConstantLengthFragSet(l) for pose in poses: prc.fragment.steal_constant_length_frag_set_from_pose(pose, fragment_sets[i+1]) threader = ... threader.build_loops(True) threader.randomize_loop_coords(True) # default threader.frag_libs(fragment_sets) threader.apply(...)
As may be apparent, the ThreadingMover
does not like non-base residues. So for ligand residues, one can just use append_subpose_to_pose
function.
prc: ModuleType = pyrosetta.rosetta.core pr_res: ModuleType = pyrosetta.rosetta.core.select.residue_selector def steal_ligands(donor_pose, acceptor_pose) -> None: """ Steals non-Protein residues from donor_pose and adds them to acceptor_pose Do not use with nucleic acid polymers. :param donor_pose: :param acceptor_pose: :return: """ PROTEIN = prc.chemical.ResidueProperty.PROTEIN prot_sele = pr_res.ResiduePropertySelector(PROTEIN) not_sele = pr_res.NotResidueSelector(prot_sele) rv = pr_res.ResidueVector(not_sele.apply(donor_pose)) # if it were DNA... # for from_res, to_res in ph.rangify(rv): # prc.pose.append_subpose_to_pose(acceptor_pose, donor_pose, from_res, to_res, True) for res in rv: prc.pose.append_subpose_to_pose(acceptor_pose, donor_pose, res, res, True)
Now back to threading. As mentioned, I just use my helper function to deal with the above:
import pyrosetta_help def thread(template_block, target_seq, target_name, template_name, temp_folder='/tmp'): # load template template = pyrosetta.Pose() prc.import_pose.pose_from_pdbstring(template, template_block) # thread aln_filename = f'{temp_folder}/{template_name}-{target_name}.grishin' ph.write_grishin(target_name=target_name, target_sequence=target_seq, template_name=template_name, template_sequence=template.sequence(), outfile=aln_filename ) aln: prc.sequence.SequenceAlignment = prc.sequence.read_aln(format='grishin', filename=aln_filename)[1] threaded: pyrosetta.Pose threader: prp.comparative_modeling.ThreadingMover threadites: pru.vector1_bool threaded, threader, threadites = ph.thread(target_sequence=target_seq, template_pose=template, target_name=target_name, template_name=template_name, align=aln ) # no need to superpose. It is already aligned # superpose(template, threaded) # fix pdb info n = threaded.total_residue() pi = prc.pose.PDBInfo(n) for i in range(1, n + 1): pi.number(i, i) pi.chain(i, 'A') threaded.pdb_info(pi) return threaded
A pythonic parenthesis. Names of callable objects are, by good practice, verbs, however stuff gets confusing when a word can be both a noun and a verb: I have tripped up so often with thread and fragment. Hence the use of the past particle in the above. For fragment, which has initial-stress derivation I naughtily use fràgment
and fragm
. è
nt
Relax
The workhorse of Rosetta is without question the FastRelax mover. It is very customisable, so most often there is no near to play with packers, Monte Carlo samplers etc. This is true for design mode. I won’t go into detail on the basics as this is covered in every tutorial.
scorefxn: pr_scoring.ScoreFunction = pyrosetta.get_fa_scorefxn() # don't forget to set weights! scorefxn.set_weight(pr_scoring.ScoreType.atom_pair_constraint, 3) scorefxn.set_weight(pr_scoring.ScoreType.coordinate_constraint, 5) cycles = 3 # more than 5 is probably not going to do much relax = prp.relax.FastRelax(scorefxn, cycles) # set movemap to freeze residues # say repack the sidechains of chain A and its neighbours # allow backbone movements for chain A only movemap = pyrosetta.MoveMap() rs: ModuleType = prc.select.residue_selector chainA_sele: rs.ResidueSelector = rs.ChainSelector('A') chainA: pru.vector1_bool = chainA_sele.apply(pose) neigh_sele: rs.ResidueSelector = rs.NeighborhoodResidueSelector(chainA_sele, True, distance) neighs: pru.vector1_bool = neigh_sele.apply(pose) movemap.set_chi(neighs) movemap.set_bb(chainA) movemap.set_jump(False) relax.set_movemap(movemap) # go! replax.apply(pose)
In the case of RFdiffusion the template ought to have been minimised otherwise either the scores will be for which case was minimised better if the other chains are also minimised (don’t) or the scores will be for cases where the designed binder interacts or not with a spuriously strained part of the binding partner.
In the above example the neighbourhood was moved as an example. If you do repack these sidechains, which is reasonable, do make sure to keep the complex and the original binding partner will no longer work.
A cool thing once can do with PyRosetta is use electron design as a constraint (tutorial) —theoretically one could used RFdiffusion via blender API, gemmi and numpy manupulations to make custom shaped protein, but that would be for fun and not actually relevant here.
If there are covalent bond lengths that need fixing, the we have do it cartesian:
constraint_weight = 5 scorefxn: pr_scoring.ScoreFunction = pyrosetta.create_score_function('ref2015_cart') scorefxn.set_weight(pr_scoring.ScoreType.coordinate_constraint, constraint_weight) scorefxn.set_weight(pr_scoring.ScoreType.angle_constraint, constraint_weight) scorefxn.set_weight(pr_scoring.ScoreType.atom_pair_constraint, constraint_weight) relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, 3) relax.cartesian(True) relax.minimize_bond_angles(True) relax.minimize_bond_lengths(True) relax.apply(pose)
In both example I allude to constraints. These are well covered and are set up in various way for example, given the pose index and atom name of the atoms we can do:
def constrain_distance(pose, fore_idx, fore_name, aft_idx, aft_name, x0_in=1.334, sd_in=0.2, tol_in=0.02, weight=1): AtomPairConstraint = pr_scoring.constraints.AtomPairConstraint # noqa fore = pyrosetta.AtomID(atomno_in=pose.residue(fore_idx).atom_index(fore_name), rsd_in=fore_idx) aft = pyrosetta.AtomID(atomno_in=pose.residue(aft_idx).atom_index(aft_name), rsd_in=aft_idx) fun = pr_scoring.func.FlatHarmonicFunc(x0_in=x0_in, sd_in=sd_in, tol_in=tol_in) if weight != 1: fun = pr_scoring.func.ScalarWeightedFunc(weight, fun) con = AtomPairConstraint(fore, aft, fun) pose.add_constraint(con) return con def constrain_angle(pose, fore_idx, fore_name, mid_name, mid_idx, aft_idx, aft_name, x0_in=109/180*3.14, sd_in=10/180*3.14, weight=1): AngleConstraint = pr_scoring.constraints.AngleConstraint fore = pyrosetta.AtomID(atomno_in=pose.residue(fore_idx).atom_index(fore_name), rsd_in=fore_idx) mid = pyrosetta.AtomID(atomno_in=pose.residue(mid_idx).atom_index(mid_name), rsd_in=mid_idx) aft = pyrosetta.AtomID(atomno_in=pose.residue(aft_idx).atom_index(aft_name), rsd_in=aft_idx) fun = pr_scoring.func.CircularHarmonicFunc(x0_radians=x0_in, sd_radians=sd_in) if weight != 1: fun = pr_scoring.func.ScalarWeightedFunc(weight, fun) con = AngleConstraint(fore, mid, aft, fun) pose.add_constraint(con) return con def constrain_position(pose: pyrosetta.Pose, target_index: int, ref_index: int, x0_in=0., sd_in=0.01): ref_ca = pyrosetta.AtomID(atomno_in=pose.residue(ref_index).atom_index('CA'), rsd_in=ref_index) target_ca = pyrosetta.AtomID(atomno_in=pose.residue(target_index).atom_index('CA'), rsd_in=frozen_index) target_xyz = pose.residue(frozen_index).xyz(target_ca.atomno()) fun = pr_scoring.func.HarmonicFunc(x0_in=x0_in, sd_in=sd_in) if weight != 1: fun = pr_scoring.func.ScalarWeightedFunc(weight, fun) con = pr_scoring.constraints.CoordinateConstraint(a1=target_ca, fixed_atom_in=ref_ca, xyz_target_in=target_xyz, func=fun, scotype=pr_scoring.ScoreType.coordinate_constraint) pose.add_constraint(con) def constrain_peptide_gap(pose, chain_break, x0_in=1.334, sd_in=0.2, tol_in=0.02): """ Close a gap between """ AtomPairConstraint = pr_scoring.constraints.AtomPairConstraint # noqa fore_c = pyrosetta.AtomID(atomno_in=pose.residue(chain_break).atom_index('C'), rsd_in=chain_break) aft_n = pyrosetta.AtomID(atomno_in=pose.residue(chain_break + 1).atom_index('N'), rsd_in=chain_break + 1) fun = pr_scoring.func.FlatHarmonicFunc(x0_in=x0_in, sd_in=sd_in, tol_in=tol_in) con = AtomPairConstraint(fore_c, aft_n, fun) pose.add_constraint(con)
One thing to keep an eye out for is how healthy are the constraints after minimisation:
def show_cons(pose): """ print the score for the various constraints """ get_atomname = lambda atomid: pose.residue(atomid.rsd()).atom_name(atomid.atomno()).strip() get_description = lambda atomid: f'{pose.residue(atomid.rsd()).name3()}{pi.pose2pdb(atomid.rsd()).strip().replace(" ", ":")}.{get_atomname(atomid)}' for con in pose.constraint_set().get_all_constraints(): if con.type() == 'AtomPair': print(con.type(), get_description(con.atom1()), get_description(con.atom2()), con.score(pose)) elif con.type() == 'Angle': print(con.type(), get_description(con.atom1()), get_description(con.atom2()), get_description(con.atom3()), con.score(pose)) else: print(con.type(), con.score(pose))
Design
One can do design with FastRelax. Due to the way the mover samples, it fixes sub-optimal residues one by one as opposed to large scale remodelling —for that Remodel might be better suitable (tutorial). Design is really useful for RFdiffusion as the latter can in essence be thought as coarse-grain so there may be some residue combinations that are subpar. Obviously, when FastDesign choses a residue it is based on the scorefunction, which is imperfect: a quirky example of this is its proclivity to change certain residues in Streptag when bound to streptavidin even though the former is the product of rounds and rounds of phage display.
To run design one needs to create a task factory for it onto which restrictions (operations) are pushed. Confusingly, the default is to design everything. The ProhibitSpecifiedBaseResidueTypes
is really handy for preventing cysteine from making life hard and also for correcting surface charges —protein are not very soluble if the isoelectric point is close to pH 7.4 or aggregate if there are exposed patches. On the latter note, proteinMPNN when using the vanilla weights potentially might make inadvertently a transmembrane helix, so that is something to keep an eye out for.
def create_design_tf(pose:pyrosetta.Pose, design_sele: pr_res.ResidueSelector, distance:int) -> prc.pack.task.TaskFactory: """ Create design task factory for relax. Designs the ``design_sele`` and repacks around ``distance`` of it. Remember to do ... code-block:: python relax.set_enable_design(True) relax.set_task_factory(task_factory) """ # design is default, so this is not done: # residues_to_design = design_sele.apply(pose) # design_ops = prc.pack.task.operation.OperateOnResidueSubset(????, residues_to_design) no_cys = pru.vector1_std_string(1) no_cys[1] = 'CYS' no_cys_ops = prc.pack.task.operation.ProhibitSpecifiedBaseResidueTypes(no_cys) # No design, but repack repack_sele = pr_res.NeighborhoodResidueSelector(design_sele, distance, False) residues_to_repack = repack_sele.apply(pose) repack_rtl = prc.pack.task.operation.RestrictToRepackingRLT() repack_ops = prc.pack.task.operation.OperateOnResidueSubset(repack_rtl, residues_to_repack) # No repack, no design frozen_sele = pr_res.NotResidueSelector(pr_res.OrResidueSelector(design_sele, repack_sele)) residues_to_freeze = frozen_sele.apply(pose) prevent_rtl = prc.pack.task.operation.PreventRepackingRLT() frozen_ops = prc.pack.task.operation.OperateOnResidueSubset(prevent_rtl, residues_to_freeze) # pyrosetta.rosetta.core.pack.task.operation.RestrictAbsentCanonicalAASRLT # pyrosetta.rosetta.core.pack.task.operation.PreserveCBetaRLT task_factory = prc.pack.task.TaskFactory() task_factory.push_back(no_cys_ops) task_factory.push_back(repack_ops) task_factory.push_back(frozen_ops) return task_factory
When I used RFdiffusion I gave the algorithm a range in the contigmap, thus allowing the best length of an insertion to win. This does however make the pipeline a bit complicated. First I had to adapt the helper scripts of ProteinMPNN to be used as Python functions (see footnote). There are then two options for the design step: refer to the original skeleton sequence, where stretches of glycines mark the designed residues, or do a sequence alignment. The former requires being organised, while the latter could use the alignment code discussed above.
def design_different(pose: pyrosetta.Pose, ref: pyrosetta.Pose, cycles = 5, scorefxn=None): ref = ref.split_by_chain(1) ref2pose: dict = align_for_atom_map(pose.split_by_chain(1), ref) conserved = list(ref2pose.values()) idx_sele = pr_res.ResidueIndexSelector() for i in range(1, len(pose.chain_sequence(1))): if i not in conserved: idx_sele.append_index(i) print(idx_sele.apply(pose)) task_factory: prc.pack.task.TaskFactory = create_design_tf(pose, design_sele=idx_sele, distance=0) if scorefxn is None: scorefxn = pyrosetta.get_fa_scorefxn() relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, cycles) relax.set_enable_design(True) relax.set_task_factory(task_factory) relax.apply(pose)
The degree of change during a FastDesign run depends on how strained the residues are. So it important that the design be FastRelax properly first and not all at once.
Actually I often use FastRelax in a first pass with the backbone fixed to prevent the model from blowing up.
This degree of change _I think_ makes a good metric for scoring. If the pose.sequence() changed by more than 20%, then that is hard no.
Scoring
RFdiffusion will be made to generate lots of designs. The reason is because they are not all amazing. So they need to be scored.
Coarse-grain clashes
First, we need to determine if the polyglycine skeleton will clash with the extended neighbourhood when the protein is part of a greater machinery that could not be feasibly fed into RFdiffusion. An example is a membrane, say one were designing a binder to a transmembrane protein, the design should not cross the membrane, so one could get the transmembrane protein from OPM database, superpose the skeleton and use the layers of O / N atoms. For a project of mine, I had a multidimensional polymeric lattice. Scoring via a scorefunction and using only the Lenard–Jones repulsion is a tad overkill and prone to errors. So simply doing distance based calculations works:
skeleton: pyrosetta.Pose = ... clasher: pyrosetta.Pose = ... # where the skeleton combined with this pose, it should not have clashes xyz_skeleton = extract_coords(skeleton) xyz_clasher = extract_coords(clasher) all_distances = np.sqrt(np.sum((xyz_clasher[:, np.newaxis, :] - xyz_skeleton[np.newaxis, :, :]) ** 2, axis=-1)) n_clashing = np.count_nonzero(all_distances < 1.)
If designing interactions, the above but with < 3. will given backbone hydrogen bonding to any atom of the reference.
Say the skeleton passed, we then generate the sequence with proteinMPNN, thread it, relax it and tune the sequence. Now we can do some proper scoring.
∆∆G
The first is the predicted Gibbs free energy from a scorefunction. The catch is that this absolute energy, which should not be taken too seriously. The reason being is that they are optimised for ∆∆G calculations, so have empirical terms (i.e. fudge factors) to make them work, including the REF value (a per residue type weight, ref2015 values shown below), which means that were one to make a linear string of peptide it will not have a predicted absolute energy of folding of zero. However, in terms of magnitude is there is a screaming high value then the model is for the bin.
{'ALA': 1.32468, 'ARG': -0.09474, 'ASN': -1.34026, 'ASP': -2.14574, 'CYS': 3.25479, 'GLN': -1.45095, 'GLU': -2.72453, 'GLY': 0.79816, 'HIS': -0.30065, 'ILE': 2.30374, 'LEU': 1.66147, 'LYS': -0.71458, 'MET': 1.65735, 'PHE': 1.21829, 'PRO': -1.64321, 'SER': -0.28969, 'THR': 1.15175, 'TRP': 2.26099, 'TYR': 0.58223, 'VAL': 2.64269}
Per residue score
The scorefunction is a sum of the scores of each particle/residue, and the per-residue scores are obtainable. These are a common way to determine how accurate the model likely is. The per-residue score should not >+10 kcal/mol, and most residues in minimised PDB protein are under 5 kcal/mol, with one or two at the +8 kcal/mol mark likely due to removed ligands and crystalline waters (implicit solvent is an approximation of bulk solvent, not frozen waters). Two things are worth remembering, the ScoreFunction
does not store any data, but Pose.energies
does. But both can access this data. Personally I prefer the latter, but the former has weights which is handy.
prc: ModuleType = pyrosetta.rosetta.core pr_scoring: ModuleType = pyrosetta.rosetta.core.scoring scorefxn = pyrosetta.ScoreFunction = pyrosetta.get_fa_scorefxn() # showcase weights: weights: Dict[str, float] = {w.name: scorefxn.get_weight(w) for w in scorefxn.get_nonzero_weighted_scoretypes()} print(weights) # get per residue scores: res_scores = [] for i in range(1, monomer.total_residue() + 1): v = pru.vector1_bool(pose.total_residue()) v[i] = 1 res_scores.append(scorefxn.get_sub_score(pose, v)) # some magic can be done with the annoying EMapVector v: pr_scoring.EMapVector = scorefxn.weights() print( v[pr_scoring.total_score], v[pr_scoring.fa_sol] )
Where using pose.energies
is a lot easier
scorefxn(pose) # fills pose.energies energies: pr_scoring.Energies = pose.energies() typed_energies: Dict[str, float] = energies.active_total_energies() a: npt.NDArray = energies.residue_total_energies_array() energy = pd.DataFrame(a) # amazingly civilised # alternative a specific residue can be queried via EMapVector as above v: pr_scoring.EMapVector = energies.residue_total_energies(1)
Do note that due to the fact that it takes two particles for a two body score (e.g. hydrogen bond), there’s a function to toggle whether these are halved or not (cf. this post).
scorefxn.weights() # Create the EnergyMap emopts = pr_scoring.methods.EnergyMethodOptions(scorefxn.energy_method_options()) emopts.hbond_options().decompose_bb_hb_into_pair_energies(True) scorefxn.set_energy_method_options(emopts)
Parenthetically, if the pd.DataFrame(pose.energies().residue_total_energies_array())
was not saved, but is needed again, here’s a function to read the junk on a scored dumped PDB.
def rosetta_pdb_to_df(pdbblock: str) -> pd.DataFrame: parsable = False _data = [] for line in pdbblock.split('\n'): if '#BEGIN_POSE_ENERGIES_TABLE' in line: parsable = True continue elif '#END_POSE_ENERGIES_TABLE' in line: break elif not parsable: continue parts = line.strip().split() if parts[0] == 'label': _data.append(parts) elif parts[0] == 'weights': _data.append([parts[0]] + list(map(float, parts[1:-1])) + [float('nan')]) else: _data.append([parts[0]] + list(map(float, parts[1:]))) data = pd.DataFrame(_data) data.columns = data.iloc[0] data = data.iloc[1:].copy() return data
Interface energy
The thing we want to have is a nice interface. There are two ways to this. One is civilised, but prone to segfaults: the InterfaceAnalyzerMover
mover.
def score_interface(complex: Union[pyrosetta.Pose, Sequence[pyrosetta.Pose]], interface: str): if isinstance(complex, Sequence): _complex = complex[0].clone() for c in complex[1:]: add_chain(_complex, c) complex = _complex ia = pyrosetta.rosetta.protocols.analysis.InterfaceAnalyzerMover(interface) ia.apply(complex) return {'complex_energy': ia.get_complex_energy(), 'separated_interface_energy': ia.get_separated_interface_energy(), 'complexed_sasa': ia.get_complexed_sasa(), 'crossterm_interface_energy': ia.get_crossterm_interface_energy(), 'interface_dG': ia.get_interface_dG(), 'interface_delta_sasa': ia.get_interface_delta_sasa()}
Like many other movers, no ligands are allowed.
Interface is a chain letter form, say A_BC
. This crashes out above 3 chains I believe.
The other way is scoring the monomers in isolation (pose.split_chain
) or translating a silly distance after having reset the pose.energies
but preferably repacked the surface layer.
Changed sequence
As said above, how many residues changed when tweaked can be a proxy for something being not quite right.
Number of neighbours
pose2pdb = oligomer.pdb_info().pose2pdb chainA_sele = pr_res.ChainSelector('A') # boolean vector: v = pr_res.NeighborhoodResidueSelector(chainA_sele, distance, False).apply(oligomer) # vector of pose idxs: close_residues = prc.select.get_residues_from_subset(v) # array of PDB numbers print( [pose2pdb(r) for r in close_residues] )
In another blog post I go through the two different residue neighbourhood selectors, but briefly, NeighborhoodResidueSelector
is centroid (roughly the beta carbon), while CloseContactResidueSelector
is closest atom.
cc_sele = pr_res.CloseContactResidueSelector() cc_sele.central_residue_group_selector(chainA_sele) cc_sele.threshold(3) # Å
Hydrophobicity & co.
In a tool developed in OPIG, Therapeutic Antibody Profiler, an set of metrics where characterised, such as patches of surface hydrophobicity etc. which most likely would be well worth using in protein binder design. (I’ve not used these yet so don’t have a snippet to share).
Movement
One may be after rigidity. An MD run is best for that, but pyrosetta can do (discussed in detail in another blog post) via the BackrubMover
mover.
def movement(original: pyrosetta.Pose, trials: int = 100, temperature: int = 1.0, replicate_number: int = 20) -> List[float]: scorefxn = pyrosetta.get_fa_scorefxn() backrub = pyrosetta.rosetta.protocols.backrub.BackrubMover() monégasque = pyrosetta.rosetta.protocols.monte_carlo.GenericMonteCarloMover(maxtrials=trials, max_accepted_trials=trials, # gen.max_accepted_trials() = 0 task_scaling=5, # gen.task_scaling() mover=backrub, temperature=temperature, sample_type='low', drift=True) monégasque.set_scorefxn(scorefxn) # find most deviant rs = [] for i in range(replicate_number): variant = original.clone() monégasque.apply(variant) if monégasque.accept_counter() > 0: rs.append(pr_scoring.bb_rmsd(variant, original)) else: rs.append(float('nan')) return rs
Conclusion
Hopefully this exhaustive, and exhausting, overview of the movers and functions in PyRosetta, has adequately showcases what operations can be done in PyRosetta to better tune and analyse variants from RFdiffusion.
Footnote: ProteinMPNN helper as API functions
As mentioned, my designs had different lengths and I wanted to create the definitions for ProteinMPNN via my script as it was reading from tar.gz archives and other details. But briefly, I parsed the RFdiffusion skeletons thusly:
no_fix = [] definitions = [] global_fixed_chains: Dict[str, List[List[str]]] = {} global_fixed_positions = {} for path in paths: # ## Load data name = path.stem pdbblock = get_ATOM_only(path.read_text()) # ## Parse definition (i.e. coordinates) definition = parse_PDBblock(pdbblock, name) definitions.append(definition) # ## Fixed chain fixed_chains = define_fixed_chains(definition, designed_chain_list='A') # chain A is designed global_fixed_chains[name] = fixed_chains # ## Fixed pos sequence = get_chainA_sequence(pdbblock) # assuming the sequence to change is chain A # The skeleton will have stretches of glycines where the design was made: masked = re.sub(r'(G{3,})', lambda match: '_' * len(match.group(1)), sequence) fixed_list = [[i for i, r in enumerate(masked) if r == '_']] fixed_positions = define_unfixed_positions(definition, ['A'], fixed_list) global_fixed_positions[name] = fixed_positions # write out definitions, global_fixed_chains, global_fixed_positions with open(chains_definitions_path, 'w') as fh: # only chains_definitions.jsonl is a JSONL for definition in definitions: fh.write(json.dumps(definition) + '\n') with open(fixed_chains_path, 'w') as fh: json.dump(global_fixed_chains, fh) with open(fixed_positions_path, 'w') as fh: json.dump(global_fixed_positions, fh)
Where the functions are:
""" This is a minor refactoring of the original code, with the following changes: `parse_PDBblock(pdbblock: str, name: str, chain_alphabet, ca_only=False)` accepts a PDB block ``pdbblock`` and a name ``name`` and returns a dictionary that forms a line for the JSONL. The code is mostly kept as was, with minor chances. """ import json from typing import List, Sequence, Dict, Any import numpy as np alpha_1 = list("ARNDCQEGHILKMFPSTWYV-") states = len(alpha_1) alpha_3 = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'GAP'] aa_1_N = {a: n for n, a in enumerate(alpha_1)} aa_3_N = {a: n for n, a in enumerate(alpha_3)} aa_N_1 = {n: a for n, a in enumerate(alpha_1)} aa_1_3 = {a: b for a, b in zip(alpha_1, alpha_3)} aa_3_1 = {b: a for a, b in zip(alpha_1, alpha_3)} def get_ATOM_only(pdbblock: str) -> str: """ This gets all ATOM, regardless of name and chain """ return '\n'.join([line for line in pdbblock.splitlines() if line.startswith('ATOM')]) def get_chainA_sequence(pdbblock: str) -> str: sequence = '' residues_seen = set() for line in pdbblock.splitlines(): if line.startswith("ATOM") and " CA " in line and " A " in line: res_info = line[17:26] # Residue name and number for uniqueness if res_info not in residues_seen: residues_seen.add(res_info) res_name = line[17:20].strip() sequence += three_to_one.get(res_name, '?') return sequence three_to_one = { 'ALA': 'A', 'CYS': 'C', 'ASP': 'D', 'GLU': 'E', 'PHE': 'F', 'GLY': 'G', 'HIS': 'H', 'ILE': 'I', 'LYS': 'K', 'LEU': 'L', 'MET': 'M', 'ASN': 'N', 'PRO': 'P', 'GLN': 'Q', 'ARG': 'R', 'SER': 'S', 'THR': 'T', 'VAL': 'V', 'TRP': 'W', 'TYR': 'Y' } def AA_to_N(x): # ["ARND"] -> [[0,1,2,3]] x = np.array(x); if x.ndim == 0: x = x[None] return [[aa_1_N.get(a, states - 1) for a in y] for y in x] def N_to_AA(x): # [[0,1,2,3]] -> ["ARND"] x = np.array(x); if x.ndim == 1: x = x[None] return ["".join([aa_N_1.get(a, "-") for a in y]) for y in x] def inner_parse_PDBblock(pdbblock, atoms=['N', 'CA', 'C'], chain=None) -> tuple: ''' input: pdb_filename = PDB filename atoms = atoms to extract (optional) output: (length, atoms, coords=(x,y,z)), sequence ''' xyz, seq, min_resn, max_resn = {}, {}, 1e6, -1e6 for line in pdbblock.split("\n"): if line[:6] == "HETATM" and line[17:17 + 3] == "MSE": line = line.replace("HETATM", "ATOM ") line = line.replace("MSE", "MET") if line[:4] == "ATOM": ch = line[21:22] if ch == chain or chain is None: atom = line[12:12 + 4].strip() resi = line[17:17 + 3] resn = line[22:22 + 5].strip() x, y, z = [float(line[i:(i + 8)]) for i in [30, 38, 46]] if resn[-1].isalpha(): resa, resn = resn[-1], int(resn[:-1]) - 1 else: resa, resn = "", int(resn) - 1 # resn = int(resn) if resn < min_resn: min_resn = resn if resn > max_resn: max_resn = resn if resn not in xyz: xyz[resn] = {} if resa not in xyz[resn]: xyz[resn][resa] = {} if resn not in seq: seq[resn] = {} if resa not in seq[resn]: seq[resn][resa] = resi if atom not in xyz[resn][resa]: xyz[resn][resa][atom] = np.array([x, y, z]) # ^^ end of xyz loop # convert to numpy arrays, fill in missing values seq_, xyz_ = [], [] try: for resn in range(min_resn, max_resn + 1): if resn in seq: for k in sorted(seq[resn]): seq_.append(aa_3_N.get(seq[resn][k], 20)) else: seq_.append(20) if resn in xyz: for k in sorted(xyz[resn]): for atom in atoms: if atom in xyz[resn][k]: xyz_.append(xyz[resn][k][atom]) else: xyz_.append(np.full(3, np.nan)) else: for atom in atoms: xyz_.append(np.full(3, np.nan)) return np.array(xyz_).reshape(-1, len(atoms), 3), N_to_AA(np.array(seq_)) except TypeError: return 'no_chain', 'no_chain' def get_chain_alphabet(): init_alphabet = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z'] extra_alphabet = [str(item) for item in list(np.arange(300))] return init_alphabet + extra_alphabet def parse_PDBblock(pdbblock: str, name: str, ca_only=False): chain_alphabet = get_chain_alphabet() my_dict = {} s = 0 concat_seq = '' for letter in chain_alphabet: if ca_only: sidechain_atoms = ['CA'] else: sidechain_atoms = ['N', 'CA', 'C', 'O'] xyz, seq = inner_parse_PDBblock(pdbblock=pdbblock, atoms=sidechain_atoms, chain=letter) if type(xyz) != str: concat_seq += seq[0] my_dict['seq_chain_' + letter] = seq[0] coords_dict_chain = {} if ca_only: coords_dict_chain['CA_chain_' + letter] = xyz.tolist() else: coords_dict_chain['N_chain_' + letter] = xyz[:, 0, :].tolist() coords_dict_chain['CA_chain_' + letter] = xyz[:, 1, :].tolist() coords_dict_chain['C_chain_' + letter] = xyz[:, 2, :].tolist() coords_dict_chain['O_chain_' + letter] = xyz[:, 3, :].tolist() my_dict['coords_chain_' + letter] = coords_dict_chain s += 1 my_dict['name'] = name my_dict['num_of_chains'] = s my_dict['seq'] = concat_seq if s < len(chain_alphabet): return my_dict else: raise Exception('Too many chains') def define_fixed_chains(chains_definition: Dict[str, Any], designed_chain_list=('A',)): """ The fixed chain definition file is a JSON dictionary of name to tuple/list of two: designed_chain_list and fixed_chain_list """ all_chain_list = [item[-1:] for item in list(chains_definition) if item[:9] == 'seq_chain'] # ['A','B', 'C',...] fixed_chain_list = [letter for letter in all_chain_list if letter not in designed_chain_list] # fix/do not redesign these chains return list(designed_chain_list), fixed_chain_list def define_global_fixed_chains(chains_definitions, global_designed_chain_list=('A',)): return {chains_definition['name']: define_fixed_chains(chains_definition, global_designed_chain_list) for chains_definition in chains_definitions} def define_fixed_positions(chains_definition: Dict[str, Any], designed_chain_list: Sequence[str], fixed_list: Sequence[Sequence[int]]): all_chain_list = [item[-1:] for item in list(chains_definition) if item[:9] == 'seq_chain'] fixed_position_dict = {} for i, chain in enumerate(designed_chain_list): fixed_position_dict[chain] = fixed_list[i] for chain in all_chain_list: if chain not in designed_chain_list: fixed_position_dict[chain] = [] return fixed_position_dict def define_unfixed_positions(chains_definition, designed_chain_list: Sequence[str], unfixed_list: Sequence[Sequence[int]]): """ This will be an entry in a dictionary passed to ``chain_id_jsonl`` (Misnomer: it's not a JSONL, but a JSON, only ``jsonl_path`` is) :param chains_definition: :param designed_chain_list: :param unfixed_list: :return: """ all_chain_list = [item[-1:] for item in list(chains_definition) if item[:9] == 'seq_chain'] fixed_position_dict = {} for chain in all_chain_list: seq_length = len(chains_definition[f'seq_chain_{chain}']) all_residue_list = (np.arange(seq_length) + 1).tolist() if chain not in designed_chain_list: fixed_position_dict[chain] = all_residue_list else: idx = np.argwhere(np.array(designed_chain_list) == chain)[0][0] fixed_position_dict[chain] = list(set(all_residue_list) - set(unfixed_list[idx])) return fixed_position_dict