Recipe: Quickfix 1.5

created by spvincent

Profile

• Name: Quickfix 1.5
• ID: 44443
• Shared with: Public
• Parent: None
• Children:
  • Quickfix 1.5.2
  • Quickfix 1.5 save bands
  • Quickfix 1.5 With Mutate
• Created on: October 30, 2012 at 03:56 AM UTC
• Updated on: October 30, 2012 at 03:56 AM UTC
• Description:

  See first comment

Code

-- Constants

kIBackboneThresholdPass1 = -2000
kIterationsBetweenShakesPass2 = 3
kIterationsBetweenShakesPass3 = 4
kGainThresholdPass3 = 100

kShakeDistanceThreshold = 14

-- Globals

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

score = 0
best_score = 0
local_wiggle_length = 4


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

    for i = 1 , n_residues - local_wiggle_length + 1 do
        
        residue_sequence_scores [ i ] = 0
        
        for j = i , i + local_wiggle_length - 1 do
            residue_sequence_scores [ i ] = residue_sequence_scores [ i ] + residue_scores [ j ]
        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 main ()

    print (  "Quickfix 1.5" )
    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 )

    print ( "lw = Local Wiggle" )
    print ( "rb = Rebuild" )
    print ( "ws = Wiggle Sidechains" )
    print ( "sa = Shake All" )
    print ( "" )

    -- Pass 1

    -- Look for consecutive segments whose backbone score is each less than -1000, such
    -- as arise from cut point closure.

    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 ] < kIBackboneThresholdPass1 ) and 
              ( backbone_scores [ i ] < kIBackboneThresholdPass1 ) ) 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

    -- Pass 2

    -- 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_wiggle_length = 3

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

        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 ]  < 5000 ) )
        
     -- Pass 3

    ns = 0
    local_wiggle_length = 4

    while true do

        save.Quickload ( kOriginalStructureOrNewBest )

        GetScores ()
        GetSequenceScores ()

        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

        ns = ns + 1

        if ( ns % kIterationsBetweenShakesPass3 == 0 ) then

             print ( "sa" )

             save.Quickload ( kOriginalStructureOrNewBest )

             selection.SelectAll ( )
             structure.ShakeSidechainsSelected ( 1 )

             score = GetScore ()
            
             record_improvement ()
 
        end

     end        -- while true

end

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

    band.EnableAll ()

end

--main ()
xpcall ( main , cleanup )

Comments

[spvincent]

A new and hopefully significantly improved version of Quickfix.

The script is designed to rapidly recover from extremely low scores such as those that occur after actions like cutpoint closure and register shift, 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 using wiggle won't cause the protein to fly apart.