Recipe: UnBUN_1.1

created by mwm64

Profile

• Name: UnBUN_1.1
• ID: 103691
• Shared with: Public
• Parent: None
• Children: None
• Created on: July 09, 2020 at 14:09 PM UTC
• Updated on: July 09, 2020 at 14:09 PM UTC
• Description:

  Script that tries to make bonds to the target protein and reduce BUNS.

Code

-- unbun
-- Try to get rid of  BUNS by finding hydrophilic segments in the
-- target protein  that are near segments in the user protein
-- and mutating the user protein segment  to something that can 
-- bond with that target; keeping changes that improve our BUN
-- score

iterations=2    -- iterations for shake and wiggle
bond_radius=9.6   -- how close segments have to be to try to make bonds

function GetBondableAcids(acid)
    -- if this is a "red" acid, return the list of blue ones
    -- and vice versa; if its a red-blue one return both
    -- and if its a M, return the same...
    if acid == 'k' or acid == 'r' or acid == 'h' or acid == 'w' then
        return {'e','d','s','t','y'}
    end
    if acid == 'e' or acid== 'd' or acid =='s' or acid == 't' or acid == 'y' then
        return {'k','r','h','w'}
    end
    if acid == 'm' then
        return {'m'}
    end
    -- otherwise pretend its a red-blue, and try lots..
    return {'k','r','h','w','e','d','s','t','y'}
end

