Recipe: Hiding Rebuild V2

created by spvincent

Profile

• Name: Hiding Rebuild V2
• ID: 37872
• Shared with: Public
• Parent: None
• Children:
  • My_541_Hiding Rebuild V2
  • Hiding Rebuild V2 MOD
  • Bonding Rebuild V2
  • Clashing Rebuild V2
  • Packing Rebuild V2
  • Rebuild Hide Phobes
  • Packing Rebuild V2
  • Hiding Rebuild V2.1
  • Hiding Remix V2
• Created on: January 24, 2012 at 23:05 PM UTC
• Updated on: January 24, 2012 at 23:05 PM UTC
• Description:

  Optimizes for a specific score component (default is hiding) . See first comment for more details.

Code

max_rebuild_length =  6                    -- Maximum rebuild length.          
min_rebuild_length =  6                     --  Must be less than or equal to max_rebuild_length. Suggest not going below 5.

number_of_rebuilds = 5                      -- It'll look for the  best score from this number of rebuilds.
 
component_to_optimize = "hiding"    -- Change to something else if you like

hours_to_run = 1                                    -- Change to taste

ss_control = 3                                 

        -- ss_control == 1                    Only rebuilds loops
        -- ss_control == 2                    Rebuilds everything. Helices and sheets will be temporarily converted to loops
        -- ss_control == 3                    Rebuild loops but allows a single non-loop residue in the rebuild segment, which will be temporarily converted to loops

min_residue = 1                                       -- Change these to rebuild part of a protein only. For example, to 
                                                                    --  rebuild  53-64, set min_residue to 53 and max_residue to 64. To restore, set to 
                                                                    -- 1 and 999  (or any large number) respectively
max_residue =  999

kWorstAllowableAngle = 82                  -- don't allow backbone angles after rebuild to be less than this (degrees). The idea is that
                                                                    -- severely kinked structures aren't worth the trouble of post processing, Note: the angle calculated
                                                                    -- is that between amino acids (their alpha carbons one assumes): in reality there are 3 bonds between each residue.

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

        -- Use of quicksave slots

kOriginalStructureOrNewBest = 1
kBestPostRebuild = 2
kTemp = 3

original_secondary_structure = {}

n_residues = 0
best_score = 0

start_time = 0
time_exceeded = false

cos_worst_allowable_angle = 0

function GetScore ()

    if ( component_to_optimize == ""  ) then
      score = current.GetEnergyScore () 
    else
      score = 0
      for i = 1 , n_residues do
         score = score + current.GetSegmentEnergySubscore ( i , component_to_optimize  )
      end
    end

    return score

end

function Coprime ( n )

-- returns a number that has no common factors with n
-- an ugly hack

    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 IsSegmentALoop ( first , last )

  -- Returns number of non-loop  residues (helices or sheets)

   n_non_loops = 0
   for i = first , last do
       if ( structure.GetSecondaryStructure ( i ) ~= "L" ) then
          n_non_loops = n_non_loops + 1
        end
     end

   return  n_non_loops

end

function ConvertToLoop ( first , last )

    for  k = first , last do
        if ( original_secondary_structure [ k ] ~= "L" ) then
            selection.DeselectAll ( )
            selection.Select ( k )
            structure.SetSecondaryStructureSelected ( "L" )
        end
    end
end

function RestoreStructure ( first , last )

    for k = first, last do
        selection.DeselectAll ()
        selection.Select ( k )
        structure.SetSecondaryStructureSelected  ( original_secondary_structure [ k ] )
    end
end

function roundto3 ( i )

-- printing convenience

   return i - i % 0.001

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 ()

    score_in_kTemp = GetScore()
    save.Quicksave ( kTemp)

    behavior.SetClashImportance ( 0.4  + 0.2 * pseudorandom ( score ) )
    structure.WiggleSelected ( 1 )
    behavior.SetClashImportance ( 1 )
    structure.WiggleSelected ( 6 )

    if ( GetScore() < score_in_kTemp ) then
      save.Quickload ( kTemp )
    end

end

