Recipe: Sepsis-2 Bonding 1.01 - Brow42

created by brow42

Profile

• Name: Sepsis-2 Bonding 1.01 - Brow42
• ID: 45709
• Shared with: Public
• Parent: None
• Children:
  • Vancomycine Bonding 1.02 - Brow42
  • Ligand Sepsis-2 Bonding 1.01 - Brow42 (COPY)
• Created on: June 08, 2013 at 21:41 PM UTC
• Updated on: June 08, 2013 at 21:41 PM UTC
• Description:

  Find bonds or almost-bonds to the ligand. Preserves bands. Makes bands. Only for the 2nd series of sepsis puzzles 680,684....

Code

--[[
  * Ligand H-Bond Utility
  
  * Makes bands where the protein bonds (or almost bonds) the ligand
  * Uses fsl.atom library to identify protein bonding atoms
  * Ligand atoms must be hard-coded into the recipe.
  
  * Version 1.0 Mar. 16 2013 Brow42
  * Programmed with the atoms for 680 Sepsis Solo Hand Folding
  * Version 1.01 Jun 8 2013 Brow42
  * Made range for "almost bonded" slider 3x bonded slider range
--]]

options = {
    do_sidechains = false,
    bond_lower = 1.6, -- based on quick testing with bands and bonding subscore
    bond_upper = 3.0, -- see above
    bond_almost = 3.6, -- Flexible...larger number, more bands and multiple bands per atom
    keep = true, -- delete all bands at end or not
    strength = 0.1,
    reset_bonded = false, -- change lengths of bonded bands
    reset_almost = false -- change lengths of almost bonded bands
}
options.bonded_length = (options.bond_lower + options.bond_upper) / 2.0
options.almost_length = options.bond_upper * 0.9

title = "Sepsis-2 Bonding Utility v. 1.01"

-- data specific to Sepsis series #2 (680, 684...)
valid_puzzle_id = { 994686, 995294 } -- all identical puzzles
valid_puzzle_name = "680 Solo Hand-Folding Sepsis Puzzle"
-- first_aa = { 39 } -- the first segment of each unlocked range or leave blank to do all unlocked
-- dolist = [] -- list of segments to do, or leave blank to autodetect unlocked segments
ligand_acceptors = { 1, 4, 6, 8, 10, 12, 16, 19, 21, 22, 24, 31 }
ligand_polarH = { 27, 29, 31, 35, 41, 44, 48 }


-- ================================== Begin New Dialog Library Functions
-- Add a wall of text from a table
function dialog.AddLabels(d,msg,nlabels) -- pass in # of existing autolabels
    local nlabels = nlabels or #(d._Order or {}) -- default, valid if never delete dialog elements
    if type(msg) == 'string' then
        msg = { msg }
    end
    for i = 1,#msg do
        d['autolabel'..tostring(i+nlabels)] = dialog.AddLabel(msg[i])
    end
end
-- Create but don't display a wall of text and 1 or 2 buttons
function dialog.CreateMessageBox(msg,title,buttontext1,buttontext0)
    title = title or ''
    local d = dialog.CreateDialog(title)
    dialog.AddLabels(d,msg)
    buttontext1 = buttontext1 or 'Ok'
    d.button = dialog.AddButton(buttontext1,1)
    if buttontext0 ~= nil then d.button0 = dialog.AddButton(buttontext0,0) end
    return d
end
-- Display a dialog box
function dialog.ShowMessageBox(msg,title,buttontext1,buttontext0)
    return dialog.Show(dialog.CreateMessageBox(msg,title,buttontext1,buttontext0))
