Code
----
---- For each rotamer position of each sidechain, does a flip, freeze, shake, unfreeze, wiggle, shake, wiggle.
---- Skips the rest of the operation sequence and continues with the next rotamer position if the first shake does not improve the score.
----
behaviour = behavior
-- Amino acids. Hydrophobics first.
single_letter_codes = { "g","a","v","l","i","m","f","w","p","s","t","c","y","n","q","d","e","k","r","h"}
three_letter_codes = { "Gly","Ala","Val","Leu","Ile","Met","Phe","Trp","Pro","Ser","Thr","Cys","Tyr","Asn","Gln","Asp","Glu","Lys","Arg","His"}
---- option default values ----
options = {
-- post_scramble = false: no sidechain scrambling; won't let walk juice.
-- post_scramble = true: brief scramble iff close to the best score
["post_scramble"] = true,
["scrambleThreshold"] = 2,
["min_snap_loss"] = 60,
-- not a valid segment index => pseudorandom start segment computed from the current score
["start_segment"] = 0,
["quicksaves"] = {
-- stores the recent best structure
["best"] = 1,
-- stores a second best structure
["closest"] = 10,
},
-- verbosity level:
-- 0: no progress output (not recommended)
-- 1: prints the segment whose sidechain is being flipped
-- 2: also prints the rotamer snap indices
["verbosity"] = 2,
}
defaults = {
-- Mode 1 (fast): Flip only Leucine and Isoleucine sidechains. Flipping these tends to be the most productive, although Glutamate is a close third.
{
["xle_only"] = true,
["max_attempts"] = 3,
["min_shake_gain"] = 0.5,
["max_shake_loss"] = 3000,
["skip_freeze_shake_unfreeze"] = false,
},
-- Mode 2 (default): good trade-off between speed and potential gain
{
["xle_only"] = false,
["max_attempts"] = 3,
["min_shake_gain"] = 0.5,
["max_shake_loss"] = 3000,
["skip_freeze_shake_unfreeze"] = false,
},
-- Mode 3 (extensive): Tries more rotamer positions. Does not require an additional sidechain to be flipped on the first shake.
{
["xle_only"] = false,
["max_attempts"] = 6,
["min_shake_gain"] = -0.5,
["max_shake_loss"] = 100000,
["skip_freeze_shake_unfreeze"] = false,
},
-- Mode 4 (disruptive): Foregoes the freeze-shake-unfreeze. This can result in severe disruption to the protein structure. Not terribly useful except perhaps in the opening.
{
["xle_only"] = false,
["max_attempts"] = 3,
["min_shake_gain"] = -0.5,
["max_shake_loss"] = 100000,
["skip_freeze_shake_unfreeze"] = true,
},
}
function SetOptions( new_options )
for name,value in pairs(new_options) do
options[name] = value
end
end
function Dialog_Options()
-- assemble the option dialog
local ask = dialog.CreateDialog("Options")
ask.level_1 = dialog.AddButton( "fast", -1 )
ask.level_2 = dialog.AddButton( "default", -2 )
ask.level_3 = dialog.AddButton( "extensive", -3 )
ask.level_4 = dialog.AddButton( "disruptive", -4 )
ask.label_1b = dialog.AddLabel( "Chooses a pseudorandom start segment, " )
ask.label_1c = dialog.AddLabel( "unless you set a valid segment index:" )
ask.start_segment = dialog.AddTextbox( "start segment", "???" )
ask.separator1 = dialog.AddLabel("")
ask.leu_ile_only = dialog.AddCheckbox( "flip only Leucine and Isoleucine sidechains", false )
ask.separator2 = dialog.AddLabel("")
ask.label_3d = dialog.AddLabel( "maximal number of rotamers to try for each sidechain:" )
ask.max_attempts = dialog.AddTextbox( "attempts", "???" )
ask.label_3b = dialog.AddLabel( "Rotamers too close to the current sidechain position" )
ask.label_3c = dialog.AddLabel( "are skipped automatically and not counted." )
ask.separator3 = dialog.AddLabel("")
ask.label_4a = dialog.AddLabel( "If the freeze-shake-unfreeze gains fewer points," )
ask.label_4b = dialog.AddLabel( "the algorithm skips to the next rotamer:" )
ask.min_shake_gain = dialog.AddTextbox( "minimal gain", "???" )
ask.label_5a = dialog.AddLabel( "Same if the score after the freeze-shake-unfreeze" )
ask.label_5b = dialog.AddLabel( "is further away from the best than:" )
ask.max_shake_loss = dialog.AddTextbox( "maximal gap", "???" )
ask.separator5 = dialog.AddLabel("")
ask.post_scramble = dialog.AddCheckbox( "do a brief scramble iff close to the best score", true )
ask.scramble_gap = dialog.AddTextbox( "threshold:", "???" )
ask.separator7 = dialog.AddLabel("")
ask.label_8a = dialog.AddLabel( "The following option can result in severe disruption" )
ask.label_8b = dialog.AddLabel( "to the protein structure:" )
ask.skip_first_shake = dialog.AddCheckbox( "skip freeze-shake-unfreeze", false )
ask.separator8 = dialog.AddLabel("")
ask.label_9 = dialog.AddLabel("Choose a mode to load pre-defined options:")
ask.ok = dialog.AddButton( "START", 1 )
while( true ) do
-- fill in the options
ask.start_segment.value = tostring( options.start_segment )
ask.leu_ile_only.value = options.xle_only
ask.max_attempts.value = tostring( options.max_attempts )
ask.min_shake_gain.value = tostring( options.min_shake_gain )
ask.max_shake_loss.value = tostring( options.max_shake_loss )
ask.post_scramble.value = options.post_scramble
ask.scramble_gap.value = tostring( options.scrambleThreshold )
ask.skip_first_shake.value = options.skip_freeze_shake_unfreeze
-- show dialog
response = dialog.Show(ask)
print( "Setting options...\n" )
SetOptions({
["start_segment"] = tonumber( ask.start_segment.value ),
["max_attempts"] = tonumber( ask.max_attempts.value ),
["min_shake_gain"] = tonumber( ask.min_shake_gain.value ),
["max_shake_loss"] = tonumber( ask.max_shake_loss.value ),
["scrambleThreshold"] = tonumber( ask.scramble_gap.value ),
["xle_only"] = ask.leu_ile_only.value,
["skip_freeze_shake_unfreeze"] = ask.skip_first_shake.value,
["post_scramble"] = ask.post_scramble.value,
})
if( response == 1 ) then
-- start recipe
return true
elseif( response < 0 ) then
-- set mode
SetOptions( defaults[ -response ] )
else
-- cancel recipe
return false
end -- if
end -- while
end
function get_aa3 ( i )
k = structure.GetAminoAcid ( i )
for j = 1, 20 do
if ( k == single_letter_codes [ j ] ) then
return three_letter_codes [ j ]
end
end
return "Unk"
end
--
-- returns a positive integer which has no common prime factor with the given integer
-- more mathematically: for each integer n: Coprime(n) > 0 such that gcd(n,Coprime(n)) == 1
-- something of a hack, but it seems unlikely we'll ever have puzzles with more than 1146 residues
--
function Coprime( n )
if( n == 0 ) then
-- every integer is a factor of 0
-- but we only need gcd(0,Coprime(0)) == 1
return 1
elseif( n < 0 ) then
n = -n
end
if( n%31 ~= 0 ) then
return 31
elseif( n%37 ~= 0 ) then
return 37
else
return 1
end
end
function pseudorandom ( score )
-- returns a pseudorandom number >= 0 and < 1
-- based on the fractional part of the score
if ( score >= 0 ) then
return score % 1
else
return (-score ) % 1
end
end
function nudge( score )
selection.SelectAll ()
recentbest.Save ()
behavior.SetClashImportance ( 0.2 + 0.2 * pseudorandom ( score ) )
structure.WiggleAll ( 1 )
behavior.SetClashImportance ( 1 )
structure.WiggleAll ( 8 )
recentbest.Restore ()
end
function quick_nudge( score )
selection.SelectAll ()
recentbest.Save ()
behavior.SetClashImportance ( 0.1 + 0.1 * pseudorandom ( score ) )
structure.ShakeSidechainsAll ( 1 )
behavior.SetClashImportance ( 1 )
structure.WiggleAll ( 8 )
structure.ShakeSidechainsAll ( 1 )
structure.WiggleAll (1,false,true)
recentbest.Restore()
end
function FlipAllSidechains()
local segment_count = structure.GetCount()
if( nil == start_segment or start_segment < 1 or start_segment > segment_count ) then
start_segment = math.floor( 1 + segment_count * pseudorandom( best_score ) )
end
-- go through the protein in a less than regular way
increment = Coprime( segment_count )
while( increment >= segment_count ) do
increment = increment - segment_count
end
segment = start_segment
for i = 1,segment_count do
-- flip the sidechain
amino_acid = get_aa3( segment )
if( options.verbosity >= 1 ) then
print( i, "/", segment_count .. ":", amino_acid, segment )
end
if( options.xle_only == false or amino_acid == "Leu" or amino_acid == "Ile" ) then
FlipSidechain( segment )
end
-- move on to next segment
segment = segment + increment
if( segment > segment_count ) then
segment = segment - segment_count
end
end
end
function FlipSidechain( segment )
r = 1
improvement_made = false
-- Lysines and arginines can have 40 or rotamers and its wasteful to go through them all: flipping the
-- sidechains of those amino acids rarely works anyway
-- Till there's something in the API that allows identification of the residue, restrict the number of attempts
n_attempts = 0
while ( r <= rotamer.GetCount(segment) and improvement_made == false and n_attempts < options.max_attempts ) do
rotamer.SetRotamer( segment, r )
score_after_snap = current.GetEnergyScore()
snap_loss = best_score - score_after_snap
if( snap_loss < 0 ) then
-- can conceivably happen
save.Quicksave( options.quicksaves.best )
best_score = score_after_snap
improvement_made = true
print ( "Immediate improvement at segment", segment )
print ( "best score:", best_score )
n_attempts = n_attempts + 1
elseif( snap_loss > options.min_snap_loss ) then
-- if snapping changes it by less than 100, this technique is unlikely to produce a gain
n_attempts = n_attempts + 1
if( options.verbosity >= 2 ) then
print( "", "attempt " .. n_attempts .. ":", "rotamer", r, "/", rotamer.GetCount( segment ) )
end
-- freeze-shake-unfreeze --
gain_from_shake = 0
if( not options.skip_freeze_shake_unfreeze ) then
freeze.Freeze( segment, false, true )
structure.ShakeSidechainsAll( 1 )
freeze.Unfreeze( segment, false, true )
score_after_shake = current.GetEnergyScore()
gain_from_shake = score_after_shake - score_after_snap
else
score_after_shake = score_after_snap
gain_from_shake = 0
end
if( gain_from_shake >= options.min_shake_gain and best_score - score_after_shake < options.max_shake_loss ) then
-- wiggle-shake --
structure.WiggleAll( 15 )
score_after_wiggle = current.GetEnergyScore()
structure.ShakeSidechainsAll( 1 )
score_after_second_shake = current.GetEnergyScore()
if( score_after_second_shake - score_after_wiggle > 0.5 ) then
-- wiggle --
structure.WiggleAll( 10 )
score_after_second_wiggle = current.GetEnergyScore()
if( score_after_second_wiggle - score_after_second_shake < 1 and
best_score - score_after_second_wiggle < 30 ) then
nudge( score_after_second_wiggle )
score = current.GetEnergyScore()
else
score = score_after_second_wiggle
end
else
score = score_after_second_shake
end
if( best_score - score < options.scrambleThreshold and options.post_scramble ) then
quick_nudge( score )
score = current.GetEnergyScore()
end
if( score > best_score ) then
structure.ShakeSidechainsAll( 1 )
save.Quicksave( options.quicksaves.best )
best_score = current.GetEnergyScore()
--print( "best score ", best_score , " snap loss " , snap_loss , " shake gain " , gain_from_shake )
print( "new best score:", best_score )
improvement_made = true;
end
if( improvement_made == false ) then -- might also be an immediate improvement from above -- Wipf 2019-Mar-10
if( score > score_closest ) then
save.Quicksave( options.quicksaves.closest )
score_closest = score
end
end
end -- if shake improved things
end -- if not current rotamer position
r = r + 1
save.Quickload( options.quicksaves.best )
end -- loop over all rotamers
end
---- main script ----
-- tidy up the protein
freeze.UnfreezeAll()
band.DisableAll()
behaviour.SetClashImportance( 1 )
-- initialize
SetOptions( defaults[2] ) -- mode 2 is default mode
save.Quicksave( options.quicksaves.best )
best_score = current.GetEnergyScore()
score_closest = 0
print( "best score:", best_score )
-- gives the user the choice to cancel the recipe here
if( Dialog_Options() ) then
FlipAllSidechains()
end
