Code
--[[
===================================
* Bridge Wiggle
* Original Author: Brow42
* Version 1.0 Jan. 13 2012
* This uses bands-to-space to randomly move the ends
* of your cysteines around, keeping any improvements
* in the total disulfide score.
* Can also perform an initial banding of cysteines to
* create the bonds.
* Version 1.0.1 Mar. 5 2012 Brow42
* Now autodetects cysteine pairings
* Version 1.1.0 June 16 2012 Brow42
* Now handles first/last cysteine in protein properly
* Removes wiggle bands on cancel
* New option boxes: band S-S together, wiggle backbone, use sidechain score
* Fixed incorrect random band (only did hemispheres)
* Auto-selection of possible bridges expanded to 4.0 A
* Version 1.2 May 26 2013 Brow42
* Adjacent cysteines are not considered for pairing
* Initial bands are preserved
* Version 1.2.1 May 26 2013
* Fixed GetSulfurAtom for unbonded cys
* Version 1.2.2 Oct 29, 2013
* Fixed crash if cysteine is last segment
--]]
-- =================================== Begin Defaults
options={}
options.max_str = 10.0
options.max_len = 0.1
options.cysteines = '' -- '2,20 14,26, 24,32' -- in pairs, defaults for CASP ROLL 3
options.wiggles = 5
options.max_max_len = 0.5 -- for the slider
options.ssband = false
options.bbwiggle = false
options.scscore = false
-- some globals
version = '1.2.2'
title = 'Bridge Wiggle v'..version
nseg = structure.GetCount()
qsave = 3 -- quicksave to use
ndigits = 5 -- digits to print
energy = current.GetEnergyScore() -- the starting score
nSeg = structure.GetCount()
-- ================================== End Defaults
-- ================================== Begin New Dialog Library Functions
--[[
Throws up a dialog containing just text, provided as a string or table
of strings in the first argument. Title and buttons are optional. Return
value is 1 if the first button is clicked and 0 if the second button is clicked
or the window is closed.
--]]
function dialog.MessageBox(msg,title,buttontext1, buttontext0, buttontext2)
title = title or ''
local d = dialog.CreateDialog(title)
if type(msg) == 'string' then d['1'] = dialog.AddLabel(msg)
else for i = 1,#msg do
d[tostring(i)] = dialog.AddLabel(msg[i])
end
end
buttontext1 = buttontext1 or 'Ok'
d.button = dialog.AddButton(buttontext1,1)
if buttontext2 then d.button2 = dialog.AddButton(buttontext2,2) end
if buttontext0 then d.button0 = dialog.AddButton(buttontext0,0) end
return dialog.Show(d)
end
--[[
Add items to the dialog from a table in the order given. The each table entry
is the dialog control function, the control label, the data key name, and the control options.
If this function fails, it gives an error saying which table item it failed on.
d = a created dialog, fields = table of dialog controls as described, object = data table
--]]
function dialog.AddFromTable(d,fields,object)
local func, label, id, value -- renamed field parameters for clarity
-- function wrappers to catch errors, explicitly named for clarity
function _AddLabel(args) d[id] = dialog.AddLabel(label) end
function _AddCheckbox(args) d[id] = dialog.AddCheckbox(label,value) end
function _AddSlider(args) d[id] = dialog.AddSlider(label,value,args[4],args[5],args[6]) end
function _AddTextbox(args) d[id] = dialog.AddTextbox(label,value) end
for i = 1, #fields do
local rc = true
local err
rc,err = pcall( function() func,label,id = fields[i][1],fields[i][2],fields[i][3] end)
value = object[id]
if func == dialog.AddLabel then id='_AL'..tostring(i) rc,err = pcall(_AddLabel,fields[i])
else
-- this is a crasher for AddTextbox but not an error for AddLabel!
if value == nil then error('Missing data field for field #'
..tostring(i)..'\n(Did you forget to pass a string?)') end
if func == dialog.AddCheckbox then rc,err = pcall(_AddCheckbox,fields[i])
elseif func == dialog.AddSlider then rc,err = pcall(_AddSlider,fields[i])
elseif func == dialog.AddTextbox then rc,err = pcall(_AddTextbox,fields[i])
else error('Unknown control in dialog field #'..tostring(i)) end
end -- if AddLabel
if rc == false then
error('Error setting dialog field #'..tostring(i)..(err and ('\n'..err) or ''))
end
end
end
--[[
This just copies the data from the dialog to the data table.
--]]
function dialog.ReadFields(d,fields,object)
for i = 1,#fields do
if fields[i][1]~= dialog.AddLabel then
object[fields[i][3]] = d[fields[i][3]].value
end
end
end
-- ================================== End New Dialog Library Functions
-- ================================== Begin Bridge Wiggle Dialogs
function DoInfoDialog()
local msg = {
' This uses bands-to-space to randomly move the ends',
'of your cysteines around, keeping any improvements',
'in the total disulfide score.',
' Bands will have random strengths and lengths up to',
'the values you specify. You can also include an',
'unbanded wiggle periodically. Bridge wiggle runs until',
'canceled. You can get the best disulfide scoring',
'configuration from quicksave # '..tostring(qsave)..'.',
' You can also perform an initial banding (at specified',
'maximum strength) of cysteines to create the bridges.',
' This recipe will remove bands.'
}
dialog.MessageBox(msg,'About '..title,'Back')
end
function DoOptionsDialog(options,defaults)
-- Short names for the dialog controls
-- this is not an invitation to override, they are not used as functions
local label = dialog.AddLabel -- no data member for this
local box = dialog.AddCheckbox
local slider = dialog.AddSlider -- 4 control options for this
local text = dialog.AddTextbox
-- Lay out the dialog here
-- { { dialog control, textlabel, data member name, control options }, ... }
-- all data member names must be unique
local fields ={
{label,'Enter pairs of cysteines here:'},
{text, 'Cysteines:','cysteines'},
{label,''},
{slider,'Max Str:','max_str',0.1,10.0,1},
{slider,'Max Len:','max_len',0.0,options.max_max_len,2},
{label,'Unbanded'},
{slider,'Wiggle Every:','wiggles',0,20,0},
{box,'Add band between sulfurs','ssband'},
{box,'Wiggle backbone too','bbwiggle'},
{box,'Include sidechain scoring','scscore'}
}
local d = dialog.CreateDialog(title)
dialog.AddFromTable(d,fields,options) -- last arg = table to read from
-- You can add more items from a different table here as long
-- as there is no name conflict in the two tables
-- Add the button now
d.ok = dialog.AddButton('Start!',1) -- 1 = okay
d.cancel = dialog.AddButton('Cancel',0) -- 0 = quit
d.help = dialog.AddButton('About',-1) -- < 0 = show options again, optional button
d.more = dialog.AddButton('More',-3)
if defaults then d.reset = dialog.AddButton('Reset',-2) end -- if defaults provided, allow reset
local button = dialog.Show(d)
if button == 0 then return 0 end -- 0 = cancel
dialog.ReadFields(d,fields,options) -- last arg = table to save to
return button
end
-- This is for the intial banding
function DoInitDialog()
local msg = {
'You can perform an initial banding here. Click "Band"',
'to do this. This will create a squashed pyramid of',
'bands of fixed length to try to force the cysteine pairs',
'you gave into bridges with only right angles. Click',
'"B & W" to also wiggle after banding. Existing bands',
'will not be removed for this step.',
'',
'Check below to add the fourth band that will force',
'the two cysteines to be at right angles to each other',
'or clear it to allow them to rotate freely around the',
'disulfide bridge.'
}
local d = dialog.CreateDialog('Initial Banding')
for i = 1,#msg do
d[tostring(i)] = dialog.AddLabel(msg[i])
end
d.okay = dialog.AddButton('Band',1)
d.b2 = dialog.AddButton('B & W',2)
d.b0 = dialog.AddButton('Cancel',0)
d.force = dialog.AddCheckbox('Force all right-angle cysteines?',false)
local rc = dialog.Show(d)
return rc,d.force.value
end
function DupeCheck(tab) -- return 1 if dupes
local tmp = {}
for i,v in pairs(tab) do
if tmp[v] then return 1 end
tmp[v] = true
end
return 0
end
function ErrorMessage(msg)
print(msg)
dialog.MessageBox(msg,'','Oops')
end
-- This check for input errors in the text box
-- ensure pairs of cysteines now so we don't have to check later
function CheckInput(txt)
local rc = 1
cys = {}
for v in string.gfind(txt,'-?%d+') do -- find numbers
v = tonumber(v)
if (v < 1 or v > nseg) then
ErrorMessage('Segment out of range 1-'..tostring(nseg))
return 0,{}
end
if structure.GetAminoAcid(v) ~= 'c' then
ErrorMessage('Segment '..tostring(v)..' is not a Cys')
return 0,{}
end
table.insert(cys,tonumber(v))
end
if rc == 1 and (#cys == 0 or #cys % 2 == 1) then
ErrorMessage('Not an even number of Cys entered.')
rc = 0
end
if rc == 1 and DupeCheck(cys) == 1 then
ErrorMessage('Duplicate segment entered')
rc = 0
end
return rc,cys
end
-- ================================== End user input
-- =============================== Begin Banding Routines
-- Return the atom number of the sulfur in a cysteine, or nil if not cys
function GetSulfurAtom(iSeg)
if structure.GetSecondaryStructure(iSeg) == 'M' then return nil end
if structure.GetAminoAcid(iSeg) ~= 'c' then return nil end
local diff = structure.GetAtomCount(iSeg) - 11
if current.GetSegmentEnergySubscore(iSeg,'disulfides') ~= tonumber('-0') then diff = diff + 1 end
if diff == 1 or diff == 3 then return 7 end
return 6
end
-- -500 if there's no bond
function GetDisulfideScore(iSeg)
local x = current.GetSegmentEnergySubscore(iSeg,'disulfides')
local y = 0
if options.scscore then y = current.GetSegmentEnergySubscore(iSeg,'sidechain') end
return (tostring(x) == '-0') and -500 or x+y
end
-- score each pair and also the sum
function GetScores(tab)
sc = {}
s = 0
for i = 1, #tab-1,2 do
table.insert(sc, GetDisulfideScore(tab[i]) + GetDisulfideScore(tab[i+1]))
s = s + sc[#sc]
end
return s,sc
end
-- Combining 2 functions into one
function band.SetParams(iBand, len, str)
band.SetGoalLength(iBand,len)
band.SetStrength(iBand,str)
end
-- Make a random band to space
function band.RandomBand(iSeg,rho,atom)
local theta = math.acos(2 * math.random() - 1)
local phi = 2 * math.pi * math.random()
local s2,s3
if iSeg == 1 then s2,s3 = 2,3
elseif iSeg == nSeg then s2,s3 = nSeg-1,nSeg-2
else s2,s3 = iSeg-1,iSeg+1 end
return band.Add(iSeg,s2,s3,rho,theta,phi,atom)
end
-- Make a random band and randomize the goal and weight
function RandomParamBand(iSeg,atom)
rho = options.max_len * math.random()^(1/3) +0.01
w = math.random() * options.max_str
if w < 0.1 then w = 0.1 end
local iBand = band.RandomBand(iSeg,rho,atom)
band.SetParams(iBand,0.0,w)
end
-- Make a random short band at last C and S in each cysteine
function BandCys(seg1)
local s = GetSulfurAtom(seg1)
RandomParamBand(seg1,s-1)
RandomParamBand(seg1,s)
if options.bbwiggle then RandomParamBand(seg1,0) end
end
function BandSS(seg1,seg2)
local s1,s2 = GetSulfurAtom(seg1), GetSulfurAtom(seg2)
iBand = band.AddBetweenSegments(seg1,seg2,s1,s2)
if (iBand) then
band.SetStrength(iBand,options.max_str)
band.SetGoalLength(iBand,2.05)
end
end
function structure.GetCysteines()
local n = structure.GetCount()
local cys = {}
for i = 1,n do if structure.GetAminoAcid(i) == 'c' then cys[#cys+1] = i end end
return cys
end
-- Get all pairwise distances between sulfur atoms
-- returns a function
function GetPairwiseDistances()
local cys = structure.GetCysteines()
local d = {}
for i = 1,#cys-1 do
d[cys[i]] = {}
local s1 = GetSulfurAtom(cys[i])
for j = i+1,#cys do
local s2 = GetSulfurAtom(cys[j])
iBand = band.AddBetweenSegments(cys[i],cys[j],s1,s2)
if iBand > 0 then d[cys[i]][cys[j]] = band.GetLength(iBand) band.Delete(iBand) end
end
end
return
function(seg1,seg2)
if seg1 > seg2 then seg1,seg2 = seg2,seg1 end
return d[seg1] and d[seg1][seg2] or math.huge
end
end
-- Band two cysteines in a way that forces right angles everywhere
function DoInitialBanding(tab)
local rc, force
rc,force = DoInitDialog()
if rc == 0 then return end
print('Creating initial bands of strength',options.max_str)
local SS = 2.05
local SC = 1.8
local d2 = math.sqrt(SS*SS+SC*SC) -- SC x dir + SS y dir
local d3 = math.sqrt(d2*d2+SC*SC) -- SC x dir + SS y dir + SC z dir
for i=1,#tab-1,2 do
local s1,s2, iBand
s1,s2 = tab[i],tab[i+1]
local sulfur1,sulfur2 = GetSulfurAtom(s1), GetSulfurAtom(s2)
iBand = band.AddBetweenSegments(s1,s2,sulfur1,sulfur2)
band.SetParams(iBand,SS,options.max_str)
iBand = band.AddBetweenSegments(s1,s2,sulfur1,sulfur2-1)
band.SetParams(iBand,d2,options.max_str)
iBand = band.AddBetweenSegments(s1,s2,sulfur1-1,sulfur2)
band.SetParams(iBand,d2,options.max_str)
if force then
iBand = band.AddBetweenSegments(s1,s2,sulfur1-1,sulfur2-1)
band.SetParams(iBand,d3,options.max_str)
end
end
if rc == 1 then return end
local s, score
s,score = GetScores(tab)
print('Selected pairs:')
PrintCys(tab,score,s)
local nwiggle = 10
print('Wiggling for',nwiggle,'iterations...')
structure.WiggleAll(nwiggle)
s,scor2 = GetScores(tab)
PrintCys(tab,score,s)
energy = current.GetEnergyScore()
end
-- =============================== End Banding Routines
-- ================================ Begin Utility Routines
-- Recursive table printer
function PrintTable(tab,indent)
indent = indent or ''
for i,v in pairs(tab) do
if type(v) == 'table' then do print(indent,i) PrintTable(v,indent..' ') end
else print(indent,i,v) end
end
end
function round(x,n)
n = n or ndigits
return math.floor(0.5 + x * 10^n)/10^n
end
function PrintCys(tab,sc,tot)
local fmtstr = '%3d %3d %8.5f'
for i = 1, #sc do
print(fmtstr:format(tab[2*i-1],tab[2*i],round(sc[i])))
end
print('Total: ',tot)
print('Overall Score:',round(current.GetEnergyScore()))
print('Overall Change:',round(current.GetEnergyScore()-energy))
print('------------------------')
end
-- brute force optimization of pairs
-- Simple greed algorithm, not optimal
function FindPairs()
local pairs = {}
local cys = structure.GetCysteines()
local distance = GetPairwiseDistances()
local n = math.floor(#cys/2) -- # pairs
local d,p,min
for k = 1,n do
min = math.huge
p = nil
for i = 1,#cys-1 do
for j = i+1,#cys do
d = distance(cys[i],cys[j])
if d < min and math.abs(cys[i] - cys[j]) > 1 then min,p = d,{cys[i],cys[j],i,j} end
end
end
if p == nil then return pairs,cys end
pairs[#pairs+1] = {p[1], p[2], distance(p[1],p[2])}
table.remove(cys,p[4]) -- 4 is after 3 so remove first
table.remove(cys,p[3])
end
return pairs,cys
end
_NBands = band.GetCount()
function DeleteBands()
while band.GetCount() > _NBands do
band.Delete(band.GetCount())
end
end
-- ================================= End Utility Routines
-- ================================= Begin Main
--defaults={}
--for i,v in pairs(options) do defaults[i]=v end -- save the initial options as defaults, shallow copy
print('Finding cysteines...(a greedy algorithm, shortest bonds, non-optimal)...')
cpairs, unpaired = FindPairs()
print('Selected pairing:')
print('Seg1 Seg2 S-S Distance')
options.cysteines = ''
cut = 4.0
for i = 1,#cpairs do
print(cpairs[i][1],cpairs[i][2],cpairs[i][3])
if cpairs[i][3] < cut then
options.cysteines = string.format('%s %d,%d',options.cysteines,cpairs[i][1],cpairs[i][2])
end
end
print('Unpaired cys:',table.concat(unpaired,', '))
print('You should cut out the longest distances. (Already cut >',cut,')')
cys = {}
repeat -- Loop if something other than okay or cancel
local rc = DoOptionsDialog(options,defaults)
if rc == 0 then print('Canceled') return end -- cancel
local rc2
rc2,cys = CheckInput(options.cysteines)
if rc == -1 then DoInfoDialog() end
-- if rc == -2 then for i,v in pairs(defaults) do options[i] = v end end
if rc2 == 1 and rc == -3 then DoInitialBanding(cys) end
until rc == 1 and rc2 == 1-- okay
--print('Printing Option Table')
--PrintTable(options)
print(puzzle.GetName(),os.date(),current.GetScore())
if options.scscore then print('Scoring on disulfide + sidechain socres.')
else print('Scoring on disulfide score part.') end
s1,score1 = GetScores(cys)
print('Selected pairs:')
PrintCys(cys,score1,s1)
print('----------------------')
seed = (os.time() * 5779) % 100000 -- that's a prime
-- seed = 81944 -- for debugging, set to value where the error was
print('Random seed = ',seed)
math.randomseed(seed)
print('Saves are in Quicksave '..tostring(qsave))
print('Resetting recent best...')
recentbest.Save()
-- Repeat forever
function WiggleForever()
iLoop = 1
repeat
DeleteBands()
s1,score1 = GetScores(cys)
save.Quicksave(qsave)
for i = 1,#cys do
BandCys(cys[i])
end
if options.ssband then
for i = 1,#cys,2 do BandSS(cys[i],cys[i+1]) end
end
structure.WiggleAll(1)
s2,score2 = GetScores(cys)
if s2 < s1 then save.Quickload(qsave)
else print('Loop',iLoop) PrintCys(cys,score2,s2) end
iLoop = iLoop + 1
if options.wiggles > 0 and (iLoop % options.wiggles +1 == options.wiggles) then
-- print('Doing an unbanded wiggle')
DeleteBands()
structure.WiggleAll(1)
end
until false
end
-- Routine called on cancel or error
function OnError(err)
if err:find('Cancel') then print('Cancelled') DeleteBands()
else print('Error: ',err)
end
return err
end
rc,err = xpcall(WiggleForever,OnError)
