Icon representing a recipe

Recipe: makehelix

created by jeff101

Profile


Name
makehelix
ID
42111
Shared with
Public
Parent
None
Children
Created on
June 03, 2012 at 03:52 AM UTC
Updated on
June 03, 2012 at 03:52 AM UTC
Description

makes beta-sheets and different kinds of helices in a protein

Best for


Code


function Main() -- -- makehelix 5/15/12 by jeff101 -- tot=structure.GetCount() -- find total # of residues if tot<3 then print('tot='..tot..' is too short a sequence for makehelix to work.') 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("numa: ",numa,1,tot,0) ask.Label1=dialog.AddLabel(" \nnuma is amino acid number where helix starts.\n ") ask.Input2=dialog.AddSlider("numb: ",numb,1,tot,0) ask.Label2=dialog.AddLabel(" \nnumb is amino acid number where helix ends.\n ") ask.Input3=dialog.AddSlider("way: ",way,2,10,0) ask.Label3a=dialog.AddLabel(" \nway 2 does beta-sheets. ways 3-10 do helices.") ask.Label3b=dialog.AddLabel("ways 3,5,7,9 do left-handed helices.") ask.Label3c=dialog.AddLabel("ways 4,6,8,10 do right-handed helices.") ask.Label3d=dialog.AddLabel("ways 3,4 do alpha, ways 5,6 do 3-10.") ask.Label3e=dialog.AddLabel("ways 7,8 do pi (4.4-16), ways 9,10 do gamma.\n ") ask.Input4=dialog.AddSlider("qflag: ",qflag,0,1,0) ask.Label4=dialog.AddLabel(" \nqflag is 0 to continue, 1 to quit now.\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 with numa='..numa..' numb='..numb..' way='..way..' qflag='..qflag..' for '..tot..' residues ('..os.date()..').') print(' starting with '..nbold..' bands.') print(' ') if(qflag==0) then -- -- below makes sure numa <= numb if numa > numb then print('Swap numa and numb to make numa <= numb hold:') print(' At first numa='..numa..' while numb='..numb..',') numa,numb=numb,numa print(' but now numa='..numa..' while numb='..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('numa='..numa..' numb='..numb..' and tot='..tot..' do not allow a 3rd unique residue number for numy.') end -- if numb print('numa='..numa..' numb='..numb..' and tot='..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(1,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 -- 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 3-10-helix lperaa=2.0 hrad=1.9 nperturn=3.0 typestr='3-10-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.2-7-helix' 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('way='..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('way='..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()

Comments


jeff101 Lv 1

This recipe asks the user for two residue numbers and a secondary-structure type. It then adds bands to all residues between the two residue numbers to make that region of the protein change to the secondary-structure's conformation. It freezes all other parts of the protein and adds some anchor bands if the protein is more than one strand. Then it wiggles the protein. Ideally, the wiggling ends quickly, and then the recipe repeats, asking the user for new residue numbers and a new secondary-structure type.

When using this recipe, it often helps to vary the Clashing Importance to change the amount of time spent wiggling. You can also remove bands when the recipe pauses between user inputs. It also helps to run this recipe with residues colored by relative score. Then you can tell which regions are in the most unfavorable conformations (they will be red rather than green).

Right now, makehelix is coded to make beta-sheets and left- or right-handed alpha, 3-10, pi, or gamma helices. I would like to expand it to include 2.2-7 helices and different kinds of turns. If you know of references listing details (distances, angles) about these structures, please e-mail me about them.

Another feature that would help is better detection of strand breaks in the protein (like in puzzle 547). Right now, makehelix uses a distance cutoff to detect breaks. If you know of other ways to detect the ends of individual strands in the protein, please e-mail me about them.

If you have other suggestions for improvements to makehelix, please e-mail them to me. Most things can be improved, and the more complicated a piece of software is, the more likely it is to have bugs.

Bruno Kestemont Lv 1

I discover this recipe too late. I think than even after the new SS tool, it can be useful by allowing us to vary the helix structure (right handed left handed etc).
Thanks !