Icon representing a recipe

Recipe: Loop rebuild 3

created by spvincent

Profile


Name
Loop rebuild 3
ID
41540
Shared with
Public
Parent
None
Children
Created on
May 11, 2012 at 19:09 PM UTC
Updated on
May 11, 2012 at 19:09 PM UTC
Description

See comment

Best for


Code


rebuild_length = 6 -- Length of sections to rebuild. number_of_rebuilds = 8 -- It'll look for the best score from this number of rebuilds. 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 loops_only = true -- If true, will still allow a single non-loop residue in the rebuild segment -- If false, non-loop residues will be temporarily converted to loops min_residue = 1 -- Change these if you want 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 9999 (or any large number) respectively max_residue = 999 kWorstAllowableAngle = 75 -- ignore rebuilds with backbone angles less than this (degrees) kMaxSphereDistance = 14 -- Residues further than this away won't be shaken -- Use of quicksave slots kOriginalStructureOrNewBest = 1 kBestPostRebuild = 2 kTemp = 3 -- Globals original_secondary_structure = {} fixed_residues = {} --wiggle_accuracy = 0 n_residues = 0 best_score = 0 cos_worst_allowable_angle = 0 -- NudgeWiggle : adaptive thresholds. Change the thresholds based on current protein. kInitNudgeWiggleThreshold = 25 nudge_wiggle_gains = {} kLowestNudgeWiggleThreshold = 8 -- Performance monitoring. Result intended for import into Excel as comma-delimited text monit = false monit_init_worst_angle = 0 monit_worst_angle = 0 monit_score_after_rebuild = 0 monit_score_after_post_rebuild_shake = 0 monit_gain_from_possible_initial_nudge_wiggle = 0 monit_gain_from_wiggle = 0 monit_gain_from_shake = 0 monit_gain_from_nudge_wiggle = 0 monit_total_gain = 0 -- Begin function r3 ( x ) -- Round to 3 decimal places t = 10 ^ 3 return math.floor ( x*t + 0.5 ) / t end function GetScore () score = current.GetEnergyScore () 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 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 GetListOfFixedResidues () for i = 1 , n_residues do fixed_residues [ i ] = false end for i = 1 , n_residues do seg_backbone_frozen, seg_sidechain_frozen = freeze.IsFrozen ( i ) if ( seg_backbone_frozen == true ) then fixed_residues [ i ] = true end s = structure.IsLocked ( i ) if ( s == true ) then fixed_residues [ i ] = true end end end function DoesSectionContainFixedResidues ( start_idx , end_idx ) for i = start_idx , end_idx do if ( fixed_residues [ i ] == true ) then return true end end return false end function GetNudgeWiggleThreshold () n = #nudge_wiggle_gains k = math.max ( n - 5 , 1 ) local x = 0 for i = k , n do x = x + nudge_wiggle_gains [ i ] end av = x / ( n - k + 1 ) return math.max ( kLowestNudgeWiggleThreshold , 1.4*av ) end function MiniFuse () selection.SelectAll () score_in_kTemp = GetScore() save.Quicksave ( kTemp) selection.SelectAll () behavior.SetClashImportance ( 0.07 ) structure.ShakeSidechainsSelected ( 1 ) behavior.SetClashImportance ( 1 ) structure.WiggleSelected ( 8 ) score = GetScore() if ( score < score_in_kTemp ) then save.Quickload ( kTemp ) end end function NudgeWiggle ( score ) selection.SelectAll () score_in_kTemp = GetScore() save.Quicksave ( kTemp) -- behavior.SetWiggleAccuracy ( 2 ) behavior.SetClashImportance ( 0.4 + 0.2 * pseudorandom ( score ) ) structure.WiggleSelected ( 1 ) -- behavior.SetWiggleAccuracy ( wiggle_accuracy ) behavior.SetClashImportance ( 1 ) structure.WiggleSelected ( 6 ) score = GetScore() gain_from_nw = score - score_in_kTemp if ( monit ) then monit_gain_from_nudge_wiggle = gain_from_nw end if ( score < score_in_kTemp ) then save.Quickload ( kTemp ) end end function WiggleShakeNudgeWiggle () -- returns the score after doing a WSNW: short-circuiting if things don't look promising score_after_rebuild_shake = GetScore () structure.WiggleSelected ( 10 ) score_after_wiggle = GetScore () if ( monit ) then monit_gain_from_wiggle = score_after_wiggle - score_after_rebuild_shake end if ( score_after_wiggle - score_after_rebuild_shake < 10 ) then -- The wiggle got stuck and didn't achieve anything NudgeWiggle ( score_after_wiggle ) score_after_nudge_wiggle = GetScore () if ( monit ) then monit_gain_from_possible_initial_nudge_wiggle = score_after_nudge_wiggle - score_after_rebuild_shake end score_after_wiggle = score_after_nudge_wiggle end nwt = GetNudgeWiggleThreshold () if ( best_score - score_after_wiggle > 1.5*nwt ) then return score_after_wiggle end structure.ShakeSidechainsSelected ( 1 ) score_after_second_shake = GetScore () if ( monit ) then monit_gain_from_shake = score_after_second_shake - score_after_wiggle end if ( best_score - score_after_second_shake < nwt ) then -- nwt = GetNudgeWiggleThreshold () start_score = GetScore() NudgeWiggle ( score_after_second_shake ) MiniFuse () -- print ( "nwt " .. r3 ( nwt ) ) score = GetScore () if ( score - start_score > 0 ) then table.insert ( nudge_wiggle_gains , score - start_score ) end else score = score_after_second_shake end return score end function post_rebuild () save.Quickload ( kBestPostRebuild ) if ( monit ) then monit_score_after_post_rebuild_shake = GetScore () end curr_score = WiggleShakeNudgeWiggle () if ( monit ) then monit_total_gain = curr_score - best_score end if ( curr_score > best_score ) then best_score = curr_score save.LoadSecondaryStructure () save.Quicksave ( kOriginalStructureOrNewBest ) print ( "Improvement to ".. r3 ( best_score ) ) 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 get_worst_angle ( start_idx , end_idx ) -- For monitoring only cos_worst_angle = -1 for i = start_idx , end_idx do cos_worst_angle = math.max ( cos_worst_angle , get_cos_backbone_angle ( i ) ) end return math.deg ( math.acos ( cos_worst_angle ) ) 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 return false end end return true end function ConstructListOfResiduesWithinRange ( SphereSegList , start_idx , end_idx ) mid_idx = math.floor ( ( end_idx + start_idx ) / 2 ) -- Look at start_idx, end_idx, and mid_idx only to save time at the expense of a small amount of accuracy for i = 1 , n_residues do if ( ( i >= start_idx and i <= end_idx ) or ( structure.GetDistance ( i , mid_idx ) < kMaxSphereDistance ) or ( structure.GetDistance ( i , end_idx ) < kMaxSphereDistance ) or ( structure.GetDistance ( i , start_idx ) < kMaxSphereDistance ) ) then table.insert ( SphereSegList , i ) end end -- print ( "n side chains " , #SphereSegList ) end function SelectSphere ( SphereSegList ) selection.DeselectAll () for i = 1 , #SphereSegList do selection.Select ( SphereSegList [ i ] ) end 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 SphereSegList = {} -- Speed things up when shaking ConstructListOfResiduesWithinRange ( SphereSegList , start_idx , end_idx ) for i = 1 , number_of_rebuilds do save.Quickload ( kOriginalStructureOrNewBest ) ConvertToLoop ( start_idx , end_idx ) selection.DeselectAll () selection.SelectRange ( start_idx , end_idx ) if ( monit ) then monit_init_worst_angle = get_worst_angle ( start_idx , end_idx ) monit_worst_angle = 0 monit_score_after_rebuild = 0 monit_score_after_post_rebuild_shake = 0 end -- 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 math.abs ( score_after_rebuild - best_score ) >= 1e-5 or k > 2 if ( monit ) then monit_score_after_rebuild = score_after_rebuild monit_worst_angle = get_worst_angle ( start_idx , end_idx ) end if ( ( math.abs ( score_after_rebuild - best_score ) >= 1e-5 ) and ( is_backbone_reasonably_angled ( start_idx , end_idx ) ) ) then good_structure_found = true SelectSphere ( SphereSegList ) structure.ShakeSidechainsSelected ( 1 ) selection.SelectAll () score_after_shake = GetScore() if ( monit ) then monit_score_after_post_rebuild_shake = score_after_shake end 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 if ( monit ) then print ( r3 ( best_score ) .. " , " .. start_idx .. " , " .. end_idx .. " , " .. r3 ( monit_score_after_rebuild ) .. " , " .. r3 ( monit_score_after_post_rebuild_shake ) .. " , " .. " , " .. " , " .. " , " .. " , " .. " , " .. r3 ( monit_init_worst_angle ) .. " , " .. r3 ( monit_worst_angle ) .. " , " .. r3 ( get_worst_angle ( start_idx , end_idx ) ) ) end end -- for i if ( good_structure_found == true ) then return 1 else return 0 end end function rebuild_all_at_specified_length ( 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 end_idx = start_idx + rebuild_length - 1 if ( ( start_idx >= min_residue and end_idx <= max_residue ) and ( DoesSectionContainFixedResidues ( start_idx , end_idx ) == false ) and ( loops_only == false or ( loops_only == true and IsSegmentALoop ( start_idx , end_idx ) < 2 ) ) ) then if ( monit ) then monit_init_worst_angle = get_worst_angle ( start_idx , end_idx ) monit_worst_angle = 0 monit_score_after_rebuild = 0 monit_score_after_post_rebuild_shake = 0 monit_gain_from_possible_initial_nudge_wiggle = 0 monit_gain_from_wiggle = 0 monit_gain_from_shake = 0 monit_gain_from_nudge_wiggle = 0 monit_total_gain = 0 end rb = rb + 1 if ( monit == false ) then print ( "rb ".. rb .. " " .. start_idx .."-" .. end_idx .. " (" ..i .. "/" .. n_possible_segs .. ")" ) end if ( rebuild_seg ( start_idx , end_idx ) == 1 ) then post_rebuild () end if ( monit ) then print ( r3 ( best_score ) .. " , " .. start_idx .. " , " .. end_idx .. " , " .. " , " .. r3 ( monit_score_after_post_rebuild_shake ) .. " , " .. r3 ( monit_gain_from_wiggle ) .. " , " .. r3 ( monit_gain_from_possible_initial_nudge_wiggle ) .. " , " .. r3 ( monit_gain_from_shake ) .. " , " .. r3 ( monit_gain_from_nudge_wiggle ) .. " , " .. r3 ( monit_total_gain ) .. " , " .. r3 ( monit_init_worst_angle ) .. " , " .. r3 ( monit_worst_angle ) .. " , " .. r3 ( get_worst_angle ( start_idx , end_idx ) ) .. " , " .. r3 ( GetNudgeWiggleThreshold () ) ) end save.LoadSecondaryStructure () 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 function get_parameters () local dlog = dialog.CreateDialog ( "Loop rebuild 3" ) dlog.rb_length = dialog.AddSlider ( "Rebuild length" , rebuild_length , 3 , 10 , 0 ) dlog.n_rebuilds = dialog.AddSlider ( "Number of rebuilds" , number_of_rebuilds , 1 , 50 , 0 ) dlog.min_residue = dialog.AddSlider ( "Min residue" , 1 , 1 , n_residues , 0 ) dlog.max_residue = dialog.AddSlider ( "Max residue" , n_residues , 1 , n_residues , 0 ) dlog.loops = dialog.AddCheckbox ( "Loops only" , loops_only ) dlog.ok = dialog.AddButton ( "OK" , 1 ) dlog.cancel = dialog.AddButton ( "Cancel" , 0 ) if ( dialog.Show ( dlog ) > 0 ) then rebuild_length = dlog.rb_length.value number_of_rebuilds = dlog.n_rebuilds.value min_residue = dlog.min_residue.value max_residue = dlog.max_residue.value loops_only = dlog.loops.value return true else return false end end function main() band.DisableAll () n_residues = structure.GetCount () for i = 1, n_residues do original_secondary_structure [ i ] = structure.GetSecondaryStructure ( i ) end save.SaveSecondaryStructure () rb = 0 save.Quicksave ( kOriginalStructureOrNewBest ) behavior.SetClashImportance ( 1 ) --wiggle_accuracy = behavior.GetWiggleAccuracy () --if ( wiggle_accuracy < 3 ) then --wiggle_accuracy = 3 --behavior.SetWiggleAccuracy ( 3 ) --print ( "Wiggle accuracy set to 3" ) --end if ( get_parameters () == false ) then return -- graceful exit end best_score = GetScore () print ( "Start score : " .. r3 ( best_score ) ) if ( loops_only == true ) then print ( "Loops only" ) else print ( "Loops, sheets, and helices" ) end max_residue = math.min ( n_residues , max_residue ) min_residue = math.max ( 1 , min_residue ) if ( max_residue < min_residue + rebuild_length - 1 ) then print ( "Rebuild range error." ) print ( "Max residue (" .. max_residue .. ") must be greater than or equal to Min residue (" .. min_residue .. ")" ) print ( "plus Rebuild Length (" .. rebuild_length .. ")" ) return end print ( "Rebuild range " .. min_residue .. " to " .. max_residue ) GetListOfFixedResidues () 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 ) ) table.insert ( nudge_wiggle_gains , kInitNudgeWiggleThreshold ) if ( monit ) then print ( "Start score" .. " , " .. "Start idx" .. " , " .. "End idx" .. " , " .. "Score after rebuild" .. " , " .. "Score after shake" .. " , " .. "Gain from wiggle" .. " , " .. "Gain from possible nudge wiggle" .. " , " .. "Gain from shake" .. " , " .. "Gain from nudge wiggle" .. " , " .. "Total gain" .. " , " .. "Initial worst angle" .. " , " .. "Post rebuild worst angle" .. " , " .. "Final worst angle" .. " , " .. "Nw Threshold" ) end while true do rebuild_all_at_specified_length ( rebuild_length ) end save.Quickload ( kOriginalStructureOrNewBest ) band.EnableAll () end function cleanup () print ( "Cleaning up" ) behavior.SetClashImportance ( 1.0 ) save.Quickload ( kOriginalStructureOrNewBest ) save.LoadSecondaryStructure () band.EnableAll () end xpcall ( main , cleanup ) rebuild_length = 6 -- Length of sections to rebuild. number_of_rebuilds = 8 -- It'll look for the best score from this number of rebuilds. 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 loops_only = true -- If true, will still allow a single non-loop residue in the rebuild segment -- If false, non-loop residues will be temporarily converted to loops min_residue = 1 -- Change these if you want 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 9999 (or any large number) respectively max_residue = 999 kWorstAllowableAngle = 75 -- ignore rebuilds with backbone angles less than this (degrees) kMaxSphereDistance = 14 -- Residues further than this away won't be shaken -- Use of quicksave slots kOriginalStructureOrNewBest = 1 kBestPostRebuild = 2 kTemp = 3 -- Globals original_secondary_structure = {} fixed_residues = {} --wiggle_accuracy = 0 n_residues = 0 best_score = 0 cos_worst_allowable_angle = 0 -- NudgeWiggle : adaptive thresholds. Change the thresholds based on current protein. kInitNudgeWiggleThreshold = 25 nudge_wiggle_gains = {} kLowestNudgeWiggleThreshold = 8 -- Performance monitoring. Result intended for import into Excel as comma-delimited text monit = false monit_init_worst_angle = 0 monit_worst_angle = 0 monit_score_after_rebuild = 0 monit_score_after_post_rebuild_shake = 0 monit_gain_from_possible_initial_nudge_wiggle = 0 monit_gain_from_wiggle = 0 monit_gain_from_shake = 0 monit_gain_from_nudge_wiggle = 0 monit_total_gain = 0 -- Begin function r3 ( x ) -- Round to 3 decimal places t = 10 ^ 3 return math.floor ( x*t + 0.5 ) / t end function GetScore () score = current.GetEnergyScore () 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 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 GetListOfFixedResidues () for i = 1 , n_residues do fixed_residues [ i ] = false end for i = 1 , n_residues do seg_backbone_frozen, seg_sidechain_frozen = freeze.IsFrozen ( i ) if ( seg_backbone_frozen == true ) then fixed_residues [ i ] = true end s = structure.IsLocked ( i ) if ( s == true ) then fixed_residues [ i ] = true end end end function DoesSectionContainFixedResidues ( start_idx , end_idx ) for i = start_idx , end_idx do if ( fixed_residues [ i ] == true ) then return true end end return false end function GetNudgeWiggleThreshold () n = #nudge_wiggle_gains k = math.max ( n - 5 , 1 ) local x = 0 for i = k , n do x = x + nudge_wiggle_gains [ i ] end av = x / ( n - k + 1 ) return math.max ( kLowestNudgeWiggleThreshold , 1.4*av ) end function MiniFuse () selection.SelectAll () score_in_kTemp = GetScore() save.Quicksave ( kTemp) selection.SelectAll () behavior.SetClashImportance ( 0.07 ) structure.ShakeSidechainsSelected ( 1 ) behavior.SetClashImportance ( 1 ) structure.WiggleSelected ( 8 ) score = GetScore() if ( score < score_in_kTemp ) then save.Quickload ( kTemp ) end end function NudgeWiggle ( score ) selection.SelectAll () score_in_kTemp = GetScore() save.Quicksave ( kTemp) -- behavior.SetWiggleAccuracy ( 2 ) behavior.SetClashImportance ( 0.4 + 0.2 * pseudorandom ( score ) ) structure.WiggleSelected ( 1 ) -- behavior.SetWiggleAccuracy ( wiggle_accuracy ) behavior.SetClashImportance ( 1 ) structure.WiggleSelected ( 6 ) score = GetScore() gain_from_nw = score - score_in_kTemp if ( monit ) then monit_gain_from_nudge_wiggle = gain_from_nw end if ( score < score_in_kTemp ) then save.Quickload ( kTemp ) end end function WiggleShakeNudgeWiggle () -- returns the score after doing a WSNW: short-circuiting if things don't look promising score_after_rebuild_shake = GetScore () structure.WiggleSelected ( 10 ) score_after_wiggle = GetScore () if ( monit ) then monit_gain_from_wiggle = score_after_wiggle - score_after_rebuild_shake end if ( score_after_wiggle - score_after_rebuild_shake < 10 ) then -- The wiggle got stuck and didn't achieve anything NudgeWiggle ( score_after_wiggle ) score_after_nudge_wiggle = GetScore () if ( monit ) then monit_gain_from_possible_initial_nudge_wiggle = score_after_nudge_wiggle - score_after_rebuild_shake end score_after_wiggle = score_after_nudge_wiggle end nwt = GetNudgeWiggleThreshold () if ( best_score - score_after_wiggle > 1.5*nwt ) then return score_after_wiggle end structure.ShakeSidechainsSelected ( 1 ) score_after_second_shake = GetScore () if ( monit ) then monit_gain_from_shake = score_after_second_shake - score_after_wiggle end if ( best_score - score_after_second_shake < nwt ) then -- nwt = GetNudgeWiggleThreshold () start_score = GetScore() NudgeWiggle ( score_after_second_shake ) MiniFuse () -- print ( "nwt " .. r3 ( nwt ) ) score = GetScore () if ( score - start_score > 0 ) then table.insert ( nudge_wiggle_gains , score - start_score ) end else score = score_after_second_shake end return score end function post_rebuild () save.Quickload ( kBestPostRebuild ) if ( monit ) then monit_score_after_post_rebuild_shake = GetScore () end curr_score = WiggleShakeNudgeWiggle () if ( monit ) then monit_total_gain = curr_score - best_score end if ( curr_score > best_score ) then best_score = curr_score save.LoadSecondaryStructure () save.Quicksave ( kOriginalStructureOrNewBest ) print ( "Improvement to ".. r3 ( best_score ) ) 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 get_worst_angle ( start_idx , end_idx ) -- For monitoring only cos_worst_angle = -1 for i = start_idx , end_idx do cos_worst_angle = math.max ( cos_worst_angle , get_cos_backbone_angle ( i ) ) end return math.deg ( math.acos ( cos_worst_angle ) ) 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 return false end end return true end function ConstructListOfResiduesWithinRange ( SphereSegList , start_idx , end_idx ) mid_idx = math.floor ( ( end_idx + start_idx ) / 2 ) -- Look at start_idx, end_idx, and mid_idx only to save time at the expense of a small amount of accuracy for i = 1 , n_residues do if ( ( i >= start_idx and i <= end_idx ) or ( structure.GetDistance ( i , mid_idx ) < kMaxSphereDistance ) or ( structure.GetDistance ( i , end_idx ) < kMaxSphereDistance ) or ( structure.GetDistance ( i , start_idx ) < kMaxSphereDistance ) ) then table.insert ( SphereSegList , i ) end end -- print ( "n side chains " , #SphereSegList ) end function SelectSphere ( SphereSegList ) selection.DeselectAll () for i = 1 , #SphereSegList do selection.Select ( SphereSegList [ i ] ) end 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 SphereSegList = {} -- Speed things up when shaking ConstructListOfResiduesWithinRange ( SphereSegList , start_idx , end_idx ) for i = 1 , number_of_rebuilds do save.Quickload ( kOriginalStructureOrNewBest ) ConvertToLoop ( start_idx , end_idx ) selection.DeselectAll () selection.SelectRange ( start_idx , end_idx ) if ( monit ) then monit_init_worst_angle = get_worst_angle ( start_idx , end_idx ) monit_worst_angle = 0 monit_score_after_rebuild = 0 monit_score_after_post_rebuild_shake = 0 end -- 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 math.abs ( score_after_rebuild - best_score ) >= 1e-5 or k > 2 if ( monit ) then monit_score_after_rebuild = score_after_rebuild monit_worst_angle = get_worst_angle ( start_idx , end_idx ) end if ( ( math.abs ( score_after_rebuild - best_score ) >= 1e-5 ) and ( is_backbone_reasonably_angled ( start_idx , end_idx ) ) ) then good_structure_found = true SelectSphere ( SphereSegList ) structure.ShakeSidechainsSelected ( 1 ) selection.SelectAll () score_after_shake = GetScore() if ( monit ) then monit_score_after_post_rebuild_shake = score_after_shake end 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 if ( monit ) then print ( r3 ( best_score ) .. " , " .. start_idx .. " , " .. end_idx .. " , " .. r3 ( monit_score_after_rebuild ) .. " , " .. r3 ( monit_score_after_post_rebuild_shake ) .. " , " .. " , " .. " , " .. " , " .. " , " .. " , " .. r3 ( monit_init_worst_angle ) .. " , " .. r3 ( monit_worst_angle ) .. " , " .. r3 ( get_worst_angle ( start_idx , end_idx ) ) ) end end -- for i if ( good_structure_found == true ) then return 1 else return 0 end end function rebuild_all_at_specified_length ( 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 end_idx = start_idx + rebuild_length - 1 if ( ( start_idx >= min_residue and end_idx <= max_residue ) and ( DoesSectionContainFixedResidues ( start_idx , end_idx ) == false ) and ( loops_only == false or ( loops_only == true and IsSegmentALoop ( start_idx , end_idx ) < 2 ) ) ) then if ( monit ) then monit_init_worst_angle = get_worst_angle ( start_idx , end_idx ) monit_worst_angle = 0 monit_score_after_rebuild = 0 monit_score_after_post_rebuild_shake = 0 monit_gain_from_possible_initial_nudge_wiggle = 0 monit_gain_from_wiggle = 0 monit_gain_from_shake = 0 monit_gain_from_nudge_wiggle = 0 monit_total_gain = 0 end rb = rb + 1 if ( monit == false ) then print ( "rb ".. rb .. " " .. start_idx .."-" .. end_idx .. " (" ..i .. "/" .. n_possible_segs .. ")" ) end if ( rebuild_seg ( start_idx , end_idx ) == 1 ) then post_rebuild () end if ( monit ) then print ( r3 ( best_score ) .. " , " .. start_idx .. " , " .. end_idx .. " , " .. " , " .. r3 ( monit_score_after_post_rebuild_shake ) .. " , " .. r3 ( monit_gain_from_wiggle ) .. " , " .. r3 ( monit_gain_from_possible_initial_nudge_wiggle ) .. " , " .. r3 ( monit_gain_from_shake ) .. " , " .. r3 ( monit_gain_from_nudge_wiggle ) .. " , " .. r3 ( monit_total_gain ) .. " , " .. r3 ( monit_init_worst_angle ) .. " , " .. r3 ( monit_worst_angle ) .. " , " .. r3 ( get_worst_angle ( start_idx , end_idx ) ) .. " , " .. r3 ( GetNudgeWiggleThreshold () ) ) end save.LoadSecondaryStructure () 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 function get_parameters () local dlog = dialog.CreateDialog ( "Loop rebuild 3" ) dlog.rb_length = dialog.AddSlider ( "Rebuild length" , rebuild_length , 3 , 10 , 0 ) dlog.n_rebuilds = dialog.AddSlider ( "Number of rebuilds" , number_of_rebuilds , 1 , 50 , 0 ) dlog.min_residue = dialog.AddSlider ( "Min residue" , 1 , 1 , n_residues , 0 ) dlog.max_residue = dialog.AddSlider ( "Max residue" , n_residues , 1 , n_residues , 0 ) dlog.loops = dialog.AddCheckbox ( "Loops only" , loops_only ) dlog.ok = dialog.AddButton ( "OK" , 1 ) dlog.cancel = dialog.AddButton ( "Cancel" , 0 ) if ( dialog.Show ( dlog ) > 0 ) then rebuild_length = dlog.rb_length.value number_of_rebuilds = dlog.n_rebuilds.value min_residue = dlog.min_residue.value max_residue = dlog.max_residue.value loops_only = dlog.loops.value return true else return false end end function main() band.DisableAll () n_residues = structure.GetCount () for i = 1, n_residues do original_secondary_structure [ i ] = structure.GetSecondaryStructure ( i ) end save.SaveSecondaryStructure () rb = 0 save.Quicksave ( kOriginalStructureOrNewBest ) behavior.SetClashImportance ( 1 ) --wiggle_accuracy = behavior.GetWiggleAccuracy () --if ( wiggle_accuracy < 3 ) then --wiggle_accuracy = 3 --behavior.SetWiggleAccuracy ( 3 ) --print ( "Wiggle accuracy set to 3" ) --end if ( get_parameters () == false ) then return -- graceful exit end best_score = GetScore () print ( "Start score : " .. r3 ( best_score ) ) if ( loops_only == true ) then print ( "Loops only" ) else print ( "Loops, sheets, and helices" ) end max_residue = math.min ( n_residues , max_residue ) min_residue = math.max ( 1 , min_residue ) if ( max_residue < min_residue + rebuild_length - 1 ) then print ( "Rebuild range error." ) print ( "Max residue (" .. max_residue .. ") must be greater than or equal to Min residue (" .. min_residue .. ")" ) print ( "plus Rebuild Length (" .. rebuild_length .. ")" ) return end print ( "Rebuild range " .. min_residue .. " to " .. max_residue ) GetListOfFixedResidues () 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 ) ) table.insert ( nudge_wiggle_gains , kInitNudgeWiggleThreshold ) if ( monit ) then print ( "Start score" .. " , " .. "Start idx" .. " , " .. "End idx" .. " , " .. "Score after rebuild" .. " , " .. "Score after shake" .. " , " .. "Gain from wiggle" .. " , " .. "Gain from possible nudge wiggle" .. " , " .. "Gain from shake" .. " , " .. "Gain from nudge wiggle" .. " , " .. "Total gain" .. " , " .. "Initial worst angle" .. " , " .. "Post rebuild worst angle" .. " , " .. "Final worst angle" .. " , " .. "Nw Threshold" ) end while true do rebuild_all_at_specified_length ( rebuild_length ) end save.Quickload ( kOriginalStructureOrNewBest ) band.EnableAll () end function cleanup () print ( "Cleaning up" ) behavior.SetClashImportance ( 1.0 ) save.Quickload ( kOriginalStructureOrNewBest ) save.LoadSecondaryStructure () band.EnableAll () end xpcall ( main , cleanup )

