Recipe: Quickfix 2.0

created by spvincent

Profile

• Name: Quickfix 2.0
• ID: 47725
• Shared with: Public
• Parent: None
• Children: None
• Created on: January 11, 2014 at 03:20 AM UTC
• Updated on: January 11, 2014 at 03:20 AM UTC
• Description:

  See first comment

Code

-- Constants

kIterationsBetweenShakesPass2 = 3
kIterationsBetweenShakesPass3 = 4
kGainThresholdPass3 = 100

kShakeDistanceThreshold = 14

-- Globals

p2gains = {}
residue_scores = {}
residue_sequence_scores = {}
backbone_scores = {}
ids = {}
secondary_structures = {}

score = 0
best_score = 0

pass_1 = true
pass_1_backbone_threshold = -2000

pass_2 = true
pass_2_threshold = 1000

pass_3 = true

leave_sheets_alone = false
mutate_not_shake = false

-- Quicksave slots

kOriginalStructureOrNewBest = 1	  -- the starting structure, or any subsequent improvement, will be stored in this quicksave slot


function r3 ( i )

-- printing convenience

   return i - i % 0.001

end

function IsPuzzleMutable ()

    for i = 1 , n_residues do
       if ( structure.IsMutable ( i ) ) then
           return true
       end
    end

    return false

end


function GetScore ()

    score = current.GetEnergyScore () 

    if ( score < -999998 ) then

        score = 8000

        for i = 1 , n_residues do
            score = score + current.GetSegmentEnergyScore ( i   )
       end

    end

    return score

end

function GetScores ()

    for i = 1 , n_residues do
         residue_scores [ i ] = current.GetSegmentEnergyScore ( i   )
   end

end

function GetSequenceScores ( local_wiggle_length )

    for i = 1 , n_residues - local_wiggle_length + 1 do
        
        residue_sequence_scores [ i ] = 0
        
        sheet_found = false
        for j = i , i + local_wiggle_length - 1 do
          
            residue_sequence_scores [ i ] = residue_sequence_scores [ i ] + residue_scores [ j ]
            if ( secondary_structures [ j ] == "E"  ) then
                sheet_found = true
            end

        end
    
       if ( ( sheet_found == true ) and ( leave_sheets_alone == true ) ) then
           residue_sequence_scores [ i ] = residue_sequence_scores [ i ] + 999999
       end
 
    end

end

function  ShellSort ( ids , sequence_scores , n )

    -- Adapted from Numerical Recipes in C

    local  inc = 1

    repeat
      inc = inc * 3 + 1
    until inc > n

    repeat

        inc = inc / 3
        inc = inc - inc % 1

        for i = inc + 1 , n do

            v = sequence_scores [ i ]
            w = ids [  i ]

            j = i
            flag = false

            while ( flag == false and sequence_scores [ j - inc ]  > v ) do

                sequence_scores [ j ] = sequence_scores [ j - inc ]
                ids [ j ] = ids [ j - inc ]

                j = j - inc

                if ( j <= inc ) then
                  flag = true
                end

            end

          sequence_scores [ j ] = v
           ids [ j ] = w

        end

    until inc <= 1

end

function rebuild ( start_idx , end_idx )

    selection.DeselectAll ( )
    selection.SelectRange ( start_idx , end_idx )
    -- structure.SetSecondaryStructureSelected ( "L" )

    -- Rebuild

    structure.RebuildSelected ( 1 )
    recentbest.Save ()

    structure.RebuildSelected ( 3 )
    recentbest.Restore ()

    -- Wiggle sidechains

    selection.DeselectAll ( )

    mid_x = ( start_idx + end_idx ) / 2

    for i = 1 , n_residues do
        d = structure.GetDistance ( i , mid_x )
        if ( d < kShakeDistanceThreshold ) then
           selection.Select ( i )
        end
    end

    structure.WiggleAll ( 1, false , true )
 
   -- Wiggle

    recentbest.Save ()
    selection.DeselectAll ( )
    selection.SelectRange ( start_idx , end_idx )

    structure.LocalWiggleSelected ( 12 )
    recentbest.Restore ()


end

function record_improvement ()

    if ( score > best_score ) then
        gain = score - best_score
        best_score = score
        save.LoadSecondaryStructure () 
        save.Quicksave ( kOriginalStructureOrNewBest )
        print ( "  Improvement to " .. r3 ( score ) )
      else
        gain = 0
    end

    return gain

end

