Recipe: H Bonds 532

created by brow42

Profile

• Name: H Bonds 532
• ID: 40035
• Shared with: Public
• Parent: None
• Children:
  • H Bonds 539
  • H Bonds 539
• Created on: March 25, 2012 at 09:25 AM UTC
• Updated on: March 25, 2012 at 09:25 AM UTC
• Description:

  Bonds donor/acceptor atoms nearest to the ligand. Experimental.

Code

--[[
 * H bond 532
 
--]]
-- Options
selected_segments = {229, 230} -- ligand segments
radius = 10.0  -- search radius
trim = true -- delete all but one band for each ligand atom
strength = 1.0 -- band strength
bondlength = 3-- nominal bond length (guess)

fsl = fsl or {} -- make sure fsl namespace exists
unpack = unpack or table.unpack -- for 5.2
-- ========================== Begin AA Atom Database ================
fsl._nul = {} -- don't need so many CONSTANT empty tables
fsl.atom = {}
-- 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['w'] = { {15,20}, {9},  fsl._nul, 24 } -- tryptophan
fsl.atom.db['v'] = { {8}, fsl._nul, fsl._nul, 16 } -- valine
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

-- Pass in a segment number of a cysteine
function fsl.atom._IsDisulfideBonded(iSeg)
    local s = current.GetSegmentEnergySubscore(iSeg,'disulfides')
    return tostring(s) ~= '-0'
end

fsl.atom._lastatomerr = 'atom number out of bounds'
-- Find out how many atoms in a segment. SLOW!
-- We do this by trying to band atoms until we get an error
function fsl.atom.CountAtoms(iSeg)
    local function LastAtom(errmsg)
        return string.find(errmsg,fsl.atom._lastatomerr) ~= nil
    end
    local s2,s3,n
    n = structure.GetCount()
    if iSeg == 1 then s2,s3 = 2,3
    elseif iSeg == n then s2,s3 = n-1, n-2
    else s2,s3 = iSeg-1,iSeg+1
    end
    local function zlb(iSeg, iAtom)
        local rc,iBand = pcall(band.Add,iSeg,s2,s3, 0.01,0,0, iAtom)
        if rc == false then error(iBand) end
        if iBand == 0 then error(string.format('Unknown Band Add Error on seg %d atom %d',iSeg,iAtom)) end
        band.Delete(iBand)
    end
    local iAtom = 0
    while true do
        iAtom = iAtom + 1
        rc, err = pcall(zlb,iSeg, iAtom)
        if rc == false then break end -- out of atoms, OR some other error
    end -- found the correct atom
    if not LastAtom(err) then error(err) end
    return iAtom - 1
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

-- iSeg is a segment number or amino acid letter code
-- if iSeg is a ligand it will manually count the atoms
function fsl.atom.GetCount(iSeg, isFirst, isLast, isDisulfideBonded)
    local db = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return fsl.atom.CountAtoms(iSeg) end
    local c = db.count
    if isFirst then c = c + 2 end
    if isLast then c = c + 1 end
    if isDisulfideBonded then c = c - 1 end -- S-H + H-S -> S-S + 2H
    return c, aa
end

-- if you already have the count saved, call this instead of IsTerminal()
function fsl.atom._IsTerminalTest(aa,count,disulfide)
    local db = fsl.atom.db[aa]
    if db == nil then error('Bad argument to an fsl.atom function (not an amino acid code)') end
    local diff = count - db.count
    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
-- SLOW
function fsl.atom.IsTerminal(iSeg)
    if structure.GetSecondaryStructure(iSeg) == 'M' then return false,false end
    local aa = structure.GetAminoAcid(iSeg)
    local count = fsl.atom.CountAtoms(iSeg)
    return fsl.atom._IsTerminalTest(aa,count,fsl.atom._IsDisulfideBonded(iSeg))
end

-- Get a list of backbone heavy atoms always (1,4) or (1,5)
function fsl.atom.GetBackboneHeavyAtoms(iSeg, isFirst, isLast)
    local db = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return end
    return isLast and {1, 2, 3, 4, 5} or { 1, 2, 3, 4}
end

