RMSD with jump in one pose

[NooraAz]

I want to get RMSD between a native pose from PDB, and a pose made with AF2. The native has a jump in the middle, missing 17 residues. So When I use this code the get RMSD, it fails with the following error:

rmsd = rosetta.core.simple_metrics.metrics.RMSDMetric(native)
rmsd.set_run_superimpose(True)
pose_rmsd = rmsd.calculate(pose)

File: /Volumes/scratch/b3.r/rosetta/commits/main/source/src/core/pose/Pose.cc:975
[ ERROR ] UtilityExitException
ERROR: Pose::residue_type( Size const seqpos ): variable seqpos is out of range!

How should I get RMSD in these cases?

[roccomoretti]

RMSDMetric can take a series of ResidueSelectors (In PyRosetta it would be the set_residue_selector(), set_residue_selector_reference(), set_residue_selector_super() and set_residue_selector_super_reference() methods.

If you want to calculate the rmsd over a subset of residues in the Pose(s), then you can set the residue selector to use to select the residues you want to use for the calculation. By default the selector passed to set_residue_selector() is used for both the primary and the reference/native structure, but in your case I'd probably recommend setting separate residue selectors for each. (The super variants are there if you want to superimpose on a different set of residues than you calculate the rmsd over.)

The only "trick" is that the sets of residues have to be ordered similarly in the pose (once the set of residues is selected, they're matched up one-to-one in the order in their respective poses). For example, if you have an ABCD structure and an ADBC structure, and you want to calculate the rmsd on just chains B & D, this way won't work properly. -- But in your case where you're just missing some residues but the poses are otherwise the same, things should work.

[NooraAz]

Thanks, Rocco.
I ended up using this script, for future reference:

selector = rosetta.core.select.residue_selector.ResidueIndexSelector('1-227,245-497')
whole_pose_selector = rosetta.core.select.residue_selector.TrueResidueSelector()

rmsd = rosetta.core.simple_metrics.metrics.RMSDMetric(pose)

rmsd.set_residue_selector(whole_pose_selector)
rmsd.set_residue_selector_reference(selector)

rmsd.set_run_superimpose(True)

pose_rmsd = rmsd.calculate(native)

I switched pose and native, because the native is the protein with the gap, and I could only add a second selector for the "reference" pose.

The result is slightly different from pymol though (~0.01), which I don't know why yet.

[roccomoretti]

Note that PyMol sometimes does a "filtered" rmsd, where it throws out residues which it feels are outliers, so you might not be comparing the rmsd across the same set. Also 0.01 Angstrom sounds like it's potentially due to slight differences in calculating the alignment (even over the same residues.)

[NooraAz]

I use cycles=0 to avoid throwing out residues. But yeah I also guess the algorithm is a bit different (?). I'm just curious to know the difference.

[ajasja]

I would not worry about a 0.01 A difference. The standard pymol align does a sequence match first, maybe if you used super you would get the exact same number. Also perhaps one or the other just takes backbone atoms vs sidechain atoms etc.