Profile
- Name
- V2 Sidechain Flipper 2.2
- ID
- 31000
- Shared with
- Public
- Parent
- None
- Children
- Created on
- August 11, 2011 at 18:11 PM UTC
- Updated on
- August 11, 2011 at 18:11 PM UTC
- Description
A working V2 version
Best for
Code
start_idx = 0 -- if this is a residue id, use it as the start point. If not
-- use a semi-random number computed from the score
start_from_best = false -- set to false if you want to start from the current configuration, not the best.
-- For each rotamer postion of each side chain, does a flip, freeze, shake, unfreeze, wiggle, shake, wiggle.
-- Most generally, if the first shake fails to produce an improvement in the score, stops working on that rotamer position and proceeds to the next.
level = 2
-- level 1. Flip Leu and Ile only. Flipping these tend to be the most productive, although Glu is a close third.
-- level 2. Default. Best trade-off between potential gain and speed.
-- level 3. Tries more positions. Unlike the lower 2 levels, does not require that on the first shake an additional side chain be flipped.
-- level 4. Foregoes the freeze-shake-unfreeze bit of the algorithm entirely. This can result in severe disruption to the protein structure.
-- Not terribly useful except perhaps in the opening.
post_scramble = 2
-- post_scramble = 1. No sidechain scramble. Won't leak walk juice.
-- post_scramble = 2. Does a brief scramble iff a flip gets within two (see kScrambleThreshold) points of the original.
kScrambleThreshold = 2
-- 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"}
-- Use of save.Quicksave slots: change if desired
kBest = 1 -- Used to hold the initial structure and any improved structures
kClosest = 10 -- Holds the best scoring structure that does not exceed that in kBest
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
function Coprime ( n )
-- returns a number that has no common factor with n. Something of a hack but
-- it seems unlikely we'll ever have puzzles with less than 31 or more than 1146 residues
if ( n < 31 ) then
return 1
elseif ( n >= 31*37 ) then
return 1
elseif ( n%31 ~= 0 ) then
return 31
else
return 37
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 Go ( n )
inc = Coprime ( n )
idx = start_idx
-- use this to go through the protein in a less than regular way
for i = 1 , n do
n_rotamers = rotamer.GetCount ( idx )
aa = get_aa3 ( idx )
if ( leu_ile_only == false or aa == "Leu" or aa == "Ile" ) then
print ( aa , " " , idx , " (", i, "/", n,")", " n_rotamers " , n_rotamers )
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 <= n_rotamers and improvement_made == false and n_attempts < max_attempts ) do
save.Quickload ( kBest )
selection.SelectAll ()
rotamer.SetRotamer ( idx , r )
score_after_snap = current.GetEnergyScore()
--print ( "rotamer ", r , " n_attempts ", n_attempts )
snap_loss = best_score - score_after_snap
if ( snap_loss < 0 ) then
-- Can conceivably happen
save.Quicksave ( kBest )
best_score = score_after_snap
improvement_made = true
print ( "Immediate improvement at ", idx )
print ( "best score ", best_score )
n_attempts = n_attempts + 1
elseif ( snap_loss > min_snap_loss ) then
-- If snapping changes it by less than 100, either we are already at the current rotamer position
-- or this technique is unlikely to produce a gain
n_attempts = n_attempts + 1
gain_from_shake = 0
if ( level < 4 ) then
selection.DeselectAll ()
selection.Select ( idx )
freeze.FreezeSelected ( false , true )
selection.SelectAll ()
structure.ShakeSidechainsAll ( 1 )
score_after_shake = current.GetEnergyScore ()
gain_from_shake = score_after_shake - score_after_snap
freeze.UnfreezeAll ()
else
score_after_shake = score_after_snap
gain_from_shake = 0
end
if ( gain_from_shake > lowest_allowable_shake_gain and best_score - score_after_shake < max_shake_loss ) then
--print ( "score after shake " , best_score - score_after_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
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 < kScrambleThreshold and post_scramble == 2 ) then
quick_nudge ( score )
score = current.GetEnergyScore ()
end
if ( score > best_score ) then
structure.ShakeSidechainsAll ( 1 )
save.Quicksave ( kBest )
best_score = current.GetEnergyScore ()
print ( "Improvement at " , idx )
--print ( "best score ", best_score , " snap loss " , snap_loss , " shake gain " , gain_from_shake )
print ( "best_score ", best_score )
improvement_made = true;
end
if ( improvement_made == false ) then
if ( score > score_in_qs10 ) then
save.Quicksave ( kClosest )
score_in_qs10 = score
--print ( "Save.Quicksave 10 " , score )
end
end
end -- if shake improved things
end -- if not current rotamer position
r = r + 1
end
--if ( n_attempts > max_attempts and improvement_made == false ) then
--print ( "R " , r )
--end
end
idx = idx + inc
if ( idx > n ) then
idx = idx - n
end
end -- for i
end
n_residues = structure.GetCount ()
behavior.SetClashImportance ( 1 )
if ( start_from_best == true ) then
absolutebest.Restore ()
end
-- Tidy up the protein
freeze.UnfreezeAll ()
band_count = band.GetCount ()
for i = band_count , 1 , -1 do
band.Disable ( i )
end
save.Quicksave ( kBest )
best_score = current.GetEnergyScore ()
score_in_qs10 = 0
print ( "best score ", best_score )
if ( level == 1 ) then
leu_ile_only = true
else
leu_ile_only = false
end
if ( level == 1 or level == 2 ) then
lowest_allowable_shake_gain = 0.5
else
lowest_allowable_shake_gain = -0.5
end
if ( level == 1 or level == 2 or level == 4 ) then
max_attempts = 3
else
max_attempts = 6
end
min_snap_loss = 60
if ( level == 1 or level == 2 ) then
max_shake_loss = 3000
else
max_shake_loss = 100000
end
if ( start_idx < 1 or start_idx > n_residues ) then
r = pseudorandom ( best_score )
start_idx = 1 + n_residues * r
start_idx = start_idx - start_idx % 1
end
Go ( n_residues )