Comments


spvincent Lv 1

Changes since 2.0

The principle change from previous versions is that the script selects the best rebuild in a different way. Rather than selecting the best scoring solution from a number of rebuilds, it takes the best score from the same number of rebuild/shakes. While this is likely to give a "better" rebuilt structure in the sense that it is more likely to result in a score improvement, there's considerable overhead in doing shakes so it will look at fewer structures than previously in the same amount of time.

It's not necessarily clear this is an improvement on earlier versions of this script. My impression, unconfirmed by comparative testing, is that it is better for rebuild lengths of 6 or more, not so good for 4 or fewer, and for rebuild lengths of 5 it's probably a wash.

Algorithm outline

1) Pick a section to rebuild
2) For how ever many rebuilds you wish to do, rebuild once and, unless the backbone looks badly angled in which case ignore it, shake.
3) Take the best result from 2) and do a wiggle/shake
4) If the resulting score is reasonably close to the score of the starting structure, do some limited fuzing and hope for an improvement

Dialog settings

Rebuild length: The size of segments to be rebuilt. 5, 6, or 7 are probably your best bets.

Number of rebuilds: How many rebuild/shake attempts to make.

Min residue/Max residue: The defaults are 1 and the number of residues in the protein respectively, but changing these lets you focus on a particular region.

Loops only: If unchecked, it'll pay no attention to the secondary structure. If checked, it'll rebuild sections that are loops only and sections containing a single non-loop residue.

tim1079 Lv 1

This program gets stuck on the detect secondary structure function call.

Thoughts?