how to use the PyRosetta to obtain protein docking results with contact map constraint？

[onlyonewater]

Hi, PyRosetta Team. Thanks for releasing this awesome tool for the community. I am a newcomer to the protein computing field. Now, I am working on protein docking. I have an initial protein pose with a slightly accurate contact map. I want to use the contact map as a constraint on the initial protein pose to obtain a more accurate protein structure. How can I use PyRosetta to achieve this goal? Could you give me some advice? Thanks!!!

[roccomoretti]

Contacts are usually input into PyRosetta as Rosetta constraints. That is, you modify the energy potential such that you add a penalty when a potential docked structure violates the contact map you have.

In particular, for most contacts, you're looking for AtomPair constraints. You probably want to use the FlatHarmonic function type, which allows you to specify a range (e.g. 0-8 Ang) over which there is no penalty, and then a quadraticly-increasing penalty outside that range.

For protocols like protein-protein docking, which go through both a full atom and low-resolution ("centroid") stage, you'll want to base your constraints off of backbone heavyatoms (typically Calpha), as the centroid stage doesn't represent sidechains atomistically.

Using constraints in (Py)Rosetta is a two-step process. The first is to add the constraint definitions to the Pose itself. The second is to make sure you're using a scorefunction with the appropriate constraint scoreterms turned on.

The easiest way to add constraints to a pose is to write a Rosetta constraint file with your contact map information. A constraint file is nothing more than a text file with one line per constraint you want. For example,

AtomPair CA 34A CA 123B FLAT_HARMONIC 0 1.5 8.0

defines a constraint which looks at the distance between the Calpha atoms on residues 34A and 123B and applies a scoring penalty if they're greater than 8.0 Ang apart. (The 1.5 changes how steeply the penalty rises. The units are in Angstrom, and smaller values are a stiffer penalty.)

You can then apply that constraint file to the pose with core.scoring.constraints.ConstraintIO.read_constraints(cstfilename, pose.constraint_set(), pose)

Alternatively, you can set up the constraints programmatically:

func = core.scoring.constraints.func.FlatHarmonicFunc( 0.0, 1.5, 8.0 )
cst = core.scoring.constraints.AtomPairConstraint( core.id.AtomID( atmno1, resi1), core.id.AtomID( atmno2, resi2), func )
pose.add_constraint( cst )

where atomno & resi are the (integer) residue id and atom number of the atoms you want to constrain.

[onlyonewater]

Thanks!! I will have a try!!

[onlyonewater]

It seems that if there is no contact in the protein, i.e., the contact map is zero, using the command like core.scoring.constraints.ConstraintIO.read_constraints(cstfilename, pose.constraint_set(), pose) would lead to the program stuck, like this:

how to solve this situation?

[roccomoretti]

What's the constraint file look like? Is it an empty file? -- I don't know why things are hanging for you, but one workaround is to doublecheck the constraint file before loading to see if it's empty or not. (If you don't have any contact constraints, don't bother loading it -- not loading constraints functions identically to loading an empty set of constraints.)

[onlyonewater]

yes, the constraint file is empty since there is no contact in the protein complex, and you are right, I skipped loading the constraint file and the program can run normally.

[onlyonewater]

and now I am using the PyRosetta to adjust the protein complex structure based on the ground-truth contact map, however, the adjusted structures are worse than the original one, I do not know why, the following is my demo:

    contact_map = np.load(contact_map_file)

    pose = pyrosetta.pose_from_pdb(pdb_file)

    create_constraints_file(contact_map, cons_file)
    ConstraintIO.read_constraints(cons_file, pose.constraint_set(), pose)

    scorefxn = get_fa_scorefxn()
    scorefxn.set_weight(pyrosetta.rosetta.core.scoring.atom_pair_constraint, 2.0)
    initial_score = scorefxn(pose)
    print(f"Initial Score: {initial_score}")

    mm = pyrosetta.MoveMap()
    mm.set_bb(True)
    mm.set_chi(True)

    min_mover = MinMover()
    min_mover.movemap(mm)
    min_mover.score_function(scorefxn)
    min_mover.apply(pose)

    pose.dump_pdb("{}_output.pdb".format(pdb_name))

    optimized_score = scorefxn(pose)
    print(f"Optimized Score: {optimized_score}")

could you give me some advice?

[roccomoretti]

The general suggestion with debugging constraints is to adjust the weight of the constraints and see how that affects things. Since constraints are a penalty added to the scorefunction, there's often a balance between satisfying the constraints and making the structure look good according to other parameters. It may take a bit in order to reshape the energy landscape enough such that Rosetta sampling recognizes the gradient to the constrained structure. Conversely, having the constraints on too high will shift things such that Rosetta will only try to satisfy the constraints, and will ignore the other score terms. Generally there's a constraint weight which balances the two opposing extremes, but you may have to guess and check a number of them until you find the right balance.

The other thing to do during this process is to double check your constraint settings, your constraint energies, and how well your structures satisfy the constraints. There's a pose.constraint_set().show_violations(...) function which potentially would help debugging. Though it might be hard to call in PyRosetta, as it needs a pyrosetta.rosetta.std.ostream object as the first parameter. -- Another way to check violations is to just look at the structure (e.g. in PyMol of ChimeraX) and measure the same distances the constraints are measuring. Does everything look like you expect? Are the constrained distances getting better like you expect?

Another thing to keep in mind when investigating is to clarify in your mind what you mean by "the structures are worse" -- are you talking about binding in the correct location? Satisfying contact maps? How well it obeys the constraint definitions? Is it improving in energy? (These all may be related but are subtly distinct.) Or is it non-contact related, like awkward sidechains, bad backbone structures or missing hydrogen bonds. Each might have a subtly different approach for investigating and fixing.

[onlyonewater]

Thanks for your detailed reply. I use the DockQ and its corresponding success rate (DockQ > 0.23) as the metric to measure the final structure quality, the structures get worse means that although the DockQ of the beginning structure is low (almost 0.12), after using the contact map as constraint, the final structure gets lower DockQ than the original structure (around 0.09). So, I do not know the reason.

After reading your answer, I think increasing/decreasing the constraint weight might not be the optimal solution since now I have adjusted the weight to 20.0/1.0 in the score function, and the DockQ does not get higher (almost the same as the original value). I am not sure if the above demo exists bug or if I missed some steps.