Recipe: Quickfix 3.0

created by spvincent

Profile

• Name: Quickfix 3.0
• ID: 100031
• Shared with: Public
• Parent: None
• Children:
  • Quickfix 3.0freezedbanded
  • Quickfix 3.1.7a
  • Quickfix 3.6.3
  • Quickfix 3.0 BP
• Created on: October 02, 2014 at 02:06 AM UTC
• Updated on: October 02, 2014 at 02:06 AM UTC
• Description:

  See first comment

Code

-- Constants

kWiggleBackbonePairsThreshold = -2000
kIdealityTripletsThreshold_a = -100
kIdealityTripletsThreshold_b = -30

kWiggleBackboneTripletsThreshold = 200
kWiggleBackboneQuadrupletsThreshold = 100
kRebuildQuadrupletsThreshold = 100

kNWorstClashesPt1 = 10

kShakeDistanceThreshold = 12


-- Globals

leave_sheet_backbones_alone = false
leave_helix_backbones_alone = false

score = 0
best_score = 0

residue_scores = {}
backbone_scores = {}
ideality_scores = {}
clashing_scores = {}
sequence_scores = {}
ids = {}
secondary_structures = {}

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 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 IsSequenceValid ( start_idx , end_idx )

    n_illegitimate_residues = 0

    for i = start_idx , end_idx do

        if ( ( ( secondary_structures [ i ] == "E"  ) and ( leave_sheet_backbones_alone == true ) ) or 
             ( ( secondary_structures [ i ] == "H"  ) and ( leave_helix_backbones_alone == true ) ) ) then
            n_illegitimate_residues = n_illegitimate_residues + 1
        end

    end

    if ( n_illegitimate_residues <= 1 ) then
        return true
      else
        return false
    end

end

function GetSequenceScores ( input_scores , length , respect_ss_restrictions )

    idx = 0

    for i = 1 , n_residues - length + 1 do

        if ( ( respect_ss_restrictions == false ) or ( IsSequenceValid ( i , i + length - 1 ) ) ) then

            idx = idx + 1
            total_scores = 0

            for j = i , i + length - 1 do
              total_scores = total_scores + input_scores [ j ]
            end

            sequence_scores [ idx ] = total_scores
          ids [ idx ] = i

        end

    end

    ShellSort ( ids , sequence_scores , idx )

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

    -- Shake or Mutate the n worst-clashing sidechains

    n = math.min ( n , n_residues )

    if ( n == n_residues ) then
        selection.SelectAll ()
      else

        for i = 1 , n_residues do
            clashing_scores [ i ] = current.GetSegmentEnergySubscore ( i , "clashing" )
            ids [ i ] = i
        end

        ShellSort ( ids , clashing_scores , n_residues )

        selection.DeselectAll ()
  
        for i = 1 , n do
            selection.Select ( ids [ i ]  )
        end
    end

    if ( mutate_not_shake == true ) then
         print ( "mu"  )
         structure.MutateSidechainsSelected ( 1 )
      else
         print ( "sh"  )
         structure.ShakeSidechainsSelected ( 1 )
    end

    score = GetScore ()
            
    record_improvement ()

end

function WiggleSidechains ()

    print ( "ws" )

    save.Quickload ( kOriginalStructureOrNewBest )

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

    score = GetScore ()
            
    record_improvement ()
 
end

function WiggleBackbonePairs ()


    -- Look for consecutive segments each of whose backbone score is less than kWiggleBackbonePairsThreshold.
    -- Default value is such that this will address mostly cutpoint closure issues.
  
    save.Quickload ( kOriginalStructureOrNewBest )

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

    GetSequenceScores ( backbone_scores , 2 , false )
   
    for i = 1 , n_residues - 1 do
        if ( sequence_scores [ i ] < kWiggleBackbonePairsThreshold ) then
           
            save.Quickload ( kOriginalStructureOrNewBest )
   
            selection.DeselectAll ( )
            selection.SelectRange ( ids [ i ]  , ids [ i ]  + 1 )
          
            print ( "lw " .. i - 1 .. "-" .. i )

            structure.LocalWiggleSelected ( 2 )

            score = GetScore ()
        
            record_improvement ()

        end
    end

end