function find_filter()
    -- first find the BUN filter so we can check for improvement
    --
    satfilt = -1
    filt = filter.GetNames ()
    print ( #filt .. " conditions or filters" )
    for ii = 1, #filt do
        print( "checking filter " .. ii .. ": " .. filt[ii] )
        if string.match(filt[ii],".*sat") then
            print( "found sat filter " )
            satfilt = ii
        end
    end
end

function delay2() 
   local xyzzy
   for i=1,1000 do
      xyzzy=5
   end
end

function delay()
   for i=1,10  do
       delay2()
   end
end


function do_unbun()
    if satfilt >= 0 then
        nsegs = structure.GetCount()
        -- guess that the boundary between target and user segments is halfway
        -- we'll move it as we see mutable / non mutables...
        boundary = math.floor( nsegs/2 )
        orig_nsat_bonus = filter.GetBonus(filt[satfilt])
        print("estimating effort:")
        wiggles_estimated = math.floor(nsegs/10)
        wiggles_done = 0
        saw_mutable = false
        move_boundary = 0
        for tps = 1, boundary do
            pest = math.floor(tps * 100 / boundary)
            if tps % 15 == 1 then
                print("estimating: ["..pest.."%]")
            end
            if not structure.IsHydrophobic(tps) then
                taa = structure.GetAminoAcid(tps)
                bonding_acids = GetBondableAcids(taa)
                for ups = boundary, nsegs do
                    if not saw_mutable and not structure.IsMutable(ups) then
                        move_boundary = move_boundary + 1
                    end
                    if not saw_mutable and structure.IsMutable(ups) then
                        saw_mutable = true
                    end 
                    if structure.IsMutable(ups) and structure.GetDistance(tps, ups) < bond_radius and string.len(taa) == 1 then
                       wiggles_estimated = wiggles_estimated + #bonding_acids
                    end
                end
            end
        end

        if boundary + move_boundary < nsegs then
            boundary = boundary + move_boundary
        end
        print("Trying to find bonds for Hydrophobics in  first " .. boundary .. " segments...")
        print("Estimating " .. wiggles_estimated .. " wiggles of " .. iterations .. " iterations")

        for tps = 1, boundary do
            if not structure.IsHydrophobic(tps) then
                -- look up amino acids to try against this protein
                taa = structure.GetAminoAcid(tps)
                bonding_acids = GetBondableAcids(taa)
                pdone = math.floor(100 * wiggles_done / wiggles_estimated)
            
                for ups = boundary, nsegs do
                    maxbonusacid = structure.GetAminoAcid(ups)
                    if structure.IsMutable(ups) and structure.GetDistance(tps, ups) < bond_radius and string.len(maxbonusacid) == 1 then
                        print("["..pdone.."%] Segment " .. ups .. " acid " .. maxbonusacid .. " is within " .. bond_radius ..  " of polar segment " .. tps)
                        segments_tried = segments_tried + 1
                        -- keep max bonus and corresponding protein
                        -- max (so far) is where we are now
                        max_nsat_bonus = filter.GetBonus(filt[satfilt])
                        recentbest.Save()
                        save.Quicksave(99)
                        selection.DeselectAll()
                        selection.Select(tps)
                        selection.Select(ups)
                        selection.Select(ups+1)
                        selection.Select(ups+2)
                       
                        print( "["..pdone.."%] Trying each of: ".. table.concat(bonding_acids,',') )
                        for bai = 1,#bonding_acids do
                            wiggles_done = wiggles_done + 1
                            pdone = math.floor(100 * wiggles_done / wiggles_estimated)
                            structure.SetAminoAcid(ups, bonding_acids[bai])
                            -- sometimes just mutating it adds a bond, so skip shake and or wiggle if
                            -- the score is already better...
                            bonus = filter.GetBonus(filt[satfilt])
                            if false and  bonus <= max_nsat_bonus then
                                -- maybe if we shake it it will bond?
                                structure.ShakeSidechainsSelected(iterations)
                                bonus = filter.GetBonus(filt[satfilt])
                                
                            end
                            if bonus <= max_nsat_bonus then
                                --  see if wiggling with the band will bond it...
                                -- try to band the sidechains so wiggle will try to make the bond
                                tcount = structure.GetAtomCount(tps)
                                ucount = structure.GetAtomCount(ups)

                                myband = band.AddBetweenSegments(tps, ups, tcount, ucount)
                                band.SetGoalLength(myband,1.6)
                                band.SetStrength(myband, 1.3)
                                -- wiggle just sidechains first, so we get as
                                -- close as possible without bending the 
                                -- backbone
                                
                                structure.LocalWiggleSelected(iterations,false,true)
                                bonus = filter.GetBonus(filt[satfilt])
                                if bonus <= max_nsat_bonus then
                                   structure.LocalWiggleSelected(iterations,true,true)
                                   bonus = filter.GetBonus(filt[satfilt])
                                end
                                if bonus <= max_nsat_bonus then
                                   structure.WiggleAll(iterations)
                                end
                                -- Foldit keeps hanging or missing band deletes when we do it right away, so
                                delay()
                                -- THEN delete all -- that way if it misses
                                -- one we get it next time.
                                band.DeleteAll()
                                myband = -1

                                bonus = filter.GetBonus(filt[satfilt])
                            end

                            if bonus > max_nsat_bonus then
                                -- we have a new winner!!!
                                maxbonusacid = bonding_acids[bai]
                                max_nsat_bonus = bonus
                                recentbest.Save()
                                save.Quicksave(99)
                                print("new best!  BUNS score: ".. bonus)
                                break
                            else
                                -- our bonus is equal, but... 
                                -- if we improved  the overall score , keep the wiggle anyway
                                if current.GetEnergyScore() > recentbest.GetEnergyScore() then
                                    recentbest.Save()
                                    save.Quicksave(99)
                                else
                                    save.Quickload(99)
                                end
                            end
                        end
                    end
                end
            end
        end
        if segments_tried == 0 then
            print("No segments found within bond_raius of " .. bond_radius .. ".")
            print("edit recipe and increase bond_radius.")
        end

    else
        print( "Couldn't find sat filter, quitting" )
    end
end

function putback()
    -- if anything happens, or we cancel, fall back to our recentbest
    -- and get rid of the band we may have left.
    
    print("Cleaning up: putback()")
    -- Foldit keeps hanging when we do stuff too soon after
    -- Wiggle-ing so
    delay()
    -- THEN
    recentbest.Restore()
    if myband > -1 then
        band.Delete(myband)
    end
    undo.SetUndo(true)
end

function unbun_main()
    print("unbun: starting with iterations=" .. iterations .. " and bond_radius=" .. bond_radius )
    myband = -1
    segments_tried = 0
    undo.SetUndo(false)
    recipe.SectionStart("unbun")
    recentbest.Save()
    save.Quicksave(99)
    find_filter()
    -- should be done as below, but this crashes for me
    -- xpcall(do_unbun, putback)
    -- instead:
    do_unbun()
    undo.SetUndo(true)
    recipe.ReportStatus()
    recipe.SectionEnd()
end

unbun_main()

Comments

[Bruno Kestemont]

puzzle 1862:

error line 117: index out of bound