What is the best way to make a residue-specific SASA constraint?

[Anonymous]

Hi everyone,

I am trying to model (not design) a protein in an alternate state, where a phenylalanine is burried in the native state, and solvent exposed in the one I am modeling. 

I would like to use a constraint, or some or other method, that specifies the solvent accessibility of that PHE, while running a relax or remodel application. I would prefer to not use a list of nearest-neighbor distance constraints, since there are quite a few of them, the residues and distances vary depending on the exposure of the PHE, and it is the identity and distances of the nearest neighbors that I am trying to figure out based on the exposure of the PHE.

Any suggestions would be greatly appreciated!

Thanks,

Jan

[Anonymous]

I know you said remodel/relax application and not script, but the SASA selection is available in scripting.

https://www.rosettacommons.org/docs/latest/scripting_documentation/RosettaScripts/TaskOperations/taskoperations_pages/LayerDesignOperation

The bad news is that there are few cascading effects to deal with:

The <Layer> elements go within <RESIDUE_SELECTORS> element. In the <MoveMap> of <FastRelax> (or similar mover)  there is written that you can add <ResidueSelector>. This in my hands is not true. No idea why the documentation says you can, but you get a “Error: Element ‘ResidueSelector’: This element is not expected. Expected is one of ( Jump, Chain, Span ).” if you do. So you have to use a movemap factory, which can.

So you have to add:

<MOVE_MAP_FACTORIES>

        <MoveMapFactory name=”moove” bb=”0″ chi=”0″>

            <Backbone residue_selector=”layerz” />

            <Chi residue_selector=”layerz” />

        </MoveMapFactory>

</MOVE_MAP_FACTORIES>

Which is called by <FastRelax name=”chillax” repeats=”1″ movemap_factory=”moove” cartesian=”1″ />.

Also, cartesian=”1″ is needed or the constraint gets ignored. But to get cartesian mode to work you have to add -beta_nov16_cart to your run. 

$ROSETTAPATH/source/bin/rosetta_scripts.static.macosclangrelease -s something.pdb -parser:protocol yourscript.xml -beta_nov16_cart

[Anonymous]

Sorry to double post but I remembered that I had a snippets of code that could something similar to what you asked* and then I re-read the question and realised you know a priori what the phenylalanine is….

Rosetta Relax has an option `-in:file:movemap`, where you can specify explicitly which residues you want to move both backbone (bb) or sidechain (chi). A movemap is also an element in Rosetta scripts and in Pyrosetta.  The movemap file will have to be made external and not in Rosetta.

Note that the values are in pose numbering and not PDB and using the 23A kind of notation does not work IIRC.

Pyrosetta has a really nice function called pdb2pose that does these conversion for you. Here is a chunk of code as an example (same thing as using the application).

import pyrosetta



pyrosetta.init()



pose = pyrosetta.pose_from_pdb("my_protein.pdb")

pdb2pose = pose.pdb_info().pdb2pose



scorefxn = pyrosetta.get_fa_scorefxn()

# or whatever.

# scorefxn = pyrosetta.rosetta.core.scoring.ScoreFunctionFactory.create_score_function("ref2015_cart.wts")



relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, 3)

movemap = pyrosetta.MoveMap()

for chain, resi in targets:

posed = pdb2pose(chain=chain, res=resi)

print(chain, res, '-->', posed)

movemap.set_bb(posed, True)

movemap.set_chi(posed, True)

relax.set_movemap(movemap)

relax.apply(pose)

print('Score', scorefxn(pose))

pose.dump_pdb('output.pdb')

*) Not that it matters but the code I was going to post was a python script using pymol that I quickly adapted:

import pymol2



#Wikipedia Tien et al. 2013 (emp.)

maxASA = {'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}





with pymol2.PyMOL() as pymol:

pymol.cmd.set('dot_solvent', 1)

pymol.cmd.set('dot_density', 3)



def is_surface(atom) -> bool:

sasa = pymol.cmd.get_area(f'resi {atom.resi} and chain {atom.chain} and name CA')

rsa = sasa/maxASA[self.mutation.from_residue]

return True if rsa < 0.2 else False #classic



pymol.cmd.load('my_PDB')

atoms = pymol.cmd.get_model(f"resn PHE")

for atom in atoms.atom:

if is_surface(atom):

resi = int(re.match('d+', resi).group())

targets.append((resi, atom.chain))



print(targets)

N

[Author]

Posts