Code
-- HBonds1 begun 4/17/2021 by jeff101.
-- Renamed to ListNetworks1 4/18/2021.
-- Last updated 4/18/2021 1246pm noon.
--
-- Below lists all global variables.
aas={} -- will hold 3-letter codes
-- and full names for each aa
aas.a={'ala','alanine'}
aas.c={'cys','cysteine'}
aas.d={'asp','aspartate'}
aas.e={'glu','glutamate'}
aas.f={'phe','phenylalanine'}
aas.g={'gly','glycine'}
aas.h={'his','histidine'}
aas.i={'ile','isoleucine'}
aas.k={'lys','lysine'}
aas.l={'leu','leucine'}
aas.m={'met','methionine'}
aas.n={'asn','asparagine'}
aas.p={'pro','proline'}
aas.q={'gln','glutamine'}
aas.r={'arg','arginine'}
aas.s={'ser','serine'}
aas.t={'thr','threonine'}
aas.v={'val','valine'}
aas.w={'trp','tryptophan'}
aas.x={'unk','unknown'}
aas.y={'tyr','tyrosine'}
bbstr='ncco' -- atoms 1-4 for most AA's.
scstr={} -- will hold single-letter atom names
-- for atom #'s 5 and up (not including H's)
-- on non-terminal amino acids.
-- atoms 1-4 for all AA's should be ncco
-- (atom 1 is nitrogen, 2 is the alpha-carbon)
-- (3 is the carbonyl carbon in C=O, 4 is oxygen).
-- atom 5 for C-terminal (highest seg #) is O,
-- so the ones listed below become 6 and up.
-- seg #1 (N-terminal) has 2 extra H's, but
-- they don't matter here.
scstr.a='c'
scstr.c='cs'
scstr.d='ccoo'
scstr.e='cccoo'
scstr.f='ccccccc'
scstr.g='' -- has backbone + atom 5 = h
scstr.h='ccnccn' -- correct for both his types?
scstr.i='cccc'
scstr.k='ccccn'
scstr.l='cccc'
scstr.m='ccsc'
scstr.n='ccon'
scstr.p='ccc'
scstr.q='cccon'
scstr.r='cccncnn'
scstr.s='co'
scstr.t='coc'
scstr.v='ccc'
scstr.w='ccccnccccc'
scstr.x='?' -- for unknown sc's. not really used
scstr.y='ccccccco'
-- Above lists all global variables.
--
-- Below lists all functions.
-- Main() below was inspired by the sample code given at:
-- https://foldit.fandom.com/wiki/Foldit_Lua_Function_structure.GetHBonds
-- For more details about other Foldit Lua Functions see:
-- https://foldit.fandom.com/wiki/Foldit_Lua_Functions
--
function Main()
local nmon,nseg,aa1,ss1,aa2,ss2
local seg1,seg2,seg2o,atom1,atom2,mon1,mon2
local isdis,score1,score2,dist,width
local atype1,atype2,btype1,btype2,bway,listit
local fnams,fison,nf,ff,scoreA,scoreB
local mssg,bonds,nbonds,bnum,tots
-- nmon is number of monomers
nmon=structure.GetSymCount()
-- nseg is number of segments/monomer
nseg=structure.GetCount()
bonds=structure.GetHBonds()
nbonds=(#bonds) -- # of bonds
--
print('Running some tests:')
print(' round2( 3.445)='..round2( 3.445)..'.')
print(' round2( 3.444)='..round2( 3.444)..'.')
print(' round2(-3.444)='..round2(-3.444)..'.')
print(' round2(-3.445)='..round2(-3.445)..'.')
print(' round3( 3.4475)='..round3( 3.4475)..'.')
print(' round3( 3.4474)='..round3( 3.4474)..'.')
print(' round3(-3.4474)='..round3(-3.4474)..'.')
print(' round3(-3.4475)='..round3(-3.4475)..'.')
print(' ')
--
-- below records settings of filters
-- then enables them all. can it record
-- which filters are showing too?
fnams=filter.GetNames() -- get names of filters
nf=(#fnams) -- # of filters
fison={} -- 1 for enabled filters, 0 for disabled ones
for ff=1,nf do
if filter.IsEnabled(fnams[ff]) then
fison[ff]=1
else -- filter is disabled
fison[ff]=0
end -- if filter
end -- for ff
filter.EnableAll()
-- below gets score with all bonuses
scoreA=current.GetEnergyScore()
-- below disables all filters
filter.DisableAll()
-- below gets score with no bonuses
scoreB=current.GetEnergyScore()
-- below restores to original filter settings.
-- can it restore which filters were shown too?
-- 4/18/2021 seems disabling a filter quits
-- displaying it, but re-enabling a filter
-- does not re-display it.
for ff=1,nf do
if fison[ff]==1 then
filter.Enable(fnams[ff])
end -- if fison
end -- for ff
--
print(string.format('Score=%.3f Bonus=%.3f=%.3f HBN=%.3f has',
trunc3(scoreA),
round3(scoreA-scoreB),round3(filter.GetBonusTotal()),
round3(filter.GetBonus('Hydrogen Bond Network'))))
print(string.format(' %d monomers (original + %d copies), each w/%d segments.',
nmon+1,nmon,nseg))
print('There are '..nbonds..' bonds involving copy 0')
print(' (the original monomer), as listed below.')
print('Width (wid) = 0.25 A (bad) to 0.75 A (good).')
print(' {')
tots={}
for bway=1,3 do
tots[bway]=0
if bway==1 then
print('Ones between 2 backbone atoms on the same monomer (bway=1):')
elseif bway==2 then
print('Ones involving 1-2 sidechain atoms on the same monomer (bway=2):')
else
print('Ones between atoms on 2 different monomers (bway==3):')
print('Some bonds between monomers are repeats by symmetry.')
end
print(' {')
for bnum=1,nbonds do
seg1=bonds[bnum]["res1"]
seg2o=bonds[bnum]["res2"]
-- 4/17/2021 It seems like seg1 <= seg2o always holds.
atom1=bonds[bnum]["atom1"]
atom2=bonds[bnum]["atom2"]
width=bonds[bnum]["width"] -- 0.25 A (bad) to 0.75 A (good)
mon1=0 -- This should always hold
-- since the original monomer is copy 0.
mon2,seg2=getmonsig(seg2o,nseg)
-- If mon2 > 0, seg2 < seg1 <= seg2o can hold.
-- 4/17/2021 Seems can get seg1 seg2 in any
-- order vs each other, and some pairs
-- can even be repeats of each other.
aa1=structure.GetAminoAcid(seg1)
aa2=structure.GetAminoAcid(seg2)
atype1,btype1=gettype(nseg,seg1,atom1,aa1)
atype2,btype2=gettype(nseg,seg2,atom2,aa2)
listit=0
if mon1~=mon2 then
-- bond is between atoms on 2 diff monomers
if bway==3 then
listit=1
end -- if bway
elseif btype1=='sc' or btype2=='sc' then
-- bond is between 1 or more sc atoms on the same monomer
if bway==2 then
listit=1
end -- if bway
else -- btype1=='bb' and btype2=='bb'
-- bond is between 2 bb atoms on the same monomer
if bway==1 then
listit=1
end -- if bway
end -- if mon1
if listit==1 then
tots[bway]=tots[bway]+1
bonds[bnum].bway=bway
ss1=structure.GetSecondaryStructure(seg1)
ss2=structure.GetSecondaryStructure(seg2)
dist=getdist(seg1,atom1,mon2,seg2,atom2)
bonds[bnum].dist=dist
bonds[bnum].mons={mon1,mon2}
bonds[bnum].segs={seg1,seg2}
bonds[bnum].sss={ss1,ss2}
bonds[bnum].aas={aa1,aa2}
bonds[bnum].atoms={atom1,atom2}
bonds[bnum].atypes={atype1,atype2}
bonds[bnum].btypes={btype1,btype2}
mssg=string.format(' (%d) Bond %d ',
tots[bway],bnum)
isdis=bonds[bnum]["bond_type"]
-- 0 for hbonds, 1 for disulfides
bonds[bnum].isdis=isdis
if isdis==0 then -- for hbonds
mssg=(mssg..'(hbond)')
score1=current.GetSegmentEnergySubscore(seg1,'Bonding')
score2=current.GetSegmentEnergySubscore(seg2,'Bonding')
else -- for disulfides
mssg=(mssg..'(disulfide)')
score1=current.GetSegmentEnergySubscore(seg1,'Disulfides')
score2=current.GetSegmentEnergySubscore(seg2,'Disulfides')
end -- if isdis
bonds[bnum].scores={score1,score2}
print(string.format('%s %d%s%d%s%d%s-%d%s%d%s%d%s %s-%s len %.2f A wid %.2f A w/scores %.3f %.3f.',
mssg,mon1,ss1,seg1,aas[aa1][1],atom1,atype1,
mon2,ss2,seg2,aas[aa2][1],atom2,atype2,
btype1,btype2,round2(dist),round2(width),
round3(score1),round3(score2)))
end -- if listit
end -- for bnum
print(' }')
end -- for bway
print(' }')
print(string.format('%d + %d + %d = %d should equal %d.',
tots[1],tots[2],tots[3],tots[1]+tots[2]+tots[3],nbonds))
print(' ')
getnets(bonds,nbonds,nseg)
print(' ')
print('All done.')
end -- Main()
-- Below finds monomer # mon2
-- and segment # seg2 given
-- a raw segment # seg2o
-- and nseg, the # of segments
-- per monomer.
-- For a monomer, mon2 = 0 (orig monomer).
-- For a dimer, mon2 = 0 to 1 (purple is 1).
-- For a trimer, mon2 = 0 to 2 (teal is 2).
-- For a tetramer, mon2 = 0 to 3 (brown is 3).
-- For a pentamer, mon2 = 0 to 4.
--
function getmonsig(seg2o,nseg)
local mon2,seg2,nsegp1
nsegp1=nseg+1 -- includes 1 pseudo-residue
-- at the end of each monomer.
mon2=math.ceil(seg2o/nsegp1)-1
-- seg2o=0*nsegp1+1 to 1*nsegp1 gives mon2=0
-- seg2o=1*nsegp1+1 to 2*nsegp1 gives mon2=1
-- seg2o=2*nsegp1+1 to 3*nsegp1 gives mon2=2
-- seg2o=3*nsegp1+1 to 4*nsegp1 gives mon2=3
seg2=seg2o-mon2*nsegp1
return mon2,seg2
end -- getmonsig()
-- Below finds distance dist between
-- atom atom1 of segment seg1 of monomer 0 and
-- atom atom2 of segment seg2 of monomer mon2.
-- It returns -1 if no distance is found.
--
function getdist(seg1,atom1,mon2,seg2,atom2)
local dist= -1
local btmp,bold
bold=band.GetCount() -- # bands at start of getdist
btmp=band.AddBetweenSegments(seg1,seg2,atom1,atom2,mon2)
-- In above, seg1 atom1 is on monomer mon1 = 0
-- while seg2 atom2 is on monomer mon2.
if btmp>bold then -- is btmp a new band?
dist=band.GetLength(btmp)
band.Delete(btmp)
end -- if btmp
return dist
end -- getdist()
-- Below returns btype & atype
-- where btype='bb' for backbone atoms
-- and btype='sc' for sidechain atoms
-- while atype='o' for oxygen atoms
-- and atype='n' for nitrogen atoms.
--
function gettype(nseg,seg,atom,aa)
local atype,btype
-- below will crash if the atom is an H
if aa=='x' then -- AA is unknown
atype='?'
btype='?'
else -- AA is one of the 20 standard AA's
if atom<=4 then -- it's a backbone atom
atype=getchar(bbstr,atom)
btype='bb'
elseif seg==nseg then -- doing C-terminal
if atom==5 then
-- atom 5 is an extra backbone O
atype='o'
btype='bb'
else -- it's a sidechain atom
atype=getchar(scstr[aa],atom-5)
btype='sc'
end -- if atom
else -- it's a normal sidechain's atom
atype=getchar(scstr[aa],atom-4)
btype='sc'
end -- if atom
end -- if aa
return atype,btype
end -- gettype()
-- Below returns the num-th
-- character in the string str.
--
function getchar(str,num)
return string.sub(str,num,num)
--
-- below might work too
-- return string.char(string.byte(str,num))
end -- getchar()
-- Below rounds inval to 2 digits
-- after the decimal point so that
-- 0.035-0.039 rounds up to 0.04 and
-- 0.031-0.034 rounds down to 0.03.
--
function round2(inval)
return roundn(inval,2)
end -- round2()
-- Below rounds inval to 3 digits
-- after the decimal point so that
-- 0.0075-0.0079 rounds up to 0.008 and
-- 0.0071-0.0074 rounds down to 0.007.
--
function round3(inval)
return roundn(inval,3)
end -- round3()
-- Below rounds inval to n digits
-- after the decimal point.
--
function roundn(inval,n)
local tmpval,outval
local tfac=10^n
tmpval=tfac*inval
if tmpval>0 then
outval=math.floor(tmpval)
if tmpval-outval>=0.5 then
outval=outval+1
end -- if tmpval
else -- below handles negative numbers
outval=math.ceil(tmpval)
if outval-tmpval>=0.5 then
outval=outval-1
end -- if outval
end -- if tmpval
return (outval/tfac)
end -- roundn()
-- Below is meant to round total
-- scores like Foldit does.
-- Does it work right for
-- negative scores?
--
function trunc3(inval)
return truncn(inval,3)
end -- trunc3()
-- Below truncates the number
-- inval to the nth digit after
-- the decimal place.
--
function truncn(inval,n)
local tfac=10^n
return (math.floor(inval*tfac)/tfac)
end -- truncn()
-- Below is to break up the groups
-- of bands into networks.
--
function getnets(bonds,nbonds,nseg)
local nets={}
local bondnets={} -- which nets each bond is in
local segnets={} -- which nets each segment is in
local bnum,bnum1
local nnets,net1,segn,seg,seg1,seg2
local bway,isold,mssg,jmax,ii
nnets=0
for segn=1,nseg do
segnets[segn]={}
end -- for segn
for bnum=1,nbonds do
bondnets[bnum]={}
if bonds[bnum].isdis==0 then
-- do below for hbonds not disulfide bonds
bway=bonds[bnum].bway
if bway>1 then
-- hbond is not between 2 bb atoms on the same monomer
seg1=bonds[bnum].segs[1]
seg2=bonds[bnum].segs[2]
isold=0 -- means seg1,seg2 aren't in networks 1 to nnets.
-- 1 means just found a network for bnum seg1 seg2.
-- 2 means already put bnum seg1 seg2 into a network.
for net1=1,nnets do
if isold==0 then
for segn=1,(#(nets[net1].segs)) do
if isold==0 then
seg=nets[net1].segs[segn]
if seg==seg1 or seg==seg2 then
-- will add bnum seg1 seg2 to network net1
isold=1
end -- if seg
end -- if isold
end -- for segn
end -- if isold
if isold==1 then
-- add bnum seg1 seg2 to network net1
segnets[seg1]=add2vec(segnets[seg1],net1)
segnets[seg2]=add2vec(segnets[seg2],net1)
bondnets[bnum]=add2vec(bondnets[bnum],net1)
nets[net1].bonds=add2vec(nets[net1].bonds,bnum)
nets[net1].segs=add2vec(nets[net1].segs,seg1)
nets[net1].segs=add2vec(nets[net1].segs,seg2)
-- Larger bway is better, and
-- overall bway is decided by the best bway so far.
if bway > nets[net1].bway then
nets[net1].bway=bway
end -- if bway
isold=2 -- won't put bnum seg1 seg2 in another network
end -- if isold
end -- for net1
if isold==0 then
-- add a new network to end of nets,
-- and put bnum seg1 seg2 in it.
nnets=nnets+1
nets[nnets]={}
nets[nnets].bonds={}
nets[nnets].segs={}
segnets[seg1]=add2vec(segnets[seg1],nnets)
segnets[seg2]=add2vec(segnets[seg2],nnets)
bondnets[bnum]=add2vec(bondnets[bnum],nnets)
nets[nnets].bonds=add2vec(nets[nnets].bonds,bnum)
nets[nnets].segs=add2vec(nets[nnets].segs,seg1)
nets[nnets].segs=add2vec(nets[nnets].segs,seg2)
nets[nnets].bway=bway
end -- if isold
--
-- Below checks if any nets contain the same
-- bonds or segments. If they do, it merges
-- the nets into one.
for bnum1=1,bnum do
if (#(bondnets[bnum1]))>1 then
-- bond bnum1 is in more than one network
mssg,jmax=listvec(bondnets[bnum1])
print(string.format('Bond %d is in %d networks%s.',
bnum1,jmax,mssg))
print(' so merge these networks below:')
nets,nnets=mergenets(nets,bondnets[bnum1])
-- below updates segnets & bondnets
segnets=newnets(nets,nnets,nseg,1)
bondnets=newnets(nets,nnets,bnum,0)
print(' ')
end -- if
end -- for bnum1
for segn=1,nseg do
if (#(segnets[segn]))>1 then
-- segment segn is in more than one network
mssg,jmax=listvec(segnets[segn])
print(string.format('Segment %d is in %d networks%s,',
segn,jmax,mssg))
print(' so merge these networks below:')
nets,nnets=mergenets(nets,segnets[segn])
-- below updates segnets & bondnets
segnets=newnets(nets,nnets,nseg,1)
bondnets=newnets(nets,nnets,bnum,0)
print(' ')
end -- if
end -- for segn
end -- if bway
end -- if bonds[bnum].isdis
end -- for bnum
--
-- below outputs the nets found above
print('Found '..nnets..' different networks:')
for netn=1,nnets do
print(string.format(' Network %d w/bway=%d contains',
netn,nets[netn].bway))
mssg,jmax=listvec(nets[netn].bonds)
print(string.format(' %d bonds%s and',jmax,mssg))
mssg,jmax=listvec(nets[netn].segs)
print(string.format(' %d segments%s or',jmax,mssg))
mssg,jmax=listvecNaas(nets[netn].segs)
print(string.format(' %s.',mssg))
end -- for netn
end -- getnets()
-- Below find segnets based on the
-- latest contents of nets.
-- sflag=1 does for segments,
-- sflag=0 does for bonds.
--
function newnets(nets,nnets,nseg,sflag)
local segnets={}
local seg1,net1,vlen,ii
for seg1=1,nseg do
segnets[seg1]={}
end -- for seg1
for net1=1,nnets do
if sflag==1 then -- do for segments
vlen=(#(nets[net1].segs))
else -- do for bonds
vlen=(#(nets[net1].bonds))
end -- if sflag
for ii=1,vlen do --
if sflag==1 then -- do for segments
seg1=nets[net1].segs[ii]
else -- do for bonds
seg1=nets[net1].bonds[ii]
end -- if sflag
-- below adds network net1 to segnets[seg1],
-- the list of networks containing segment seg1.
segnets[seg1]=add2vec(segnets[seg1],net1)
end -- for ii
end -- for net1
return segnets
end -- newnets()
-- Below takes all the nets listed in nlist
-- and merges them to make a new version of nets.
-- If bway differs among the nets to be merged,
-- the merged net will have the largest bway.
--
function mergenets(nets,nlist)
local mnets={} -- merged version of nets
local inds={} -- mnets[inds[ii]] will include nets[ii]
local ii,jj,lo,hi,mssg,jmax,net1
net1=nlist[1]
hi=0
for ii=2,(#nlist) do
lo=hi+1
hi=nlist[ii]
for jj=lo,hi-1 do
inds[jj]=jj-(ii-2)
end -- for jj
inds[hi]=net1
end -- for ii
lo=hi+1
hi=(#nets)
for jj=lo,hi do
inds[jj]=jj-((#nlist)-1)
end -- for jj
mssg,jmax=listvec(nlist)
print(string.format('For %d nets, nlist%s',(#nets),mssg))
mssg,jmax=listvec(inds)
print(string.format(' gives inds%s.',mssg))
hi=(#nets)-((#nlist)-1) -- hi is now desired length for mnets
for ii=1,hi do
mnets[ii]={}
mnets[ii].bonds={}
mnets[ii].segs={}
mnets[ii].bway=0
end -- for ii
for ii=1,(#nets) do
for jj=1,(#(nets[ii].bonds)) do
mnets[inds[ii]].bonds=add2vec(mnets[inds[ii]].bonds,nets[ii].bonds[jj])
end -- for jj
for jj=1,(#(nets[ii].segs)) do
mnets[inds[ii]].segs=add2vec(mnets[inds[ii]].segs,nets[ii].segs[jj])
end -- for jj
if mnets[inds[ii]].bway < nets[ii].bway then
mnets[inds[ii]].bway = nets[ii].bway
end -- if mnets
end -- for ii
return mnets,(#mnets)
end -- mergenets()
-- Below looks through the elements of ovec for val
-- and if it doesn't find val, it puts val in numerical
-- order to make nvec, which it returns.
--
function add2vec(ovec,val)
local nvec,vlen,ii,jj,qflag
qflag=0
vlen=(#ovec)
nvec={}
for ii=1,vlen do
if qflag==0 then -- don't quit yet
if val<ovec[ii] then
-- insert val into nvec
-- so nvec[ii]=val.
for jj=1,ii-1 do
nvec[jj]=ovec[jj]
end -- for jj
nvec[ii]=val
for jj=ii,vlen do
nvec[jj+1]=ovec[jj]
end -- for jj
qflag=1
elseif val==ovec[ii] then
-- val is already in ovec,
-- so don't put a 2nd val into nvec.
for jj=1,vlen do
nvec[jj]=ovec[jj]
end -- for jj
qflag=1
end -- if val
end -- if qflag
end -- for ii
if qflag==0 then
-- val must be > ovec[vlen],
-- so set nvec[vlen+1] to val.
for jj=1,vlen do
nvec[jj]=ovec[jj]
end -- for jj
nvec[vlen+1]=val
end -- if qflag
return nvec
end -- add2vec()
-- Below makes a string mssg
-- containing the vlen elements
-- of the vector vec with single
-- spaces before each element.
--
function listvec(vec)
local mssg,vlen,ii
vlen=(#vec) -- # of elements in vec
mssg=''
for ii=1,vlen do
mssg=(mssg..' '..vec[ii])
end -- for ii
return mssg,vlen
end -- listvec()
-- Below makes a string mssg
-- containing the vlen elements
-- of the vector vec with single
-- spaces before each element.
-- Assuming the elements of vec
-- are segment #'s, it adds AA
-- single-letter abbreviations
-- to mssg for each element.
--
function listvecNaas(vec)
local mssg,vlen,ii,aa
vlen=(#vec) -- # of elements in vec
mssg=''
for ii=1,vlen do
aa=structure.GetAminoAcid(vec[ii])
mssg=(mssg..' '..aa..vec[ii])
end -- for ii
return mssg,vlen
end -- listvecNaas()
-- Above lists all functions.
Main() -- start program in Main() function.