function  WiggleShakeWiggleNudge ()

     -- returns the score after doing a WSW:  short-circuiting if things don't look promising

    score_a = GetScore ()

    structure.WiggleSelected ( 10 )
    score_after_wiggle = GetScore ()

    if ( score_after_wiggle - score_a  < 10 and best_score - score_after_wiggle > 30 ) then
        --    The wiggle got stuck and didn't achieve anything
        nudge ( score_after_wiggle )
        score_after_wiggle = GetScore ()
    end
    
    structure.ShakeSidechainsSelected ( 1 )
    
    score_after_second_shake = GetScore ()

    structure.WiggleSelected ( 10 )
    score_after_second_wiggle = GetScore ()

    nudge ( score_after_second_wiggle )

    score = GetScore  ()

    return  score

end


function post_rebuild ()

  improvement_made = false

  save.Quickload ( kBestPostRebuild )

  WiggleShakeWiggleNudge ()

  curr_score = GetScore ()

  if ( curr_score > best_score ) then

       best_score = curr_score

       improvement_made = true
       save.LoadSecondaryStructure ()

       save.Quicksave ( kOriginalStructureOrNewBest )
       print ( "Improvement to "..  roundto3 ( best_score ) .. " real score : " .. roundto3 ( current.GetEnergyScore () ) )
                     
    end

end

function get_cos_backbone_angle ( i )

    if ( i == 1 or i == n_residues ) then
      return -1
     else
  
      -- Cosine rule

       a = structure.GetDistance ( i  - 1 , i )
       b = structure.GetDistance ( i , i + 1 )
       c = structure.GetDistance ( i - 1 , i + 1 )

       cos_angle = ( a*a + b*b - c*c ) /  ( 2*a*b )

       return cos_angle

    end

end

function is_backbone_reasonably_angled ( start_idx , end_idx )

  for i = start_idx , end_idx do
    if ( get_cos_backbone_angle ( i ) > cos_worst_allowable_angle ) then
      -- print (  math.deg ( math.acos ( cos_angle ) ) )
      return false
    end
  end

  return true

end

function rebuild_seg ( start_idx , end_idx )

--    return codes
--
--    0	    rebuild had no effect: seg probably locked
--    1    normal

  best_score_after_rebuild_and_shake = -99999
  good_structure_found = false

   for i = 1 , number_of_rebuilds do

      save.Quickload (  kOriginalStructureOrNewBest )

      ConvertToLoop ( start_idx ,  end_idx )
  
      selection.DeselectAll ()
      selection.SelectRange (  start_idx  ,  end_idx  )
     
      -- for some reason, structure.RebuildSelected ( 1)  often leaves  the protein unchanged

      k  = 0 
      repeat 
          k = k + 1
         structure.RebuildSelected ( k )
         score_after_rebuild = GetScore ()
      until  score_after_rebuild ~= best_score or k  > 2
  
      if ( math.abs ( score_after_rebuild - best_score ) < 1e-5 ) then
          -- Probably locked
          return 0
      elseif ( is_backbone_reasonably_angled ( start_idx , end_idx ) ) then
          good_structure_found = true
          selection.SelectAll ()
          structure.ShakeSidechainsSelected ( 1 )
          score_after_shake = GetScore()

           if ( score_after_shake > best_score_after_rebuild_and_shake ) then
              best_score_after_rebuild_and_shake = score_after_shake
              save.Quicksave ( kBestPostRebuild )
           end
       end     
    
  end     -- for i

  if ( good_structure_found == true ) then
    return 1
   else
    return 0
  end

end

