Code
function Main()
--
-- makehelix 5/15/12 by jeff101
-- MakeHelix Update 10/09/20 by Formula350
--
tot=structure.GetCount() -- find total # of residues
if tot<3 then
print('Total='..tot..' residues is too short a sequence for MakeHelix to work, sorry!')
end
print('MakeHelix may force parts of the protein into unfavorable conformations.')
print(' Color by relative score shows in red the regions in the least favorable conformations.\n ')
print('MakeHelix freezes parts of the protein and uses bands to change conformations.')
print('It also Wiggles the protein. Varying the Clashing Importance may help wiggling work better.')
print(' Changing the Clashing Importance can end the wiggling when it goes too long')
print(' and can extend the wiggling when it ends too soon.')
print(' ')
-- below are some default values
numa=1+0
numb=tot+0
way=4+0
--
qflag=0+0
while qflag==0 do
local ask=dialog.CreateDialog('~~/\\/\\/ MakeHelix /\\/\\/~~ ')
ask.Input1=dialog.AddSlider("Seg.A: ",numa,1,tot,0)
ask.Label1=dialog.AddLabel(" ^- Segment number where to START.\n ")
ask.Input2=dialog.AddSlider("Seg.B: ",numb,1,tot,0)
ask.Label2=dialog.AddLabel(" ^- Segment number where to END.\n ")
ask.Label2s=dialog.AddLabel("")
ask.Input3=dialog.AddSlider("Method: ",way,2,12,0)
--ask.Label3a=dialog.AddLabel("")
ask.Label3b=dialog.AddLabel("\n\n\n\n |- 2 = Beta Sheet\n |- 3 or 4 = Alpha Helix (3.6 seg per turn)\n |- 5 or 6 = '310' Helix (3 seg per turn)\n |- 7 or 8 = Pi Helix (4.4 seg per turn)\n |- 9 or 10 = Gamma Helix (5.1 seg per turn)\n '- 11 or 12 = '2.27 Ribbons' (2.2 seg per reversal)")
ask.Input3c=dialog.AddLabel(" ")
ask.Input3d=dialog.AddLabel(" ")
ask.Input3e=dialog.AddLabel(" ")
ask.Label3f=dialog.AddLabel(" ODDS (3,5,7,9,11) have Left-handed ends\nEVENS (4,6,8,10,12) have Right-handed ends")
ask.Input3g=dialog.AddLabel(" ")
ask.Input4=dialog.AddSlider("Quit-Flag: ",qflag,0,1,0)
ask.Label4=dialog.AddLabel(" \n 0 to Continue (loops), 1 to Quit\n ")
ask.OK = dialog.AddButton("OK", 1)
dialog.Show(ask)
numa=ask.Input1.value+0
numb=ask.Input2.value+0
way=ask.Input3.value+0
qflag=ask.Input4.value+0
nbold=band.GetCount() -- total number of bands on structure so far
print('Running MakeHelix using Start='..numa..' End='..numb..' Method='..way..' for '..tot..' residues (@'..os.date()..').')
print(' starting with '..nbold..' prior bands.')
print(' ')
if(qflag==0) then
--
-- below makes sure numa <= numb
if numa > numb then
print('Swap Seg.A and Seg.B to make A <= B hold:')
print(' At first Seg.A='..numa..' while Seg.B='..numb..',')
numa,numb=numb,numa
print(' but now Seg.A='..numa..' while Seg.B='..numb..'.')
print(' ')
end
if numb>numa+1 then
numy=numa+1
elseif numa>1 then
numy=numa-1
elseif numb<tot then
numy=numb+1
else
numy=1
print('Start='..numa..' End='..numb..' and Total='..tot..' do not allow a 3rd unique residue number for numy.')
end -- if numb
print('Start='..numa..' End='..numb..' and Total='..tot..' give numy='..numy..'.')
print(' ')
if(way>=3 and way<=12) then -- make a helix
nbmid,handstr=dohelix(numa,numb,numy,tot,way,nbold)
else -- way==2 -- make a beta-sheet
nbmid,handstr=dosheet(numa,numb,numy,tot,way,nbold)
end
-- below adds some anchor bands so frozen parts won't flop around so much
if numa>1 then
addtband(1,1)
end -- if numa
if numb<tot then
addtband(tot,-1)
end -- if numb
-- below puts anchors at strand breaks too
for num1=1,tot-1 do
num2=num1+1
if structure.GetDistance(num1,num2)>4 then -- break in strand between num1 and num2
if(num1-2>1 and (num1<numa or num1>numb)) then
addtband(num1,-1) -- addtband needs residues num1,num1-2,num1-1 to exist
end -- if num1
if(num2+2<tot and (num2<numa or num2>numb)) then
addtband(num2,1) -- addtband needs residues num2,num2+2,num2+1 to exist
end -- if num2
end -- if structure
end -- for num1
nbnew=band.GetCount() -- numbers of bands so far.
nbdiff=nbnew-nbold
print(' ')
print('Temporary bands '..(nbmid+1)..' to '..nbnew..' should help')
if(way>=3 and way<=12) then -- for helices
print(' establish '..handstr..'ness of helix and anchor rest of protein.')
else
print(' anchor the rest of the protein.')
end
print('Freezing everything not between residues '..numa..' and '..numb..'.')
selection.DeselectAll()
if numa-1>=1 then
selection.SelectRange(1,numa-1)
end
if numb+1<=tot then
selection.SelectRange(numb+1,tot)
end
freeze.FreezeSelected(true,true) -- freeze both backbone and sidechains
print('Start wiggling everything between residues '..numa..' and '..numb..' ('..os.date()..').')
print(' Wiggling can be slow if forcing residues into unfavorable conformations.')
print(' Vary the Clashing Importance under Behavior to make Wiggling end sooner or later.')
con=behavior.GetClashImportance()
print(' Beginning with Clashing Importance of '..con..'.')
structure.WiggleAll(2,true,true)
-- wiggle backbone and sidechains here to get everything comfortably in place
-- 1 iteration can still take a long time if bands can't get comfortable (default)
-- 5 iterations seems too much
con=behavior.GetClashImportance()
print(' Ending with Clashing Importance of '..con..'.')
freeze.UnfreezeAll()
print('Done wiggling now ('..os.date()..').')
print('Removing temporary bands '..(nbmid+1)..' to '..nbnew..' now.')
for bnum=nbnew,nbmid+1,-1 do
band.Delete(bnum)
end
print(' ')
--
end -- if qflag==0
nbnew=band.GetCount()
nbdiff=nbnew-nbold
print('MakeHelix added '..nbdiff..' (net) bands to the structure,')
print(' raising the total from '..nbold..' to '..nbnew..' by the end of the script.')
if qflag==0 then
if(way>=3 and way<=12) then -- for helices
print('Ideally, these '..nbdiff..' new bands (bands '..(nbold+1)..' to '..nbnew..') will maintain the '..handstr..' helical structure,')
print(' but if perturbed too much, the '..handstr..'ness of parts of the helix')
else
print('Ideally, these '..nbdiff..' new bands (bands '..(nbold+1)..' to '..nbnew..') will maintain the present secondary structure,')
print(' but if perturbed too much, parts of the secondary structure')
end
print(' (especially those in unfavorable conformations) can change.')
end -- if qflag
print(' ')
print('If you want to add new bands to a region that already contains bands,')
print(' now would be a good time to manually remove the old bands.')
print(' ')
end -- while qflag
print('END OF MAKEHELIX SCRIPT ('..os.date()..').')
end -- Main()
function dohelix(numa,numb,numy,tot,way,nbold)
local lperaa,hrad,nperturn,typestr,dab,mypi,hmypi
local w,handstr,hlen,xyz,num,numsh,num1,num2,nbmid,anums,rnums
if way==3 or way==4 then -- doing alpha-helix (some call it 3.6-13-helix)
--
-- For details see pp.26-7 of Lubert Stryer's "Biochemistry",
-- 3rd edition, 1988, W.H.Freeman & Company, ISBN=0-7167-1843-X.
-- For details see pp.1099-1100 of "Molecular Biology of the Cell",
-- 2nd edition, 1989, Garland Publishing, Inc., ISBN=0-8240-3695-6.
-- All N-H should point one way, all C=O should point other way
-- along the axis of the helix, here the x-axis.
-- C=O for residue n is bonded to N-H on residue n+4.
--
lperaa=1.5 -- lperaa is length helix per rise in residue number
-- 1.5 Angstroms per 100 degree rotation
-- 3.6 amino acids per turn (1 turn is 360 degrees)
-- (1.5 A/100 deg)*(360 deg/turn)*(1 turn/3.6 AA)=(1.5 A/AA)
hrad=2.3 -- radius of helix in Angstroms (1e-10 meters or 0.1 nm)
-- some give hrad=2.3, others give 2.5
nperturn=3.6 -- number of amino acids per turn in helix
typestr='Alpha Helix'
--
-- http://swift.cmbi.ru.nl/teach/aainfo/angles.shtml gives values like above and below
-- It calls lperaa the "Translation per residue"
-- and nperturn the "Residues per turn"
-- but does not give helix radii.
--
-- http://en.wikipedia.org/wiki/Protein_secondary_structure gives some info
-- http://biomedapps.curtin.edu.au/biochem/tutorials/prottute/helices.htm gives some info
-- http://www.answers.com/topic/gamma-helix gives gamma helix
-- http://www.answers.com/topic/pi-helix for pi-helix
-- http://www.answers.com/topic/3-10-helix
--
elseif way==5 or way==6 then -- doing 310-helix
lperaa=2.0
hrad=1.9
nperturn=3.0
typestr='310 Helix'
elseif way==7 or way==8 then -- doing pi-helix (one says Xi-helix, another says 4.4-16-helix)
-- below info comes from where?
lperaa=1.1 -- some give 1.1, others 1.15
hrad=2.8
nperturn=4.4 -- some give 4.3, others 4.4, others 4.1
typestr='Pi Helix'
elseif way==9 or way==10 then -- doing gamma-helix -- should be 5.03 A/turn
lperaa=0.99 -- 0.99 A/AA
hrad=3.2 -- calc'd from lperaa and nperturn 5/26/12
nperturn=5.1 -- 5.1 AA/turn -- (0.99 A/AA)x(5.1 AA/turn)=(5.03 A/turn)
typestr='Gamma Helix'
-- Seems like helices obey formulas like below:
-- N=#AA/turn=nperturn
-- L=length/AA=lperaa=distance along x-axis
-- R=radius=hrad
-- Imagine helix as a cylinder of radius R.
-- When viewed from one end (along the cylinder's axis, the x-axis here, where y=z=0)
-- the vector from the cylinder axis to each alpha-carbon is like a spoke on a wheel.
-- These spokes have length R. Their vectors have x-components of zero.
-- When viewed from one end, the angle between each spoke is theta=(2 pi)/N radians or 360/N degrees.
-- Imagine that all these spokes and their alpha-carbon endpoints project onto the yz plane.
-- Let W be the distance in the yz plane between the projections of 2 adjacent alpha-carbons.
-- These will make many triangles with sides of length W,R,R and the angle theta between the sides of length R.
-- These triangles give W = 2 R sin(theta/2) and w^2 = R^2 + R^2 - 2 R^2 cos(theta) = 2 R^2 (1-cos(theta))
-- so that W = sqrt(2) R sqrt(1-cos(theta)) holds.
-- Also, the actual distance between 2 adjacent alpha-carbons will be l=sqrt(L^2+W^2).
-- In general, this distance l should be fixed.
-- l=3.83 Angstroms is fairly consistent with many of the L=lperaa, R=hrad, N=nperturn values listed here.
-- l=3.79 Angstroms would match the beta-sheet structure used here.
elseif way==11 or way==12 then -- doing 2.2-7-helix
lperaa=0
hrad=0
nperturn=2.2
typestr='2.27 Ribbon'
end -- if way
dab=structure.GetDistance(numa,numb)
mypi=math.pi
hmypi=mypi/2.0
w=2.0*mypi/nperturn -- how many radians the angle increases per amino acid
handstr='right-handed'
if way==3 or way==5 or way==7 or way==9 or way==11 then -- reverse sign of w for left-handed helices
w= -w
handstr='left-handed'
end -- if way
print('Method='..way..' will make a '..handstr..' '..typestr..' between residues '..numa..' and '..numb..' below:')
print(' ')
hlen=(numb-numa)*lperaa -- overall length of alpha-helix
--
-- next builds helix in x,y,z coordinates in 2-dimensional table xyz
xyz={} -- stores x,y,z coordinates
anums={} -- stores atom numbers
rnums={} -- stores residue (AA) numbers
for num=numa,numb do
x=0.5*(dab-hlen)+(num-numa)*lperaa -- helix translates along x axis from numa to numb
y=hrad*math.cos(w*(num-numa)) -- circular part of helix is in yz plane
z=hrad*math.sin(w*(num-numa))
aa=structure.GetAminoAcid(num)
print('Assuming CA of '..aa..num..' is at (x,y,z)=('..x..','..y..','..z..').')
xyz[num]={x,y,z}
rnums[num]=num+0
anums[num]=0+0 -- 0 is for default which should be alpha-carbon CA (2 might work here too)
structure.SetSecondaryStructure(num,'H') -- H is for helix
end -- for num
print(' ')
for numsh=2,4 do -- which bands are really needed to maintain the structure?
--
-- ideally for alpha helices:
-- CA's on residues n and n+3 are very close
-- CA's on n and n+4 are very close
--
for num1=numa,numb-numsh do
num2=num1+numsh
setdist(rnums,anums,xyz,num1,num2)
end -- for num1
print(' ')
end -- for numsh
nbmid=band.GetCount() -- number of bands so far.
-- nbmid is the last band number for residue to residue constraints.
nbdiff=nbmid-nbold
print('MakeHelix has added '..nbdiff..' bands ('..(nbold+1)..' to '..nbmid..') to the structure for internal constraints.')
print(' ')
-- nbmid+1 will be the first band number for residue to space constraints.
-- nbmid+1 and on seem needed to enforce the handedness of a helix
print('Adding '..(numb-numa+1)..' temporary bands to enforce '..handstr..'ness.')
for num=numa,numb do
bandadd(numa,numb,numy,num,rnums,anums,xyz)
end -- for num
return nbmid,handstr
end -- dohelix()
function dosheet(numa,numb,numy,tot,way,nbold)
local nbmid,handstr,xyz,num,num1,num2,dab,anums,rnums,x,y,z,anum,ii,slen,lperaa
local h,hh,xsh,ysh,r1,r3,r4,ang1,ang3,ang4,cf,npts
-- make a beta-sheet based on distances and angles in
-- http://n0b3l1a.blogspot.com/2009/01/20-amino-acids.html
handstr=''
nbmid=0+0
print('Method='..way..' will make a Beta-Sheet between residues '..numa..' and '..numb..' below:')
print(' ')
dab=structure.GetDistance(numa,numb)
lperaa=3.5 -- distance in x direction between adjacent CA's (width in xy plane)
h=1.460959255 -- distance in y direction between adjacent CA's (height in xy plane)
hh=h/2.0 -- half the height
slen=(numb-numa)*lperaa -- overall length of beta-sheet in x-direction
xyz={}
anums={}
rnums={}
cf=math.pi/180.0 -- converts degrees into radians
--
-- below sets up relative positions within each peptide bond / amide plane
-- only x,y are needed since all atoms are in a plane
-- CA will always be at 0,0
--
-- N is next
r1=1.46
ang1=cf*37.47590055
x1= -r1*math.cos(ang1)
y1=r1*math.sin(ang1)
--
-- C in C=O bond is next
r3=1.51
ang3=cf*43.47590055
x3=r3*math.cos(ang3)
y3=r3*math.sin(ang3)
--
-- O in C=O bond is next
r4=1.24
ang4=ang3+cf*123.5
x4=x3+r4*math.cos(ang4)
y4=y3+r4*math.sin(ang4)
--
-- next lay out all the peptide bonds / amide planes side by side
npts=0+0
for num=numa,numb do
x=0.5*(dab-slen)+(num-numa)*lperaa -- translate along x axis from numa to numb
-- x is position of an alpha-carbon CA (atom 2 or 0 for each amino acid)
z=0+0 -- all of beta-sheet is in xy plane so set z to 0
for anum=1,4 do -- will position atoms 1,2,3,4 for each amino acid
ii=(num-numa)*4+anum
anums[ii]=anum+0
rnums[ii]=num+0
aa=structure.GetAminoAcid(num)
-- in below, xsh is how far in x-direction from alpha-carbon
-- and ysh is how far in y-direction from alpha-carbon position
-- as if alpha-carbon is at origin in xsh,ysh system
if anum==1 then -- N
xsh=x1
ysh=y1
astr='N'
elseif anum==2 then -- CA
xsh=0+0
ysh=0+0
astr='CA'
elseif anum==3 then -- C in C=O bond
xsh=x3
ysh=y3
astr='C'
else -- anum==4 -- O in C=O
xsh=x4
ysh=y4
astr='O'
end -- if anum
if isodd(num)==1 then
y=ysh-hh -- have odd ones point up
else
y=hh-ysh -- have even ones point down
end -- if isodd
xyz[ii]={x+xsh,y,z}
print('('..ii..') Placing '..astr..' (atom '..anum..') of '..aa..num..' at (x,y,z)=('..(x+xsh)..','..y..','..z..').')
structure.SetSecondaryStructure(num,'E') -- E is for sheet
npts=npts+1
end -- for anum
end -- for num
print(' ')
-- which bands are really needed to maintain the structure?
--
-- ideally for beta-sheets:
-- CA's on residues n and n+2 are 7.0 apart
-- C's on residues n and n+2 are 7.0 apart
-- O's on residues n and n+2 are 7.0 apart
-- N's on residues n and n+2 are 7.0 apart
-- N on n and C on n+2 are farther than 7.0 apart
-- C on n and N on n+2 are closer than 7.0 apart
for num=numa,numb-2 do
-- for anum=1,4 do -- do bands between all 4 atom types
for anum=2,2 do -- just do bands between CA's
ii=(num-numa)*4+anum
setdist(rnums,anums,xyz,ii,ii+8) -- band same atoms on residues n and n+2
end -- for anum
ii=(num-numa)*4
setdist(rnums,anums,xyz,ii+1,ii+11) -- band N on n to C on n+2
setdist(rnums,anums,xyz,ii+3,ii+9) -- band C on n to N on n+2
end -- for num
print(' ')
nbmid=band.GetCount() -- number of bands so far.
-- nbmid is the last band number for residue to residue constraints.
nbdiff=nbmid-nbold
print('MakeHelix has added '..nbdiff..' bands ('..(nbold+1)..' to '..nbmid..') to the structure for internal constraints.')
print(' ')
-- nbmid+1 will be the first band number for residue to space constraints.
print('Adding '..npts..' temporary bands to enforce positions.')
for num=1,npts do
bandadd(numa,numb,numy,num,rnums,anums,xyz)
end -- for num
return nbmid,handstr
end -- dosheet()
function addtband(num,dnum)
-- dnum=+1 uses num,num+2,num+1
-- dnum=-1 uses num,num-2,num-1
local anchorlen=0.0001 -- length of each anchor
local hmypi=math.pi/2.0
print('Adding a temporary band to residue '..num..' to anchor its end.')
local bnum=band.Add(num,num+2*dnum,num+dnum,anchorlen,hmypi,hmypi) -- along local -x axis
band.SetGoalLength(bnum,0)
band.SetStrength(bnum,10)
end -- addtband()
function bandadd(numa,numb,numy,num,rnums,anums,xyz)
local x,y,z=xyz[num][1],xyz[num][2],xyz[num][3]
local rho=math.sqrt(x^2+y^2+z^2)
local theta=math.acos(z/rho)
local phi=math.atan2(y,x)
local anum=anums[num]
local rnum=rnums[num]
if phi<0.0 then
phi=phi+2.0*math.pi
end -- if phi
local bnum1=band.Add(numa,numb,numy,rho,theta,phi)
-- previous line uses numa's position as origin, x direction is toward numb,
-- y direction is perpendicular to x direction and gives positive y on the side toward residue numy
local bnum2=band.AddToBandEndpoint(rnum,bnum1,anum)
band.SetStrength(bnum2,10) -- 10 is max, 1 is default
band.SetGoalLength(bnum2,0) -- 3.5 is default
band.Delete(bnum1)
end -- bandadd()
function setdist(rnums,anums,xyz,num1,num2)
local dist=getdist(xyz[num1],xyz[num2],3)
local anum1=anums[num1]
local anum2=anums[num2]
local rnum1=rnums[num1]
local rnum2=rnums[num2]
local bnum=band.AddBetweenSegments(rnum1,rnum2,anum1,anum2)
band.SetGoalLength(bnum,dist) -- 3.5 is default
band.SetStrength(bnum,10) -- 10 is max, 1 is default
local aa1=structure.GetAminoAcid(rnum1)
local aa2=structure.GetAminoAcid(rnum2)
print('Setting distance between atom '..anum1..' of '..aa1..rnum1..' and atom '..anum2..' of '..aa2..rnum2..' to '..dist..'.')
end -- setdist()
function getdist(p1,p2,dim)
-- finds the distance between two vectors of length dim
local res=0.0, ii
for ii=1,dim do
res=res+(p1[ii]-p2[ii])^2
end -- for ii
res=math.sqrt(res)
return res
end -- getdist()
-- below returns 1 if num is odd and 0 if num is even
function isodd(num)
local res=0
if math.fmod(num,2)==1 then
-- if num/2 has remainder 1
res=1
end -- if num
return res
end -- isodd()
-- below returns 1 if num is even and 0 if num is odd
function iseven(num)
local res=0
if math.fmod(num,2)==0 then
-- if num/2 has remainder 0
res=1
end -- if num
return res
end -- iseven()
Main()