end
-- Build a multipage dialog, returns table of dialogs
-- pages is a table contain strings or tables of strings
-- title is table or string
function dialog.CreateBook(pages,title,buttontext1,buttontext0)
    local dialogs = {}
    for i = 1,#pages do
        local t = type(title) == "string" and title or title[i]
        local d = dialog.CreateMessageBox(pages[i],t,buttontext1,buttontext0)
        if i > 1 then d.prev = dialog.AddButton("Prev",2) end
        if i < #pages then d.next = dialog.AddButton("Next",3) end
        dialogs[#dialogs+1] = d
    end
    return dialogs
end
-- Display a multipage dialog
function dialog.ShowBook(dialogs,page)
    local rc
    page = page or 1
    repeat
        rc = dialog.Show(dialogs[page])
        if rc == 2 then page = page - 1 end
        if rc == 3 then page = page + 1 end
    until rc == 1 or rc == 0
    return rc
end

unpack = unpack or table.unpack -- for 5.2
-- ========================== Begin AA Atom Database ================
fsl = fsl or {} -- make sure fsl namespace exists
fsl.atom = {}
-- Take just these 4 things if all you need is to know if the sidechain starts at 5 or 6 (tail terminal = true)
fsl.atom.atomcount = { -- shorter table just for identifying terminals
    a=10, c=11, d=12, e=15, f=20, g=7, h=17, i=19, k=22, l=19,
    m=17, n=14, p=15, q=17, r=24, s=11, t=14, v=16, w=24, y=21
}

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

function fsl.atom._IsTerminalTest(aa,count,disulfide)
    local dbvalue = fsl.atom.atomcount[aa]
    if dbvalue == nil then error('Bad argument to an fsl.atom function (not an amino acid code)') end
    local diff = count - dbvalue
    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
function fsl.atom.IsTerminal(iSeg)
    if structure.GetSecondaryStructure(iSeg) == 'M' then return false,false,false end
    local aa = structure.GetAminoAcid(iSeg)
    local count = structure.GetAtomCount(iSeg)
    return fsl.atom._IsTerminalTest(aa,count,fsl.atom._IsDisulfideBonded(iSeg))
end

--     ====== now the rest of Atom Tables =====

fsl._nul = {} -- don't need so many CONSTANT empty tables
-- 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['v'] = { {8}, fsl._nul, fsl._nul, 16 } -- valine
fsl.atom.db['w'] = { {15,20}, {9},  fsl._nul, 24 } -- tryptophan
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

function fsl.atom.CountAtoms(iSeg) return structure.GetAtomCount(iSeg) 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
function fsl.atom.GetExpectedCount(iSeg, isFirst, isLast, isDisulfideBonded)
    local db = fsl.atom._NumberOrCode(iSeg)
    if db == nil then return fsl.atom.CountAtoms(iSeg) end
    if isFirst == nil and type(iSeg) == 'number' then isFirst, isLast = fsl.atom.IsTerminal(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
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
    if isFirst == nil and type(iSeg) == 'number' then isFirst, isLast = fsl.atom.IsTerminal(iSeg) 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
    if isFirst == nil and type(iSeg) == 'number' then isFirst, isLast = fsl.atom.IsTerminal(iSeg) 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
    if isFirst == nil and type(iSeg) == 'number' then isFirst, isLast = fsl.atom.IsTerminal(iSeg) 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
    if isFirst == nil and type(iSeg) == 'number' then isFirst, isLast = fsl.atom.IsTerminal(iSeg) 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
    if isFirst == nil and type(iSeg) == 'number' then isFirst, isLast = fsl.atom.IsTerminal(iSeg) 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 ================

-- Returns 1 if the current puzzle is recognized as Sepsis #2
function CheckPuzzleID()
    for i = 1, #valid_puzzle_id do
        if puzzle.GetPuzzleID() == valid_puzzle_id[i] then return 1 end
    end
    local rc = dialog.ShowMessageBox(
        {'This recipe knows atom positions only for',valid_puzzle_name,
        'Are you sure you want to continue?','(this will only work if the puzzle is the same)'}, 'Puzzle Mismatch', 'Yes','Quit')
    return rc
end

-- Menu dialogs, 1 for continue, 0 for quit
function MainMenu()
    local slider_max, slider_min = 5.0,1.5 -- all same for clarity
    local d, rc, h
    d = dialog.CreateDialog(title)
    dialog.AddLabels(d, 'This recipe is for the 2nd series of sepsis puzzles.')
    d['dosc'] = dialog.AddCheckbox('Also check sidechain bonds',options.do_sidechains)
    d['lower'] = dialog.AddSlider('Bond Low',options.bond_lower      ,slider_min, slider_max, 2)
    d['upper'] = dialog.AddSlider('Bond Upper',options.bond_upper    ,slider_min, slider_max, 2)
    d['almost'] = dialog.AddSlider('Almost Upper',options.bond_almost,slider_min, slider_max*3, 2)
    d['keep'] = dialog.AddCheckbox('Keep bands after (dont\' delete)',options.keep)
    d['str'] = dialog.AddSlider('Band Str',options.strength,0.1, 10, 1)
    d['resetbond'] = dialog.AddCheckbox('Change bonded band lengths',reset_bonded)
    d['resetalmost'] = dialog.AddCheckbox('Change almost-bonded band lengths',reset_almost)
    d['bondlen'] = dialog.AddSlider('Bonded Len',options.bonded_length, 0.0, 6.0, 2)
    d['almostlen'] = dialog.AddSlider('Almost Len',options.almost_length, 0.0, 6.0, 2)
    
    d['okay'] = dialog.AddButton('Okay',1)
    d['cancel'] = dialog.AddButton('Cancel',0)
    d['about'] = dialog.AddButton('About',2)
    
    h = dialog.CreateMessageBox({
        'This recipe is built for the 2nd series of sepsis puzzles',
        'starting with 680. Each ligand bonding recipe has to',
        'be hand-coded with the list of ligand bonding atoms.',
        'This recipe can be used on puzzles identical to 680.',
        '',
        'This recipe will find all the donor-acceptor pairs',
        'between the protein and the ligand that fall within a',
        'certain distance. These are assumed to be bonded.',
        'Atoms a little bit farther apart are almost bonded.',
        '',
        'You can keep the bands that are bonded or almost-',
        'bonded. By default, the bands are weak and lock the',
        'atoms in place, but you can change the length to pull',
        'them closer.',
        '',
        'Banding scripts may run faster when minimized. Try',
        'making extra bands and delete ones you don\'t want.'
        }, 'About '..title, 'Back')

    repeat
        rc = dialog.Show(d)
        
        if (rc == 0) then return 0 end
        options.do_sidechains = d['dosc'].value
        options.bond_lower = d['lower'].value
        options.bond_upper = d['upper'].value
        options.bond_almost = d['almost'].value
        options.keep = d['keep'].value
        options.strength = d['str'].value
        options.reset_bonded = d['resetbond'].value
        options.reset_almost = d['resetalmost'].value
        options.bonded_len = d['bondlen'].value
        options.almost_lenght = d['almostlen'].value
        if not (options.bond_lower <= options.bond_upper) or not (options.bond_upper <= options.bond_almost) then
            dialog.ShowMessageBox({'Lower bound must be lowest and',' Almost Bonded bound must be largest'},'Input Error','Oops')
            rc = -1
        end
        if (rc == 2) then dialog.Show(h) end
    until rc == 1
    return rc
end

-- ======= Begin Recipe functions =======

UNBONDED,BONDED,ALMOST = 0,1,2
function BandClassify(iBand)
    if iBand == 0 then return UNBONDED end
    local len = band.GetLength(iBand)
    if len < options.bond_lower or len > options.bond_almost then return UNBONDED end
    if len < options.bond_upper then return BONDED end
    return ALMOST
end

-- Band a bonding pair and keep if bonded or almost bonded
-- Set strength and goal length
function BandOnePair(seg1, seg2, atom1, atom2)
    iBand = band.AddBetweenSegments(seg1,seg2,atom1,atom2)
    if iBand then
        band.SetStrength(iBand,options.strength)
        band.SetGoalLength(iBand,band.GetLength(iBand))
        class = BandClassify(iBand)
        if class == BONDED then
            NBonded = NBonded + 1
            if options.reset_bonded then
                band.SetGoalLength(iBand,options.bonded_length)
            end
        elseif class == ALMOST then 
            NAlmost = NAlmost + 1
            if options.reset_almost then
                band.SetGoalLength(iBand,options.almost_length)
            end
        else
            band.Delete(iBand)
        end
    end
end

-- band all donor-acceptor pairs, keep only bonded or almost bonded
function BandPairs(seg1, seg2, list1, list2)
    for i = 1,#list1 do
        for j = 1,#list2 do
            BandOnePair(seg1,seg2,list1[i],list2[j])
        end
    end
end

function DeleteBands()
    if NBands == 0 then band.DeleteAll() end
    while band.GetCount() > NBands do
        band.Delete(band.GetCount())
    end
end

function ErrorHandler(errmsg)
    if string.find(errmsg,"Cancelled") then
        print("User cancel.") -- print this instead of error message 
        DeleteBands()
    else -- it's a real error, not a cancel
        -- You may also want to do cleaning up stuff here
        print(errmsg) -- one place to print the error
    end
    return errmsg
end

-- Finds unlocked segments with start points given by first_aa
-- OR just finds all unlocked segments
function FindUnlockedSegments()
    local dolist = {}
    if first_aa == nil or #first_aa == 0 then
        for i = 1,structure.GetCount()-1 do
            if not structure.IsLocked(i) then dolist[#dolist + 1] = i end
        end
    else
        for iStart = 1,#first_aa do
            i = first_aa[iStart]
            while not structure.IsLocked(i) do
                dolist[#dolist + 1] = i
                i = i + 1
            end
        end
    end
    return dolist
end

-- ============ Begin Main ==============

print(puzzle.GetName(),os.date())

if CheckPuzzleID() == 0 or MainMenu() == 0 then
    print("Quit")
    return
end

print(string.format('Bond lengths are %4.2f to %4.2f',options.bond_lower,options.bond_upper))
print(string.format('Almost-bonded lengths are %4.2f to %4.2f',options.bond_upper,options.bond_almost))
print('Checking sidechains:',(options.do_sidechains and 'Yes' or 'No'))

iLigand = structure.GetCount()
NBonded,NAlmost = 0,0
NBands = band.GetCount()

dolist = dolist or {}
if #dolist == 0 then dolist = FindUnlockedSegments() end

function MainLoop()
    for i = 1,#dolist do
        iSeg = dolist[i]
     -- iSeg is never first or last peptide so okay to go by default values of residue type
        isFirst, isLast = fsl.atom.IsTerminal(iSeg)
        aa = structure.GetAminoAcid(iSeg)
        p = fsl.atom.GetPolarHydrogens(aa, is_first, is_last)
        a = fsl.atom.GetAcceptorAtoms(aa, is_first, is_last)
        if not options.do_sidechains then -- throw away sidechain atoms for just backbone
            if isLast then a = { a[1],a[2] } -- last residue has 2 O
            else aa = {a[1]}
            end
            if isFirst then p = { p[1],p[2] } -- first residue has 2 H
            else p = {p[1]}
            end
        end
        BandPairs(iSeg,iLigand,p, ligand_acceptors)
        BandPairs(iSeg,iLigand,a, ligand_polarH)
    end
end

rc = xpcall(MainLoop,ErrorHandler)

if rc == true then
    if options.keep == false then DeleteBands() end

    print("# of bonds to ligand:",NBonded)
    print("# almost bonded:",NAlmost)
end

Comments

[brow42]

I made the slider for almost bands to go up to 15, so you can keep LOTS of bands if you want.

Also, since Sepsis 3 uses the same ligand and Sepsis-2, I've taken out the warning for puzzle 727.