function rebuild ( rebuild_length  )

    n_possible_segs = n_residues - rebuild_length + 1

    inc = Coprime ( n_possible_segs  )

    --    The point of  this Coprime, start_idx business is so we go through the protein in  a
    --    nonuniform way. Rather than rebuilding 1-4, 2-5 , 3-6 for example, we might do 21 4, 53-57, 3-6  
    --    or something like  that

    for  i = 1 , n_possible_segs do
         
        improvement_made = false

        end_idx = start_idx + rebuild_length - 1

        if (  start_idx >= min_residue and end_idx <= max_residue and
                   ( ss_control == 2 or 
                   ( ss_control == 1 and IsSegmentALoop ( start_idx , end_idx  ) == 0 ) or
                   ( ss_control == 3 and IsSegmentALoop ( start_idx , end_idx  ) < 2  ) ) ) then
    
            rb =  rb + 1
               
            print (  "rb ".. rb .. " " .. start_idx  .."-" ..  end_idx ..  " ("  ..i .. "/" .. n_possible_segs .. ")"  )

           k = rebuild_seg ( start_idx, end_idx )

           if ( k == 1 ) then
             post_rebuild ()
           end

           if ( ss_control > 1 ) then
               save.LoadSecondaryStructure ()
           end

           time = os.time ()
           delta_time = os.difftime ( time , start_time )
           
           if ( delta_time > hours_to_run * 3600 ) then
              time_exceeded = true
              return 1
           end
      end

      start_idx = start_idx + inc
      if ( start_idx >  n_possible_segs ) then
          start_idx =  start_idx - n_possible_segs
      end
        
    end        --    for i

end

-- Start here.
-- Tidy up if necessary

band.DisableAll ()

n_residues = structure.GetCount ()

if ( ss_control  > 1 ) then
    for i = 1, n_residues do
        original_secondary_structure [ i ] = structure.GetSecondaryStructure ( i )
    end
end

save.SaveSecondaryStructure () 

rb = 0
save.Quicksave ( kOriginalStructureOrNewBest )

behavior.SetClashImportance ( 1 )

best_score = GetScore ()
print (  "Start " .. component_to_optimize .. " score : " .. roundto3 ( best_score ) .. " real score : " .. roundto3 ( current.GetEnergyScore () ) )

if ( ss_control == 1 ) then
    print (  "Loops only" )
  elseif ( ss_control == 2 ) then
    print (  "Loops, sheets, and helices" )
 else
   print ( "Loops with one non-loop allowed" )
end

if ( max_residue > n_residues ) then
    max_residue = n_residues
end

if ( min_residue < 1 ) then
    min_residue = 1
end

if ( max_residue <= min_residue ) then
    print ( "Rebuild  range error" )
end

print  ( "Rebuild  range " .. min_residue .. " to " .. max_residue )

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

cos_worst_allowable_angle = math.cos ( math.rad ( kWorstAllowableAngle ) )
start_time = os.time()

repeat
    for rebuild_length = max_rebuild_length , min_rebuild_length , -1  do 

      if ( time_exceeded == false ) then
        rebuild ( rebuild_length ) 
      end
  end 
until  time_exceeded == true

save.Quickload ( kOriginalStructureOrNewBest )

band.EnableAll ()

Comments

[spvincent]

This script is not intended to improve the score, although it may indeed have that happy side effect.

Rather the aim is to improve a specific component of the score like hiding or bonding at the likely expense of the total score: subsequent scripts like Quakes or Rebuilders can then be used to tidy up or you can adjust things manually. The default component to improve is hiding (other styles of script generally seem to have a problem with this component) but this can readily be changed (see component_to_optimize at start). It's probably best to confine the scripts effects to a specific region of the protein like a large loop: see min_residue and max_residue at start.

This script runs for a specific time: the default is one hour. If you stop it before then you may have to step back through the undo graph to find the end product of the script.

It looks for the "best" score after multiple rebuild/shakes: rejecting rebuild results that have unduly kinked backbones. It then takes that structure, wiggles, shakes, wiggles at lower CI, ands wiggles again.

[mottiger]

Am I right that the component option can be changed to clashing, packing, hiding, bonding and sidechain?

And how did you manage that the script runs for 1hr? Is there a way to change that?

[spvincent]

Yes to question 1. Hiding's the default because it's my general impression that that's the component hardest to improve with other types of script.

Q2: Use the os library in Lua. Something like this.

start_time = os.time ()

.. Do something. Call it X.

end_time = os.time ()

time_taken = os.difftime ( end_time , start_time )

print ( "Seconds to do X " .. time_taken * 3600 )