Code
--[[
* Sheet Atom Bander
* Original Author: Brow42
* Version 1.0 May 1 2014
* Very experimental banding script for sheets
*
* Pick 2 sheets (by selecting or freezing them).
* Script will try to guess if they are anti-parallel or
* parallel, and band the actual acceptors and hydrogens.
* (not just the center carbons)
*
* The CLOSEST pair of AAs will determine which residues
* one one strand band to which residues on the other.
* (the script will assume the closest pair should bond)
*
* It may try to bond both sides of the stand, but then
* delete bands over 7 A long.
*
* If the bands look wrong, you can examine them in atom
* view (CTRL-SHIFT-V). Return to cartoon with CTRL-SHIFT-M.
*
* Glycines and prolines are not banded. Beta bulges and cis
* peptide bonds disrupt the sheet bonding, and the script
* doesn't understand this. If you have a beta bulge or
* cis bond, make that segment not a sheet and bonde the
* sheets on each side separately.
*
* All bands are normal strength.
--]]
-- options
BandLengthCutoff = 7
BondBandLength = 1.8
CenterBandLength = 3.0
-- globals
title = "Sheet Atom Bander v 1.0"
nSeg = structure.GetCount()
-- ================================== Begin New Dialog Library Functions
-- Add a wall of text from a table
function dialog.AddLabels(d,msg,nlabels) -- pass in # of existing autolabels
local nlabels = nlabels or #(d._Order or {}) -- default, valid if never delete dialog elements
if type(msg) == 'string' then
msg = { msg }
end
for i = 1,#msg do
d['autolabel'..tostring(i+nlabels)] = dialog.AddLabel(msg[i])
end
end
-- Create but don't display a wall of text and 1 or 2 buttons
function dialog.CreateMessageBox(msg,title,buttontext1,buttontext0)
title = title or ''
local d = dialog.CreateDialog(title)
dialog.AddLabels(d,msg)
buttontext1 = buttontext1 or 'Ok'
d.button = dialog.AddButton(buttontext1,1)
if buttontext0 ~= nil then d.button0 = dialog.AddButton(buttontext0,0) end
return d
end
-- Display a dialog box
function dialog.ShowMessageBox(msg,title,buttontext1,buttontext0)
return dialog.Show(dialog.CreateMessageBox(msg,title,buttontext1,buttontext0))
end
-- Display a box AND print to the output
function ErrorMessage(msg,buttontext1,buttontext0)
if type(msg) == 'string' then
msg = { msg }
end
for i = 1,#msg do
print(msg[i])
end
return dialog.ShowMessageBox(msg,title,buttontext1,buttontext0)
end
-- ======================================= End Dialog Functions
unpack = unpack or table.unpack -- for 5.2
-- ========================== Begin AA Atom Database ================
fsl = fsl or {} -- make sure fsl namespace exists
fsl.atom = {}
-- Take just these 4 things if all you need is to know if the sidechain starts at 5 or 6 (tail terminal = true)
fsl.atom.atomcount = { -- shorter table just for identifying terminals
a=10, c=11, d=12, e=15, f=20, g=7, h=17, i=19, k=22, l=19,
m=17, n=14, p=15, q=17, r=24, s=11, t=14, v=16, w=24, y=21
}
-- Pass in a segment number of a cysteine
function fsl.atom._IsDisulfideBonded(iSeg)
local s = current.GetSegmentEnergySubscore(iSeg,'disulfides')
return tostring(s) ~= '-0'
end
function fsl.atom._IsTerminalTest(aa,count,disulfide)
local dbvalue = fsl.atom.atomcount[aa]
if dbvalue == nil then error('Bad argument to an fsl.atom function (not an amino acid code)') end
local diff = count - dbvalue
if disulfide then diff = diff + 1 end -- because count is one less because a H was removed
if diff == 0 then return false, false, disulfide
elseif diff == 1 then return false, true, disulfide
elseif diff == 2 then return true, false, disulfide
elseif diff == 3 then return true, true, disulfide
end
error('Strange atom count, report this!')
end
-- Determine if segment is start or end of polypeptide by looking for excess atoms
function fsl.atom.IsTerminal(iSeg)
if structure.GetSecondaryStructure(iSeg) == 'M' then return false,false,false end
local aa = structure.GetAminoAcid(iSeg)
local count = structure.GetAtomCount(iSeg)
return fsl.atom._IsTerminalTest(aa,count,fsl.atom._IsDisulfideBonded(iSeg))
end
-- ====== now the rest of Atom Tables =====
fsl._nul = {} -- don't need so many CONSTANT empty tables
-- This is a table of atom classifications for AAs in a 'standard' state
-- The query functions apply adjustments due to bonds that change atom count
fsl.atom.db = {} -- to save typing, these are stored as an array { polar-H, donor, acceptor, total count }
fsl.atom.db['a'] = { {6}, fsl._nul, fsl._nul, 10 } -- alanine
fsl.atom.db['c'] = { {7}, fsl._nul, fsl._nul, 11 } -- cysteine -- NB 10 if S-S bonded!
fsl.atom.db['d'] = { {9}, fsl._nul, {7,8}, 12} -- aspartate
fsl.atom.db['e'] = { {10}, fsl._nul, {8,9}, 15 } -- glutamate
fsl.atom.db['f'] = { {12}, fsl._nul, fsl._nul, 20 } -- phenylalanine
fsl.atom.db['g'] = { {5}, fsl._nul, fsl._nul, 7 } -- glycine
fsl.atom.db['h'] = { {11,17}, {10}, {7}, 17 } -- histidine
fsl.atom.db['i'] = { {9}, fsl._nul, fsl._nul, 19 } -- isoleucine
fsl.atom.db['k'] = { {10,20,21,22}, {9}, fsl._nul, 22 } -- lysine
fsl.atom.db['l'] = { {9}, fsl._nul, fsl._nul, 19 } -- leucine
fsl.atom.db['m'] = { {9}, fsl._nul, fsl._nul, 17 } -- methionine
fsl.atom.db['n'] = { {9,13,14}, {8}, {7}, 14 } -- asparagine
fsl.atom.db['p'] = { fsl._nul, fsl._nul, fsl._nul, 15 } -- proline
fsl.atom.db['q'] = { {10,16,17}, {9}, {8}, 17 } -- glutamine
fsl.atom.db['r'] = { {12,20,21,22,23,24}, {8, 10, 11}, fsl._nul, 24 } -- arginine
fsl.atom.db['s'] = { {7,11}, {6}, {6}, 11 } -- serine
fsl.atom.db['t'] = { {8,11}, {6}, {6}, 14 } -- theronine
fsl.atom.db['v'] = { {8}, fsl._nul, fsl._nul, 16 } -- valine
fsl.atom.db['w'] = { {15,20}, {9}, fsl._nul, 24 } -- tryptophan
fsl.atom.db['y'] = { {13,21}, {12}, {12}, 21 } -- tyrosine
-- create symbolic table entries
for i,v in pairs(fsl.atom.db) do
v.polarh, v.donor, v.acceptor, v.count = unpack(v)
end
-- returns nil if it's a ligand, else returns the DB entry and aa code
-- This function lets us pass either the number of an actual segment, or
-- just look up the properties by AA code. It also handles ligand segments and
-- upper/lower case input
function fsl.atom._NumberOrCode(iSeg)
local aa
if type(iSeg) == 'number' then
local type = structure.GetSecondaryStructure(iSeg)
if type == 'M' then return end
aa = structure.GetAminoAcid(iSeg)
else aa = string.lower(iSeg)
end
local c = fsl.atom.db[aa]
if c == nil then error('Bad argument to an fsl.atom function (not an amino acid code)') end
return c, aa
end
-- ==================== End Atom Tables (subset)
-- Begin Main
print(title)
ERROR_BOTH = -1
function FindSelected()
local appending = false
local f,s,e
local range
local list = {}
for i = 1, nSeg do
f,s = freeze.IsFrozen(i), selection.IsSelected(i)
e = structure.GetSecondaryStructure(i) == 'E'
if s and f then return ERROR_BOTH end
if (f or s) and e then
if appending then range[2] = i
else
range = {i,i}
list[#list+1] = range
appending = true
end
else
if appending then appending = false end
end
end
return list
end
list = FindSelected()
if list == ERROR_BOTH then
rc = ErrorMessage("Used both selection and freezing. Clear both?","Yes","No")
if rc == 1 then
selection.DeselectAll()
freeze.UnfreezeAll()
end
return
end
if #list == 0 then
ErrorMessage({"Nothing selected. Choose 2 sheets","by freezing or selecting."})
return
end
if #list ~= 2 then
ErrorMessage("Need exactly 2 sheets selected.")
return
end
for i = 1,#list do print(i,list[i][1],list[i][2]) end
d = math.huge
for i = list[1][1],list[1][2] do
for j = list[2][1],list[2][2] do
d1 = structure.GetDistance(i,j)
if d1 < d then s1,s2,d = i,j,d1 end
end
end
print("Closest pair:",s1,s2,d)
function GetSafeDistance(s1,s2)
if s1 < list[1][1] or s1 > list[1][2] then return math.huge end
if s2 < list[2][1] or s2 > list[2][2] then return math.huge end
return structure.GetDistance(s1,s2)
end
PARALLEL,ANTIPARALLEL = 1,2
d = math.huge
d1 = GetSafeDistance(s1+1,s2+1)
if d1 < d then dir,d=PARALLEL,d1 end
d1 = GetSafeDistance(s1-1,s2-1)
if d1 < d then dir,d=PARALLEL,d1 end
d1 = GetSafeDistance(s1+1,s2-1)
if d1 < d then dir,d=ANTIPARALLEL,d1 end
d1 = GetSafeDistance(s1-1,s2+1)
if d1 < d then dir,d=ANTIPARALLEL,d1 end
print("Direction =", (dir==1 and "Parallel" or "Anti-Parallel"))
list = {{list[1][1]-s1,list[1][2]-s1},{list[2][1]-s2,list[2][2]-s2}}
if dir == PARALLEL then
d1,d2 = 1,1
start = math.max(list[1][1],list[2][1])
finish = math.min(list[1][2],list[2][2])
else
d1,d2 = 1,-1
start = math.max(list[1][1],-list[2][2])
finish = math.min(list[1][2],-list[2][1])
if start % 2 == 1 then start=start+1 end
end
print("Strand 1 banded from",s1+start*d1,"to",s1+finish*d1)
print("Strand 2 banded from",s2+start*d2,"to",s2+finish*d2)
function FirstH(iSeg)
local db, aa = fsl.atom._NumberOrCode(iSeg)
local _,last = fsl.atom.IsTerminal(iSeg)
local shift = 0
if last then shift = 1 end
if aa == 'g' then return 5+shift end
if aa == 'p' then return 9 + shift end -- but one of the backbone H is missing
return db.polarh[1] + shift
end
function BandableResidue(iSeg)
aa = structure.GetAminoAcid(iSeg)
if aa == 'g' or aa =='p' then return false end
return true
end
-- band the i_th pair
function BandAntiParallel(i)
local iBandr
seg1,seg2 = s1+d1*i,s2+d2*i
if BandableResidue(seg1) and BandableResidue(seg2) then
--print(seg1,fsl.atom.GetPolarHydrogens(seg1)[1],seg2,fsl.atom.GetPolarHydrogens(seg2)[1])
iBand = band.AddBetweenSegments(seg1,seg2,4,FirstH(seg2))
band.SetGoalLength(iBand,BondBandLength)
if band.GetLength(iBand) > BandLengthCutoff then band.Delete(iBand) end
iBand = band.AddBetweenSegments(seg1,seg2,FirstH(seg1),4)
band.SetGoalLength(iBand,BondBandLength)
if band.GetLength(iBand) > BandLengthCutoff then band.Delete(iBand) end
end
print("Pair (acceptor, donor):",seg1,seg2)
print()
band.AddBetweenSegments(seg1,seg2)
band.SetGoalLength(band.GetCount(),CenterBandLength)
end
-- band the i_th pair
function BandParallel(i)
local iBand
local seg1,seg2 = s1+d1*i,s2+d2*i
local seg3 = seg2+1
local seg4 = seg2-1
--print(seg1,fsl.atom.GetPolarHydrogens(seg1)[1],seg2,fsl.atom.GetPolarHydrogens(seg2)[1])
if BandableResidue(seg1) then
if seg3 <= nSeg and BandableResidue(seg3) then
iBand = band.AddBetweenSegments(seg1,seg3,4,FirstH(seg3))
band.SetGoalLength(iBand,BondBandLength)
if band.GetLength(iBand) > BandLengthCutoff then band.Delete(iBand) end
end
if seg4 >= 1 and BandableResidue(seg4) then
iBand = band.AddBetweenSegments(seg1,seg4,FirstH(seg1),4)
band.SetGoalLength(iBand,BondBandLength)
if band.GetLength(iBand) > BandLengthCutoff then band.Delete(iBand) end
end
end
print("Pairs (acceptor,donor):",seg1,seg3," ",seg4,seg1)
print()
band.AddBetweenSegments(seg1,seg2)
band.SetGoalLength(band.GetCount(),CenterBandLength)
end
if dir == PARALLEL then
Band = BandParallel
else
Band = BandAntiParallel
end
for i = start,finish,1 do
Band(i)
end
