I modelled by own mutate residue function based on the one provided in the PyRosetta toolbox scripts. It sets up a task which specifices that only residues within a certain distance from the mutation are to be repacked, as well as the small molecule ligand which is docked in binding pocket. This is then given to the packer (PackRotamersMover) and applied to the pose. But when I inspected the decoys, of which there were thousands, it was obvious that the rotamer of the ligand in the output was identical to the input structure's rotamer: i.e. the PackRotamersMover was not actually repacking the ligand in the binding pocket.
However, if I print the task before sending it to the packer, it indicates that it should be mutating residue 108 to serine, repacking surrounding residues (i.e. methionine at 109) as well as repacking the ligand, which I've called PCZ here:
>>> resid pack? design? allowed_aas
>>> 1 FALSE FALSE
>>> : : :
>>> 108 TRUE TRUE SER
>>> 109 TRUE FALSE MET
>>> : : :
>>> 308 TRUE FALSE PCZ
My initial thoughts are that perhaps the PackRotamersMover doesn't repack non-canonical AAs unless some specific option is set for that residue, but its been hard to test this. I've also checked the params file, and it does indeed point to the right PDB_ROTAMERS file, and the packer does also print a line to the terminal saying that it has successfully loaded so and so many rotamers for the given position of the ligand whenever the function is called. Also, when I run a ligand dock protocol with the DockMCMProtocol() with exactly the same input, the decoys all have a range of different rotamers in them, so I don't think its a problem with the input files.
Any thoughts as to why this might be happening, or a potential solution? A minimal version of the script is below (pretty much identical to the Toolbox script).
def mutate_residue_1(pose, mutant_position, mutant_aa, pack_radius, pack_scorefxn): new_pose = Pose() # copy pose new_pose.assign(pose) task = standard_packer_task(new_pose) # Each number between between 1-20 correspond to a unique AA. # Get the mutant AA's code. mutant_aa = aa_from_oneletter_code(mutant_aa) # Make a vector of True's False's, where each index corresponds to an AA, # indicating whether to include it or not during design. We only want the # mutant AA. aa_bool = utility.vector1_bool() for i in xrange(1,21): aa_bool.append(i == mutant_aa) task.nonconst_residue_task( mutant_position).restrict_absent_canonical_aas(aa_bool) # Repack those within a certain pack_radius. centre = pose.residue(mutant_position).nbr_atom_xyz() for i in xrange(1,pose.total_residue() + 1) : if (i == mutant_position) or pose.residue(i).is_ligand(): continue; elif (centre.distance_squared(new_pose.residue(i).nbr_atom_xyz()) < pack_radius**2): aa_bool = utility.vector1_bool() curr_aa = pose.residue(i).aa() for j in xrange(1,21): aa_bool.append(j == curr_aa) task.nonconst_residue_task(i).restrict_absent_canonical_aas(aa_bool) else: task.nonconst_residue_task(i).prevent_repacking() # Apply the mutation and pack nearby residues. packer = PackRotamersMover(pack_scorefxn, task) packer.apply(new_pose) return new_pose