Recipe: Atom Number Identifier v1.0.1 -- Brow42

created by brow42

Profile

• Name: Atom Number Identifier v1.0.1 -- Brow42
• ID: 39989
• Shared with: Public
• Parent: None
• Children:
  • Atom Number Identifier v1.0.1 -- Brow42versuch
• Created on: March 24, 2012 at 05:45 AM UTC
• Updated on: March 24, 2012 at 05:45 AM UTC
• Description:

  1) go to view-atom mode 2) make a small number of bands on a small number of segments (< 3). Select or freeze the segments. 4) Run script to get atom numbers.

Code

--[[
 * Atom Number Identifier
 
 * Original Author: Brow42 Mar. 23 2012
 
 * Version 1.0
 
 * This script will loop through all the atoms of a segment, looking
 * for ones that have bands attached. It doesn't know which segment
 * to check, so you have to mark them using selection interface or
 * freezing them (backbone or sidechain). The more you select, the
 * more atoms will be tested for each band you have made, so test
 * 3 segments at a time.
 *
 * To identify the number of an atom, make a band from it, freeze
 * or select the segment, and run the script. You can make multiple
 * bands on multiple segments; output will be in the order the
 * bands were created.
 *
 * Program will output the segment and atom numbers for the STARTING
 * point of each band you have made (doesn't matter where the other
 * end is). Note that the first and last segments of a polypeptide
 * have different atom numbering due to extra atoms.

 * 1.0.1 bug fix on found array (had holes, wasn't an array) 2012/3/23
--]]

-- Distance from an atom to the end of a band, measured by making a test band
-- call in protected mode unless you KNOW the atom number is valid!!!!
function Distance(iSeg, iAtom, iBand)
    local jBand = band.AddToBandEndpoint(iSeg, iBand, iAtom)
    if jBand < 1 then error(string.format('Unspecified error banding from seg %d, atom %d to band %d',iSeg,iAtom,iBand)) end
    local d = band.GetLength(jBand)
    band.Delete(jBand)
    return d
end

_lastatomerr = 'atom number out of bounds'
function LastAtom(errmsg)
    return string.find(errmsg,_lastatomerr) ~= nil
end

function IsFrozen(iSeg) local a,b = freeze.IsFrozen(iSeg) return a or b end

--
-- ---------------------------- Begin Main
--

if band.GetCount() == 0 then
    print('No bands! First make bands from atoms, then mark the segments.')
    return
end

tocheck = {} -- segments to check
tofind = {} -- bands to find
found = {} -- bands we've found

-- initial list of bands
for i = 1,band.GetCount() do tofind[i] = i end

-- Find the segments to check
print('Looking for selected or frozen segments which\nwill be checked for bands')
for i=1,structure.GetCount() do
    if IsFrozen(i) or selection.IsSelected(i) then
        tocheck[#tocheck+1] = i
    end
end

-- Tell the user what segments we'll be checking
if #tocheck == 0 then
    print('No segments marked for checking. You must\nmark only a few segments, otherwise')
    print('it will take a really long time to run.')
else
    print('Found',#tocheck,'segments')
    if #tocheck > 3 then print('That\'s a lot, recommend no more than 3 segments\nunless they are really small.') end
    print('Segs:',table.concat(tocheck,', '))
end

-- while looop of list of bands in reverse order
while #tofind > 0 do
    local iAtom, iBand, iSeg
    iBand = tofind[#tofind]
    for iCheck = 1,#tocheck do
        iSeg = tocheck[iCheck]
        dist  = band.GetLength(iBand)
        iAtom = 0
        local rc, d -- test band length OR error message
        -- loop over atoms here. we don't know how many atoms there are, so catch the error
        repeat
            iAtom = iAtom + 1
            rc, d = pcall(Distance,iSeg, iAtom, iBand)
            if rc == false then break end -- out of atoms, OR some other error
        until d == dist -- found the correct atom
        if rc == true then -- found the atom
            found[#found+1] = { iBand, iSeg, iAtom }
            break -- break out of segment loop
        elseif not LastAtom(d) then print(d) -- was error message an unknown error? Print it.
        end
    end
    tofind[#tofind] = nil -- We found the atom or we didn't, either way, done looking for this band
end

print(#found,'atoms identified.')
if #found == 0 then print('Perhaps you forgot to mark the segments?') end

for i = #found,1,-1 do
    print(string.format('Band %d = seg %d, atom %d', unpack(found[i])))
end
print('Done.')

Comments

[brow42]

This script will loop through all the atoms of a segment, looking for ones that have bands attached. It doesn't know which segment to check, so you have to mark them using selection interface or freezing them (backbone or sidechain). The more you select, the more atoms will be tested for each band you have made, so test 3 segments at a time.

To identify the number of an atom, make a band from it, freeze or select the segment, and run the script. You can make multiple bands on multiple segments; output will be in the order the bands were created.

Program will output the segment and atom numbers for the STARTING point of each band you have made (doesn't matter where the other end is). Note that the first and last segments of a polypeptide have different atom numbering due to extra atoms.

This script is just a temporary measure. Eventually all the information should be in a table, like the old hydrophobe table. The script will still be useful for non-amino-acid ligands.

To see atoms, go to (advanced) View, pick Stick + H or Polar H. Set the coloring to CKY. Blue = Nitrogen (#1 on the backbone), grey/other = C, Red = Oxygen, White = Hydrogen. Polar H only shows hydrogens attached to a hydrogen bond donor.

Keyboard shortcut: Ctrl-Shift-V