function Pass1 ()

    -- Pass 1

    -- Look for consecutive segments each of whose backbone score is less than pass_1_backbone_threshold.
    -- Default value is such that this will address mostly cutpoint closure issues.
    
    print ( "Pass 1" )

    for i = 1 , n_residues do
        backbone_scores [ i ] = current.GetSegmentEnergySubscore ( i , "backbone" )
    end

    for i = 2 , n_residues do
        if ( ( backbone_scores [ i - 1 ] < pass_1_backbone_threshold ) and 
             ( backbone_scores [ i ] < pass_1_backbone_threshold ) ) then
           
            save.Quickload ( kOriginalStructureOrNewBest )
   
            selection.DeselectAll ( )
            selection.SelectRange(  i - 1 , i )
          
            print ( "lw " .. i - 1 .. "-" .. i )

            structure.LocalWiggleSelected ( 2 )

            score = GetScore ()
        
            record_improvement ()
        end
    end

end

function Pass2 ()

    -- Local wiggle with an occasional global sidechain wiggle

    -- Make sure we don't hammer precisely the same region each time

    prev_idx = -1
    prev_prev_idx = -1
    local   local_wiggle_length = 3
    
    print ( "" )
    print ( "Pass 2" )

    ns = 0

    repeat

         if ( ns % kIterationsBetweenShakesPass2 == 0 ) then

             print ( "ws" )

             save.Quickload ( kOriginalStructureOrNewBest )

             selection.SelectAll ( )
             structure.WiggleAll ( 1, false , true )

             score = GetScore ()
            
             record_improvement ()
 
        end

        ns = ns + 1

        save.Quickload ( kOriginalStructureOrNewBest )

        GetScores ()
        GetSequenceScores ( local_wiggle_length )

        for i = 1 , n_residues - local_wiggle_length + 1 do
           ids [ i ] = i
        end

        ShellSort ( ids , residue_sequence_scores , n_residues - local_wiggle_length + 1 )

        start_idx = ids [ 1 ]
        
        if ( ( start_idx == prev_idx ) or ( start_idx == prev_prev_idx ) ) then
           start_idx = ids [ 2 ]
   
           if ( ( start_idx == prev_idx ) or ( start_idx == prev_prev_idx ) ) then
              start_idx = ids [ 3 ]
           end
        end
       
        end_idx = start_idx + local_wiggle_length - 1

        prev_prev_idx = prev_idx
        prev_idx = start_idx

        print ( "lw " .. start_idx .. "-" .. end_idx )

        selection.DeselectAll ( )
        selection.SelectRange ( start_idx , end_idx )
          
        structure.LocalWiggleSelected ( 8 )

        score = GetScore ()
        
        gain = record_improvement ()

         if ( gain < 1 ) then

            save.Quickload ( kOriginalStructureOrNewBest )

            print ( "rb " .. start_idx .. "-" .. end_idx )

            rebuild ( start_idx , end_idx )

            score = GetScore ()
            
            gain = record_improvement ()
 
        end

        p2gains [ ns ] = gain

    until ( ns > 2 and ( p2gains [ ns ] + p2gains [ ns - 1 ] + p2gains [ ns - 2 ]  < pass_2_threshold ) )
  
end

function Pass3 ()

    print ( "" )
    print ( "Pass 3" )

    ns = 0
    local   local_wiggle_length = 4

    while true do

        save.Quickload ( kOriginalStructureOrNewBest )

        GetScores ()
        GetSequenceScores ( local_wiggle_length )

        for i = 1 , n_residues - local_wiggle_length + 1 do
           ids [ i ] = i
        end

        ShellSort ( ids , residue_sequence_scores , n_residues - local_wiggle_length + 1 )

        start_idx = ids [ 1 ]
        
        if ( ( start_idx == prev_idx ) or ( start_idx == prev_prev_idx ) ) then
           start_idx = ids [ 2 ]
   
           if ( ( start_idx == prev_idx ) or ( start_idx == prev_prev_idx ) ) then
              start_idx = ids [ 3 ]
           end
        end
       
        end_idx = start_idx + local_wiggle_length - 1

        prev_prev_idx = prev_idx
        prev_idx = start_idx

        print ( "lw " .. start_idx .. "-" .. end_idx )

        selection.DeselectAll ( )
        selection.SelectRange( start_idx , end_idx )
          
        structure.LocalWiggleSelected ( 9 )

        score = GetScore ()
            
        gain = record_improvement ()
 
        if ( gain < kGainThresholdPass3 ) then

            save.Quickload ( kOriginalStructureOrNewBest )

            print ( "rb " .. start_idx .. "-" .. end_idx )

            rebuild ( start_idx , end_idx )

            score = GetScore ()
            
            record_improvement ()
 
        end

        if ( ns % kIterationsBetweenShakesPass3 == 0 ) then

             save.Quickload ( kOriginalStructureOrNewBest )

             selection.SelectAll ( )
             
             if ( mutate_not_shake == true ) then
                  print ( "m" )
                  structure.MutateSidechainsSelected ( 1 )               
               else
                 print ( "sa" )
                 structure.ShakeSidechainsSelected ( 1 )
             end

             score = GetScore ()
            
             record_improvement ()
 
        end

        ns = ns + 1

     end        -- while true

