cartesian delta delta G

Hey!

I recently implemented the algorithm for the cartesian ddG (DiMaio, https://www.rosettacommons.org/docs/latest/cartesian-ddG), following the source code (https://github.com/RosettaCommons/main/blob/master/source/src/apps/public/ddg/cartesian_ddg.cc). In principle it does select the next sequential neighbors, allowing chi-relax and then the neighbors in space to do SC relax. 

A note here: I didn't implement the interface mode which is in the original code. Also, this is the legacy version (see c++ code above) - but it should definitely solve your problem with the selection ;) 

Another note: you can make the code a bit more easy and remove the repacking (which in the end is just a mutation) with this function (I just wanted to play a bit with ResfileCommandOperation: pyrosetta.toolbox.mutants.mutate_residue(repacked_pose, int(target_position), mutant, pack_radius=repack_radius, pack_scorefxn=sfxn)

Also, I did not compare the performance of it with the original C++ (this should def. be done) - so some errors might still be in the code - if you find them, I would appreciate a note!

Cheers,

Johanna

def mutate_repack_func4(pose, target_position, mutant, repack_radius, sfxn, ddg_bbnbrs=1, verbose=True, cartesian=False, max_iter=None):
    import time
    from pyrosetta.rosetta.core.pack.task import operation

    logger.warning("Interface mode not implemented (should be added!)")
    
    if cartesian:
        sfxn.set_weight(pyrosetta.rosetta.core.scoring.ScoreTypeManager.score_type_from_name('cart_bonded'), 0.5)
        #sfxn.set_weight(atom_pair_constraint, 1)#0.5
        sfxn.set_weight(pyrosetta.rosetta.core.scoring.ScoreTypeManager.score_type_from_name('pro_close'), 0)
        #logger.warning(pyrosetta.rosetta.basic.options.get_boolean_option('ex1'))#set_boolean_option( '-ex1', True )
        #pyrosetta.rosetta.basic.options.set_boolean_option( 'ex2', True )
    
    #Cloning of the pose including all settings
    working_pose = pose.clone()

    #Select mutant residue
    mutant_selector = pyrosetta.rosetta.core.select.residue_selector.ResidueIndexSelector(target_position)
    
    #Select all except mutant
    all_nand_mutant_selector = pyrosetta.rosetta.core.select.residue_selector.NotResidueSelector()
    all_nand_mutant_selector.set_residue_selector(mutant_selector)

    #Select neighbors with mutant
    nbr_or_mutant_selector = pyrosetta.rosetta.core.select.residue_selector.NeighborhoodResidueSelector()
    nbr_or_mutant_selector.set_focus(str(target_position))
    nbr_or_mutant_selector.set_distance(repack_radius)
    nbr_or_mutant_selector.set_include_focus_in_subset(True)

    #Select mutant and it's sequence neighbors
    seq_nbr_or_mutant_selector = pyrosetta.rosetta.core.select.residue_selector.PrimarySequenceNeighborhoodSelector(ddg_bbnbrs, ddg_bbnbrs, mutant_selector, False)            

    #Select mutant, it's seq neighbors and it's surrounding neighbors
    seq_nbr_or_nbr_or_mutant_selector = pyrosetta.rosetta.core.select.residue_selector.OrResidueSelector()
    seq_nbr_or_nbr_or_mutant_selector.add_residue_selector(seq_nbr_or_mutant_selector)
    seq_nbr_or_nbr_or_mutant_selector.add_residue_selector(nbr_or_mutant_selector)    

    if verbose:
        logger.info(f'mutant_selector: {pyrosetta.rosetta.core.select.residue_selector.selection_positions(mutant_selector.apply(working_pose))}')
        logger.info(f'all_nand_mutant_selector: {pyrosetta.rosetta.core.select.residue_selector.selection_positions(all_nand_mutant_selector.apply(working_pose))}')
        logger.info(f'nbr_or_mutant_selector: {pyrosetta.rosetta.core.select.residue_selector.selection_positions(nbr_or_mutant_selector.apply(working_pose))}')
        logger.info(f'seq_nbr_or_mutant_selector: {pyrosetta.rosetta.core.select.residue_selector.selection_positions(seq_nbr_or_mutant_selector.apply(working_pose))}')
        logger.info(f'seq_nbr_or_nbr_or_mutant_selector: {pyrosetta.rosetta.core.select.residue_selector.selection_positions(seq_nbr_or_nbr_or_mutant_selector.apply(working_pose))}')
        

    #Mutate residue and pack rotamers before relax
    if list(pose.sequence())[target_position-1] != mutant:
        #generate packer task
        tf = TaskFactory()
        tf.push_back(operation.InitializeFromCommandline())
        tf.push_back(operation.IncludeCurrent())
        
        #Set all residues except mutant to false for design and repacking
        prevent_repacking_rlt = operation.PreventRepackingRLT()
        prevent_subset_repacking = operation.OperateOnResidueSubset(prevent_repacking_rlt, all_nand_mutant_selector, False )
        tf.push_back(prevent_subset_repacking)
        
        #Assign mutant residue to be designed and repacked
        resfile_comm = pyrosetta.rosetta.protocols.task_operations.ResfileCommandOperation(mutant_selector, f"PIKAA {mutant}")
        resfile_comm.set_command(f"PIKAA {mutant}")
        tf.push_back(resfile_comm)

        #Apply packing of rotamers of mutant
        packer = pyrosetta.rosetta.protocols.minimization_packing.PackRotamersMover()
        packer.score_function(sfxn)
        packer.task_factory(tf)
        logger.warning(tf.create_task_and_apply_taskoperations(working_pose))
        packer.apply(working_pose)

    #allow the movement for bb for the mutant + seq. neighbors, and sc for neigbor in range, seq. neighbor and mutant
    movemap = pyrosetta.rosetta.core.select.movemap.MoveMapFactory()
    movemap.all_jumps(False)
    movemap.add_bb_action(pyrosetta.rosetta.core.select.movemap.mm_enable, seq_nbr_or_mutant_selector)
    movemap.add_chi_action(pyrosetta.rosetta.core.select.movemap.mm_enable, seq_nbr_or_nbr_or_mutant_selector)
    
    #for checking if all has been selected correctly
    if verbose:
        mm  = movemap.create_movemap_from_pose(working_pose)
        logger.info(mm)

    #Generate a TaskFactory
    tf = TaskFactory()
    tf.push_back(operation.InitializeFromCommandline())
    tf.push_back(operation.IncludeCurrent())
    #tf.push_back(operation.NoRepackDisulfides())

    #prevent all residues except selected from design and repacking
    prevent_repacking_rlt = operation.PreventRepackingRLT()
    prevent_subset_repacking = operation.OperateOnResidueSubset(prevent_repacking_rlt, seq_nbr_or_nbr_or_mutant_selector, True )
    tf.push_back(prevent_subset_repacking)

    # allow selected residues only repacking (=switch off design)
    restrict_repacking_rlt = operation.RestrictToRepackingRLT()
    restrict_subset_repacking = operation.OperateOnResidueSubset(restrict_repacking_rlt , seq_nbr_or_nbr_or_mutant_selector, False)
    tf.push_back(restrict_subset_repacking)


    #Perform a FastRelax
    fastrelax = pyrosetta.rosetta.protocols.relax.FastRelax()
    fastrelax.set_scorefxn(sfxn)
    
    if cartesian:
        fastrelax.cartesian(True)
    if max_iter:
        fastrelax.max_iter(relax_iter)
        
    fastrelax.set_task_factory(tf)
    fastrelax.set_movemap_factory(movemap)
    fastrelax.set_movemap_disables_packing_of_fixed_chi_positions(True)
    
    if verbose:
        logger.info(tf.create_task_and_apply_taskoperations(working_pose))
    fastrelax.apply(working_pose)
    return working_pose

Hey David,

good catch! I copied the functions over and forgot to edit this one!

With the iteration of 3 (from the documentation), it is meant to call the  mutate_repack_func4() function 3 times, in order to get statistics.

The max_iter value is with respect to the fastrelax cycles, so (as you mentioned) definitely more than 3! If you don't set the max_iter value, the default for fastrelax is used, which has been in our hands more than enough - but feel free to play around with it and vary it.

Yes, indeed, the github link only works if you are signed as Rosetta member in but you can see the origin of the file from the link and can check it out in your source code ;)

Hope it makes more sense now and reach out if it doesn't!

Cheers,

Johanna