-- Sidechain heavy atoms, starting with C-beta
-- Returns nil for ligands
function fsl.atom.GetSidechainHeavyAtoms(iSeg, isFirst, isLast)
    local db, aa = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return end
    local a,b
    if aa == 'g' then return {} -- g has no sidechain
    elseif aa == 'p' then a,b = 5,7 -- p no polar hydrogen to tell us this range
    else a,b = 5, db.polarh[1]-1 -- everything else has sidechain followed by a polar hydrogen in table
    end
    if isLast then a,b = a+1,b+1 end -- shift due to terminal O at 5
    local list = {}
    for i = a,b do list[#list + 1] = i end
    return list
end

-- Get list of donor atoms
-- Returns nil for ligands
function fsl.atom.GetDonorAtoms(iSeg, isFirst, isLast)
    local db,aa = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return nil end
    local list
    if aa == 'p' then list = {} -- {1} -- is this true?!
    else list = {1} -- the backbone N
    end
    local shift = 0
    if isLast then shift = 1 end -- shift due to terminal O, all remaining donors are on sidechain
    for i = 1,#db.donor do list[#list+1] = db.donor[i] + shift end
    return list
end

-- Get list of acceptor atoms
-- Returns nil for ligands
function fsl.atom.GetAcceptorAtoms(iSeg, isFirst, isLast)
    local db = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return nil end
    local list = { 4 } -- The O attached to the backbone C'
    local shift = 0
    if isLast then
        shift = 1
        list = {4 , 5} -- 2 O on the terminus
    end
    -- all remaining acceptors are on the sidechain
    for i = 1,#db.acceptor do list[#list+1] = db.acceptor[i] + shift end
    return list
end

-- Get list of polar hydrogens atoms
-- Returns nil for ligands
function fsl.atom.GetPolarHydrogens(iSeg, isFirst, isLast)
    local db,aa = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return nil end
    local list = {}
    local tmp, first -- our list of H, and where the first H will go
    if aa == 'p' then tmp, first = {}, 8 -- no table, so hardcoded
    else tmp, first = db.polarh, db.polarh[1] -- table provided, first element is first element
    end
    local shift = 0
    if isLast then shift = 1 end -- all H are after the sidechain
    first = first + shift -- Shifted by terminal O
    if isFirst then
        shift = shift + 2 -- Shift the table up after the 2 added H
        list = { first, first + 1 } -- Add 2 terminal H
    end
    for i = 1,#tmp do list[#list+1] = tmp[i] + shift end
    return list
end


-- ========================== End AA Atom Database ================
function PrintTable(tab)
    for i,v in pairs(tab) do print(i,v) end
end

function TerminalCheck(iSeg)
    return iSeg == 1, iSeg == structure.GetCount(), fsl.atom._IsDisulfideBonded(iSeg)
end

function GetAtomLists(iSeg)
    atoms = { donors = {}, acceptors = {}, polar = {} }
    local tmp = fsl.atom.GetDonorAtoms(iSeg, TerminalCheck(iSeg))
    for i = 1,#tmp do atoms.donors[#atoms.donors + 1] = tmp[i] end
    tmp = fsl.atom.GetAcceptorAtoms(iSeg, TerminalCheck(iSeg))
    for i = 1,#tmp do atoms.acceptors[#atoms.acceptors + 1] = tmp[i] end
    tmp = fsl.atom.GetPolarHydrogens(iSeg, TerminalCheck(iSeg))
    for i = 1,#tmp do atoms.polar[#atoms.polar + 1] = tmp[i] end
    return atoms
end

-- Select all segmements withing radius distance of selected segments
-- Exclude neighbors of selected segments
function SphereSelect(iSeg,radius)
    local list = {}
    local n = structure.GetCount()
    for i = 1,n do
        if structure.GetDistance(i,iSeg) < radius then
            list[i] = true
        end
    end

    list[iSeg] = nil
    list[iSeg+1] = nil
    list[iSeg-1] = nil

    local list2 = {}
    for i,v in pairs(list) do list2[#list2 + 1] = i  end
    return list2
end

function BandPairs(origin_segment,origin,targets,group1,group2,trim)
    local n_existing = band.GetCount()
    local dists
    for i = 1,#origin[group1] do
        dists = {}
        bands = {}
        for j,list in pairs(targets) do
            for k = 1,#list[group2] do
                local iBand = band.AddBetweenSegments(origin_segment,j,origin[group1][i],list[group2][k])
                if iBand > 0 then
                    d = band.GetLength(iBand)
                    dists[#dists+1] = d
                    bands[#bands+1] = iBand
                end
            end -- all of the atoms in segment j
        end -- all of the target segments
        if trim then
            local mindist = math.min(unpack(dists))
            for m = #bands,1,-1 do
                if dists[m] ~= mindist then band.Delete(bands[m]) end
            end
        end
    end -- all origin atoms
    -- set params of non-trimmed bands
    for i = n_existing + 1, band.GetCount() do
        band.SetStrength(i,strength)
        band.SetGoalLength(i,bondlength)
    end
end

for i = 1, #selected_segments do
    local iSeg = selected_segments[i]
    ligand_atoms = GetAtomLists(iSeg)
    neighbors = SphereSelect(selected_segments[i], radius) 
    table.sort(neighbors)
    print('Selected neighbors for',iSeg,table.concat(neighbors,', '))
    local neighbor_atoms = {}
    for j = 1,#neighbors do
        neighbor_atoms[neighbors[j]] = GetAtomLists(neighbors[j])
    end
    BandPairs(iSeg, ligand_atoms,neighbor_atoms,'donors','acceptors',trim)
    BandPairs(iSeg, ligand_atoms,neighbor_atoms,'acceptors','donors',trim)
end

Comments

[pedroff_1]

what am I supposed to do?