function IdealizeBackboneTriplets ( threshold )

    -- Correct really bad Idealities

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

    GetSequenceScores ( ideality_scores , 3 , false )
   
    for i = 1 , n_residues - 2 do
        if ( sequence_scores [ i ] < threshold ) then
           
            save.Quickload ( kOriginalStructureOrNewBest )
   
            start_idx = ids [ i ]
            end_idx = start_idx + 2

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

            if ( start_idx > 1 ) then
                 structure.InsertCut ( start_idx - 1 )
            end

            if ( end_idx < n_residues ) then
                structure.InsertCut ( end_idx - 1 )
            end

            selection.DeselectAll ()
            selection.SelectRange ( start_idx , end_idx  )
        
            structure.IdealizeSelected ()
 
            if ( start_idx > 1 ) then
                 structure.DeleteCut ( start_idx - 1 )
            end

            if ( end_idx < n_residues ) then
                structure.DeleteCut ( end_idx - 1 )
            end

           selection.DeselectAll ()
           selection.SelectRange ( start_idx , end_idx )
          
           structure.LocalWiggleSelected ( 2 )
        
           score = GetScore ()

           record_improvement ()
 
        end
    end
          
end

function WiggleBackboneTriplets ()

    local  p2gains = {}

    prev_idx = -1
    prev_prev_idx = -1
    ns = 0

    repeat

        ns = ns + 1

        save.Quickload ( kOriginalStructureOrNewBest )

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

        GetSequenceScores ( backbone_scores , 3 , false )
 
        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 + 2

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

        score = GetScore ()
        
        gain = record_improvement ()

        p2gains [ ns ] = gain

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

function WiggleQuadruplets ()

    local  p2gains = {}

    prev_idx = -1
    prev_prev_idx = -1
    ns = 0

    repeat

        ns = ns + 1

        save.Quickload ( kOriginalStructureOrNewBest )

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

        GetSequenceScores ( residue_scores , 4 , true )
 
        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 + 3

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

        score = GetScore ()
        
        gain = record_improvement ()

        p2gains [ ns ] = gain

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

function RebuildQuadruplets ( threshold )

    local  p2gains = {}

    prev_idx = -1
    prev_prev_idx = -1
    ns = 0

    repeat

        ns = ns + 1

        save.Quickload ( kOriginalStructureOrNewBest )

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

        GetSequenceScores ( residue_scores , 4 , true )
 
        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 + 3

        prev_prev_idx = prev_idx
        prev_idx = start_idx
  
        print ( "rb " .. start_idx .. "-" .. end_idx )

        selection.DeselectAll ( )
        selection.SelectRange ( start_idx , end_idx )
          
        rebuild ( start_idx , end_idx )

        score = GetScore ()
        
        gain = record_improvement ()

        p2gains [ ns ] = gain

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


function get_parameters ()

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

    dlog.leave_sheets_alone = dialog.AddCheckbox ( "Leave sheet backbones alone" , leave_sheet_backbones_alone )  
    dlog.leave_helices_alone = dialog.AddCheckbox ( "Leave helix backbones alone" , leave_helix_backbones_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
    
        leave_sheet_backbones_alone = dlog.leave_sheets_alone.value
        leave_helix_backbones_alone = dlog.leave_helices_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 3.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 ( "id = Idealize" )
    print ( "rb = Rebuild" )
    print ( "ws = Wiggle Sidechains" )
    if ( mutate_not_shake == true ) then
       print ( "mu = Mutate" )
      else
        print ( "sh = Shake" )
    end

    print ( "" )

    WiggleBackbonePairs ()
    WiggleSidechains ()
    ShakeSidechains ( 10 )
    WiggleBackboneTriplets ()
    IdealizeBackboneTriplets ( kIdealityTripletsThreshold_a )
    WiggleSidechains ()
    ShakeSidechains ( 20 )
    WiggleQuadruplets ()  
    ShakeSidechains ( n_residues )
    RebuildQuadruplets ( kRebuildQuadrupletsThreshold )  
    IdealizeBackboneTriplets ( kIdealityTripletsThreshold_b )
    RebuildQuadruplets ( -9999999 )  

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 shape. 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.

Changes since 2.0

Incorporates elements from MicroIdealize
More efficient and frequent use of Shake and Mutate.
Allows protection of sheets and helices from rebuild.

Generally these changes enable faster recovery.

[WBarme1234]

Hi,
rebuild seems to hang and repeating the same three rebuild areas again and again

example
rb 47-50
rb 23-26
rb 22-25

rb 47-50
rb 23-26
rb 22-25

rb 47-50
rb 23-26
rb 22-25

thanks for soon correction

regards

[Bruno Kestemont]

it's really a great recipe.

tvdl handled this kind of problem (always repeating the same zones) by allowing user to skip a number of worst zones.

I suggest something like, for several zones, a kind of list of "alreadydone". This "alreadydone" would be true only during 1 or several loops, so that the recipe still can come back to these zones after other zones were managed.
(the recipe is quite gaining when returning to worst segments several times - they should not be skipped forever)

[GaryForbis]

could you add a deselect all in the cleanup routine?