end

function get_parameters ()

    local dlog = dialog.CreateDialog ( "Quickfix 2.0" )

    dlog.pass_1= dialog.AddCheckbox ( "Pass 1" , pass_1 )  

    dlog.pass_2 = dialog.AddCheckbox ( "Pass 2" , pass_2 )  
    
    dlog.pass_3 = dialog.AddCheckbox ( "Pass 3" , pass_3 )  

    dlog.leave_sheets_alone = dialog.AddCheckbox ( "Leave sheets alone" , leave_sheets_alone )  

    design_puzzle = IsPuzzleMutable ()

    if ( design_puzzle == true ) then
        dlog.mutate= dialog.AddCheckbox ( "Mutate not shake" , false )  
    end

    dlog.ok = dialog.AddButton ( "OK" , 1 )
    dlog.cancel = dialog.AddButton ( "Cancel" , 0 )

    if ( dialog.Show ( dlog ) > 0 ) then
    
        pass_1 = dlog.pass_1.value

        pass_2 = dlog.pass_2.value

        leave_sheets_alone = dlog.leave_sheets_alone.value

        if ( design_puzzle == true ) then
           mutate_not_shake = dlog.mutate.value
        end

        return true

      else
        return false
    end

end


function main ()

    print (  "Quickfix 2.0" )
    print (  "" )

    band.DisableAll ()
    freeze.UnfreezeAll ()

    n_residues = structure.GetCount ()

    save.SaveSecondaryStructure () 

    behavior.SetClashImportance ( 1 )

    best_score = GetScore ()
    print (  "Start score : " .. r3  ( best_score )  )

    save.Quicksave ( kOriginalStructureOrNewBest )

    for i = 1 , n_residues do
        secondary_structures [ i ] = structure.GetSecondaryStructure ( i )
    end

    if ( get_parameters () == false ) then
        error ()
    end

    print ( "lw = Local Wiggle" )
    print ( "rb = Rebuild" )
    print ( "ws = Wiggle Sidechains" )
    if ( mutate_not_shake == true ) then
       print ( "m = Mutate All" )
      else
        print ( "sa = Shake All" )
    end
    print ( "" )

   if ( pass_1 == true ) then
        Pass1 ()
   end

   if ( pass_2 == true ) then
        Pass2 ()
   end

   if ( pass_3 == true ) then
        Pass3 ()
   end

      
end

function cleanup ()
   
    print ( "Cleaning up" )
    save.LoadSecondaryStructure () 
    save.Quickload ( kOriginalStructureOrNewBest )

    band.EnableAll ()

end

--main ()
xpcall ( main , cleanup )

Comments

[spvincent]

Quickfix is not intended to improve the score, although it might have that happy side effect. Rather; the script is designed to rapidly recover from extremely low scores such as those that occur after actions like cut point closure and register shifts, while making no more than minimal changes to the protein. Eventually the score should rise to a value where it's safe to kill the script and subsequent wiggling won't cause the protein to fly apart.

The default values work pretty well: still some notes may be helpful.

Pass 1

Looks for adjacent residues which have a combined backbone score less than a threshold value (default -2000): then applies a quick local wiggle. Particularly intended to help tidy up after cut point closure. Very quick so might as well always leave it on.

At the start of the script, anything frozen will be unfrozen and all bands will be disabled. On killing the script bands will be restored and the pose set to that corresponding to the most recent improvement.

Pass 2

First, wiggles sidechains. Then looks at the total scores of all 3-residue segments; selecting the worst and local wiggling it. If said local wiggle doesn't improve the score, try rebuilding that section a few times: then take the best-scoring result and local wiggle.

Repeat all the above (taking care not to hammer away unsuccessfully on one section should that prove intractable), until the total improvement for three consecutive tries is less than some threshold (default=1000), whereupon things move on to Pass 3 (if selected).

Pass 3

Similar to Pass 2 except it local wiggles/rebuilds the worst-scoring 4-residue section (as opposed to 3). Rebuilds if the gain from local wiggle is less than a certain threshold (default=100) Also at intervals does a Shake (or Mutate: see below). Does not terminate: it's up to the user to kill the script, whereupon the highest scoring solution made is restored.

Leave sheets alone

In Passes 2 and 3, don't work on sections that include a sheet (so they don't get mangled)

Mutate not shake

For design puzzles only. If checked, substitutes a mutate for a shake in Pass 3. Often better but Mutate can be quite slow.

[Bruno Kestemont]

Very useful after hand folding !

[AsDawnBreaks]

Yes, this will be helpful. So long taking eons trying to recover from a very low score. :)