How to set up an electron density scorefunction in Pyrosetta

Energy minimising structures in Rosetta/Pyrosetta is essential to avoid artifactual results. Say a mutation is introduced and in the protocol the neighbourhood is repacked: if the structure is not energy minimised properly the neighbourhood repacking step will spuriously reward the mutation a very negative ∆∆G. One worry is that the energy minimisation is not faithful to the crystal structure. This argument has two sides, on one the fudgey force fields in Rosetta do not truly model the chemical interactions while on the other crystal packing may be unnatural. Both points have merit. After all Rosetta does use implicit water, which do not behave like the stripped crystallographic waters and some residues may have non-standard protonations etc. But if one wants one can use a scorefunction that is weighted by the electron density map and here is how.

Getting a map

The first thing to do is to get a CCP4 map. RCSB PDB only provides mtz and fo-fc maps (dsn6 format) and no CCP4 maps. However, the all too often overlooked PDBe provides CCP4 maps —Pyrosetta documentation uses the name MRC maps, which is a CCP4 but for electron microscopy and is nearly identically and can read both.

scorefxnED = pyrosetta.get_fa_scorefxn()
ED = pyrosetta.rosetta.core.scoring.electron_density.getDensityMap('1ubq.ccp4')
print(ED.matchPose(pose))
sdsm = pyrosetta.rosetta.protocols.electron_density.SetupForDensityScoringMover()
sdsm.apply(pose)
elec_dens_fast = pyrosetta.rosetta.core.scoring.ScoreType.elec_dens_fast
scorefxnED.set_weight(elec_dens_fast, 30)

The scoring namespace keeps track of the DensityMap instance, hence why this is not passed to SetupForDensityScoringMover, which is rather odd. The score results in no longer in kcal/mol, but severely weighted in favour of the electron density —but after all scorefunctions with tinkered weights are for movers and not for people.

print(scorefxnED.weights()[elec_dens_fast])
print(scorefxnED(pose))

Energy minimisation in Rosetta is most often done with the FastRelax mover. The thorough option in the standolone relax application does 15 cycles. I personally do a compromise, wherein I ramp down the weight of the density map in 3 steps of 5 cycles and get very nice results.

for w in (30, 20, 10):
    scorefxnED.set_weight(elec_dens_fast, w)
    relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxnED, 5)
    relax.apply(pose)
reg_scorefxn = pyrosetta.get_fa_scorefxn()
print('How good is the density match:', ED.matchPose(pose))
print(reg_scorefxn.get_name(), reg_scorefxn(pose))

To give it a test I suggest using the ShakeStructureMover mover, which wobbles the structure by an amount specified. I am not sure what it is actually for —it appears to be a lot more complex than just shaking the protein in internal coordinate mode—, but for sure it is cool for testing little things like this, say:

shaker = pyrosetta.rosetta.protocols.simple_moves.ShakeStructureMover()
shaker.set_mc_temperature(5)
shaker.apply(pose)
pymol.pymol_name('shaken')
pymol.apply(pose)
... relax now as above
pymol.pymol_name('fixed')
pymol.apply(pose)

Footnote: LoadDensityMapMover

In various scripts, the mover LoadDensityMapMover is used to load the map.

map_mover = pyrosetta.rosetta.protocols.cryst.LoadDensityMapMover(map_filename)
map_mover.apply(pose)

This loads it globally like getDensityMap, but doing a single test (using PDB:1UBQ and 15 cycles of FastRelax (not cartesian) and a map weight of 30) it seems to give slightly worse fits than getDensityMap, with scores within error. So it seems to work fine with getDensityMap, but not as a replacement for it.

getDensityMap	SetupForDensityScoringMover	LoadDensityMapMover	ori RSMD	ED fit
✗	✗	✓	1.183	NA
✗	✓	✓	1.140	0.794
✓	✓	✗	0.751	0.769
✓	✓	✓	0.812	0.769