Recipe: TvdL MakeContactbands 1.3.0

created by Timo van der Laan

Profile

• Name: TvdL MakeContactbands 1.3.0
• ID: 49497
• Shared with: Public
• Parent: None
• Children: None
• Created on: June 04, 2016 at 13:12 PM UTC
• Updated on: June 04, 2016 at 13:12 PM UTC
• Description:

  Makes bands for contacts. BandLengthFactors are settable. Defaults are the ones Karen said to use in global. The bands are, I hope, between the contact points. ; New: Added a parameter to prevent pull from to far away.

Code

-- MakeContactBands
-- Makes the ContactBands for a selection or the whole protein
-- Strengths are settable
-- Tvdl, 28-06-2014 Free to use for noncommercial purposes

-- Handy shorts module
normal= (current.GetExplorationMultiplier() == 0)
segCnt=structure.GetCount()
segCnt2=segCnt
while structure.GetSecondaryStructure(segCnt2)=="M" do segCnt2=segCnt2-1 end

-- On request of gmn
CIfactor=1
maxCI=true
function CI(CInr)
    if CInr > 0.99 then maxCI=true else maxCI=false end
    behavior.SetClashImportance(CInr*CIfactor)
end

function CheckCI()
   local ask=dialog.CreateDialog("Clash importance is not 1")
   ask.l1=dialog.AddLabel("Last change to change it")
   ask.l2=dialog.AddLabel("CI settings will be multiplied by set CI")
   ask.continue=dialog.AddButton("Continue",1)
   dialog.Show(ask)
end

if behavior.GetClashImportance() < 0.99 then CheckCI() end
CIfactor=behavior.GetClashImportance()

-- Score functions
function Score(pose)
    if pose==nil then pose=current end
    local total= pose.GetEnergyScore()
    -- FIX for big negatives

    if normal then
        return total
    else
        return total*pose.GetExplorationMultiplier()
    end
end

function SegScore(pose)
    if pose==nil then pose=current end
    local total=8000
    for i=1,segCnt2 do
        total=total+pose.GetSegmentEnergyScore(i)
    end
    return total
end

function RBScore()
    return Score(recentbest)
end

function round3(x)--cut all afer 3-rd place
    return x-x%0.001
end

-- Module Filteractive

Filterscore=Score()
behavior.SetSlowFiltersDisabled(true)
FilterOffscore=Score()
behavior.SetSlowFiltersDisabled(false)
maxbonus=Filterscore-FilterOffscore

function Filter()
   local ask=dialog.CreateDialog("Slow filters seem to be active")
   ask.disable=dialog.AddCheckbox("Run with disabled slow filters",Filteractive)
   ask.l1=dialog.AddLabel("Current bonus is: "..maxbonus)
   ask.l2=dialog.AddLabel("If this is not the maximum bonus put in a number")
   if maxbonus < 0 then maxbonus=0 end
   ask.maxbonus=dialog.AddTextbox("Set maxbonus:",maxbonus)
   ask.l3=dialog.AddLabel("Scores will only be checked for real gains if")
   ask.l4=dialog.AddLabel("Score with filter off+maxbonus is a potential gain")
   ask.continue=dialog.AddButton("Continue",1)
   dialog.Show(ask)
   maxbonus=ask.maxbonus.value
if maxbonus=="" then maxbonus=0 end
   Filteractive=ask.disable.value
end

BetterRecentBest=false
function FilterOff()
    -- Filters off but restore a better recentbest with filter off
    behavior.SetSlowFiltersDisabled(true)
    if BetterRecentBest then
        save.Quicksave(99)
        save.Quickload(98)
        recentbest.Save()
        save.Quickload(99)
    end
end

function FilterOn()
    -- Filter on but remember recent best if better than current
    BetterRecentBest= Score(recentbest) > Score()
    if BetterRecentBest then
        save.Quicksave(99)
        recentbest.Restore()
        save.Quicksave(98)
        save.Quickload(99)
    end
    behavior.SetSlowFiltersDisabled(false)
end

Filteractive=(math.abs(maxbonus) > 0.1)
if Filteractive then
   --Filters active, give people a choice
   --And ask what the maximum bonus is.
   Filter()
end
-- End of module Filteractive

bestScore=Score()
if Filteractive then FilterOff() end
function SaveBest()
  if (not Filteractive) or
     (Score()+maxbonus>bestScore) then
     if Filteractive then FilterOn() end
     local g=Score()-bestScore
     if g>0 then
        if g>0.01 then print("Gained another "..round3(g).." pts.") end
        bestScore=Score()
        save.Quicksave(3)
     end
     if Filteractive then FilterOff() end
  end
end

-- New WiggleFactor
WF=1
-- Wiggle function
-- Optimized due to Susumes ideas
-- Note the extra parameter to be used if only selected parts must be done
function Wiggle(how, iters, minppi,onlyselected)
    --score conditioned recursive wiggle/shake
    --fixed a bug, absolute difference is the threshold now
    if how==nil then how="wa" end
    if iters==nil then iters=3 end
    if minppi==nil then minppi=0.1 end
    if onlyselected==nil then onlyselected=false end
    local wf=1
    if maxCI then wf=WF end

    --if iters>0 then
        --iters=iters-1
        local sp=Score()
        if onlyselected then
            if how == "s" then
                -- Shake is not considered to do much in second or more rounds
                structure.ShakeSidechainsSelected(1)
                return
            elseif how == "wb" then structure.WiggleSelected(2*wf*iters,true,false)
            elseif how == "ws" then structure.WiggleSelected(2*wf*iters,false,true)
            elseif how == "wa" then structure.WiggleSelected(2*wf*iters,true,true)
            end
        else
            if how == "s" then
                -- Shake is not considered to do much in second or more rounds
                structure.ShakeSidechainsAll(1)
                return
            elseif how == "wb" then structure.WiggleAll(2*wf*iters,true,false)
            elseif how == "ws" then structure.WiggleAll(2*wf*iters,false,true)
            elseif how == "wa" then structure.WiggleAll(2*wf*iters,true,true)
            end
        end
        --if math.abs(Score()-sp) > minppi then return Wiggle(how, iters, minppi,onlyselected) end
    --end
end

-- end of handy shorts module

-- Module for contact puzzles. Free to used for non commercial purposes.

-- Administration of made contactbands
contactbandlist={} -- list of all contactbands
hascontactband={} -- holds bandnumber on position seg1*segCnt2+seg2
KeepBandstrength=0.5
PullBandstrength=1.5
KeepFactor=0.9
PullFactor=0.7
MaxLengthContact=0
MinLengthNoContact=100
MinContactGap=2
MaxNoContactGap=14
PrintNonContacts=false

Nrofstartbands=band.GetCount() --Remember the number of existing bands

function AddContactBand(from,to,keepstr,pullstr)
   local atomfrom=5
   local atomto=5
   if keepstr == nil then keepstr=KeepBandstrength end
   if pullstr == nil then pullstr=PullBandstrength end
   if from>to then from,to = to,from end
   if to-from<MinContactGap then return end -- no band if to close
   if structure.GetAminoAcid(from) == 'g' then atomfrom=2 end
   if structure.GetAminoAcid(to) == 'g' then atomto=2 end
   if contactmap.GetHeat(from,to) > 0 then
      local bnr=hascontactband[from*segCnt2+to]
      if bnr == nil then --new band
         bnr=band.AddBetweenSegments(from,to,atomfrom,atomto)
         hascontactband[from*segCnt2+to]=bnr
         contactbandlist[#contactbandlist+1]=from*segCnt2+to
      end
      local l=band.GetLength(bnr)
      if contactmap.IsContact(from,to) then
         if l>MaxLengthContact then MaxLengthContact=l end
-- print("Contact: ",from,to,"Length: ",l)
         band.SetGoalLength(bnr,l*KeepFactor)
         band.SetStrength(bnr,keepstr)
      else
         if l>MaxNoContactGap then band.Delete(bnr)
         else
            if l<MinLengthNoContact then MinLengthNoContact=l end
            if PrintNonContacts then print(from.."-"..to.." distance: "..l) end
            band.SetGoalLength(bnr,l*PullFactor)
            band.SetStrength(bnr,pullstr)
         end
      end
   end
end

function AddContactsInRange(n1,n2,keepstr,pullstr)
   -- Makes all contactbands inside a range
   for i=n1,n2-2 do
      for j=i+2,n2 do
         AddContactBand(i,j,keepstr,pullstr)
      end
   end
end

function AddContactsRangeToRange(n1,n2,n3,n4,keepstr,pullstr)
   -- Makes all contactbands from one range to another
   for i=n1,n2 do
      for j=n3,n4 do
         AddContactBand(i,j,keepstr,pullstr)
      end
   end
end

function AddContactsSelection(Sel,keepstr,pullstr)
   if keepstr == nil then keepstr=KeepBandstrength end
   if pullstr == nil then pullstr=PullBandstrength end
   for i=1,#Sel-1 do
      for j=i+1,#Sel do
         AddContactsRangeToRange(Sel[i][1],Sel[i][2],Sel[j][1],Sel[j][2],keepstr,pullstr)
      end
      AddContactsInRange(Sel[i][1],Sel[i][2],keepstr,pullstr)
   end
   AddContactsInRange(Sel[#Sel][1],Sel[#Sel][2],keepstr,pullstr)
end

function AddAllContactsRange(s1,s2,keepstr,pullstr)
-- Makes all contactbands to (and inside) a range
    if s1>s2 then s1,s2=s2,s1 end
    for i=s1,s2 do
       for j=1,segCnt2 do AddContactBand(i,j,keepstr,pullstr) end
    end
end

function AddAllContactsSelection(Sel,keepstr,pullstr)
-- Makes all contactbands to (and inside) a selection
   for i=1,#Sel do AddAllContactsRange(Sel[i][1],Sel[i][2],keepstr,pullstr) end
end

function ClrContactBands()
   for i=#contactbandlist,1,-1 do
      band.Delete(hascontactband[i])
      hascontactband[i]=nil
   end
   contactbandlist={}
end

Contactscore=Score()
function UpdateContactBands(keepstr,pullstr)
   if keepstr == nil then keepstr=KeepBandstrength end
   if pullstr == nil then pullstr=PullBandstrength end
   if Contactscore~=Score() then
      print("Updating Bands")
      for i=1,#contactbandlist do
         local n=contactbandlist[i]
         local seg2=n%segCnt2
         local seg1=(n-seg2)/segCnt2
         if seg2==0 then seg2=segCnt2 seg1=seg1-1 end
         AddContactBand(seg1,seg2,keepstr,pullstr)
      end
      Contactscore=Score()
   end
end

-- Segment set and list module
-- Notice that most functions assume that the sets are well formed
-- (=ordered and no overlaps)
-- 02-05-2012 TvdL Free to use for non commercial purposes

function SegmentListToSet(list)
    local result={}
    local f=0
    local l=-1
    table.sort(list)
    for i=1,#list do
        if list[i] ~= l+1 and list[i] ~= l then
            -- note: duplicates are removed
            if l>0 then result[#result+1]={f,l} end
            f=list[i]
        end
        l=list[i]
    end
    if l>0 then result[#result+1]={f,l} end
    --print("list to set")
    --SegmentPrintSet(result)
    return result
end

function SegmentSetToList(set)
    local result={}
    for i=1,#set do
        --print(set[i][1],set[i][2])
        for k=set[i][1],set[i][2] do
            result[#result+1]=k
        end
    end
    return result
end

function SegmentCleanSet(set)
-- Makes it well formed
    return SegmentListToSet(SegmentSetToList(set))
end

function SegmentInvertSet(set,maxseg)
-- Gives back all segments not in the set
-- maxseg is added for ligand
    local result={}
    if maxseg==nil then maxseg=structure.GetCount() end
    if #set==0 then return {{1,maxseg}} end
    if set[1][1] ~= 1 then result[1]={1,set[1][1]-1} end
    for i=2,#set do
        result[#result+1]={set[i-1][2]+1,set[i][1]-1}
    end
    if set[#set][2] ~= maxseg then result[#result+1]={set[#set][2]+1,maxseg} end
    return result
end

function SegmentInvertList(list)
    table.sort(list)
    local result={}
    for i=1,#list-1 do
       for j=list[i]+1,list[i+1]-1 do result[#result+1]=j end
    end
    for j=list[#list]+1,segCnt2 do result[#result+1]=j end
    return result
end

function SegmentInList(s,list)
    table.sort(list)
    for i=1,#list do
        if list[i]==s then return true
        elseif list[i]>s then return false
        end
    end
    return false
end

function SegmentInSet(set,s)
    for i=1,#set do
        if s>=set[i][1] and s<=set[i][2] then return true
        elseif s<set[i][1] then return false
        end
    end
    return false
end

function SegmentJoinList(list1,list2)
    local result=list1
    if result == nil then return list2 end
    for i=1,#list2 do result[#result+1]=list2[i] end
    table.sort(result)
    return result
end

function SegmentJoinSet(set1,set2)
    return SegmentListToSet(SegmentJoinList(SegmentSetToList(set1),SegmentSetToList(set2)))
end

function SegmentCommList(list1,list2)
    local result={}
    table.sort(list1)
    table.sort(list2)
    if #list2==0 then return result end
    local j=1
    for i=1,#list1 do
        while list2[j]<list1[i] do
            j=j+1
            if j>#list2 then return result end
        end
        if list1[i]==list2[j] then result[#result+1]=list1[i] end
    end
    return result
end

function SegmentCommSet(set1,set2)
    return SegmentListToSet(SegmentCommList(SegmentSetToList(set1),SegmentSetToList(set2)))
end

function SegmentSetMinus(set1,set2)
    return SegmentCommSet(set1,SegmentInvertSet(set2))
end

function SegmentPrintSet(set)
    print(SegmentSetToString(set))
end

function SegmentSetToString(set)
    local line = ""
    for i=1,#set do
        if i~=1 then line=line..", " end
        line=line..set[i][1].."-"..set[i][2]
    end
    return line
end

function SegmentSetInSet(set,sub)
    if sub==nil then return true end
    -- Checks if sub is a proper subset of set
    for i=1,#sub do
        if not SegmentRangeInSet(set,sub[i]) then return false end
    end
    return true
end

function SegmentRangeInSet(set,range)
    if range==nil or #range==0 then return true end
    local b=range[1]
    local e=range[2]
    for i=1,#set do
        if b>=set[i][1] and b<=set[i][2] then
            return (e<=set[i][2])
        elseif e<=set[i][1] then return false end
    end
    return false
end

function SegmentSetToBool(set)
    local result={}
    for i=1,structure.GetCount() do
        result[i]=SegmentInSet(set,i)
    end
    return result
end
--- End of Segment Set module

-- Module Find Segment Types
function FindMutablesList()
    local result={}
    for i=1,segCnt2 do if structure.IsMutable(i) then result[#result+1]=i end end
    return result
end
function FindMutables()
    return SegmentListToSet(FindMutablesList())
end

function FindFrozenList()
    local result={}
    for i=1,segCnt2 do if freeze.IsFrozen(i) then result[#result+1]=i end end
    return result
end
function FindFrozen()
    return SegmentListToSet(FindFrozenList())
end

function FindLockedList()
    local result={}
    for i=1,segCnt2 do if structure.IsLocked(i) then result[#result+1]=i end end
    return result
end
function FindLocked()
    return SegmentListToSet(FindLockedList())
end

function FindSelectedList()
    local result={}
    for i=1,segCnt do if selection.IsSelected(i) then result[#result+1]=i end end
    return result
end
function FindSelected()
    return SegmentListToSet(FindSelectedList())
end

function FindAAtypeList(aa)
    local result={}
    for i=1,segCnt2 do
        if structure.GetSecondaryStructure(i)== aa then result[#result+1]=i end
    end
    return result
end
function FindAAtype(aa)
    return SegmentListToSet(FindAAtypeList(aa))
end

function FindAminotype(at) --NOTE: only this one gives a list not a set
    local result={}
    for i=1,segCnt2 do
        if structure.GetAminoAcid(i) == at then result[#result+1]=i end
    end
    return result
end
-- end Module Find Segment Types
-- Module AskSelections
-- 02-05-2012 Timo van der Laan, Free to use for non commercial purposes
function AskForSelections(title,mode)
    local result={{1,structure.GetCount()}} -- All segments
    if mode == nil then mode={} end
    if mode.askloops==nil then mode.askloops=true end
    if mode.asksheets==nil then mode.asksheets=true end
    if mode.askhelixes==nil then mode.askhelixes=true end
    if mode.askligands==nil then mode.askligands=false end
    if mode.askselected==nil then mode.askselected=true end
    if mode.asknonselected==nil then mode.asknonselected=true end
    if mode.askmutateonly==nil then mode.askmutateonly=true end
    if mode.askignorelocks==nil then mode.askignorelocks=true end
    if mode.askignorefrozen==nil then mode.askignorefrozen=true end
    if mode.askranges==nil then mode.askranges=true end
    if mode.defloops==nil then mode.defloops=true end
    if mode.defsheets==nil then mode.defsheets=true end
    if mode.defhelixes==nil then mode.defhelixes=true end
    if mode.defligands==nil then mode.defligands=false end
    if mode.defselected==nil then mode.defselected=false end
    if mode.defnonselected==nil then mode.defnonselected=false end
    if mode.defmutateonly==nil then mode.defmutateonly=false end
    if mode.defignorelocks==nil then mode.defignorelocks=false end
    if mode.defignorefrozen==nil then mode.defignorefrozen=false end
    local Errfound=false
  repeat
    local ask = dialog.CreateDialog(title)
    if Errfound then
        ask.E1=dialog.AddLabel("Try again, ERRORS found, check output box")
        result={{1,structure.GetCount()}} --reset start
        Errfound=false
    end
    if mode.askloops then
        ask.loops = dialog.AddCheckbox("Work on loops",mode.defloops)
    elseif not mode.defloops then
        ask.noloops= dialog.AddLabel("Loops will be auto excluded")
    end
    if mode.askhelixes then
        ask.helixes = dialog.AddCheckbox("Work on helixes",mode.defhelixes)
    elseif not mode.defhelixes then
        ask.nohelixes= dialog.AddLabel("Helixes will be auto excluded")
    end
    if mode.asksheets then
        ask.sheets = dialog.AddCheckbox("Work on sheets",mode.defsheets)
    elseif not mode.defsheets then
        ask.nosheets= dialog.AddLabel("Sheets will be auto excluded")
    end
    if mode.askligands then
        ask.ligands = dialog.AddCheckbox("Work on ligands",mode.defligands)
    elseif not mode.defligands then
        ask.noligands= dialog.AddLabel("Ligands will be auto excluded")
    end
    if mode.askselected then ask.selected = dialog.AddCheckbox("Work only on selected",mode.defselected) end
    if mode.asknonselected then ask.nonselected = dialog.AddCheckbox("Work only on nonselected",mode.defnonselected) end
    if mode.askmutateonly then ask.mutateonly = dialog.AddCheckbox("Work only on mutateonly",mode.defmutateonly) end
    if mode.askignorelocks then
        ask.ignorelocks =dialog.AddCheckbox("Dont work on locked ones",true)
    elseif mode.defignorelocks then
        ask.nolocks=dialog.AddLabel("Locked ones will be auto excluded")
    end
    if mode.askignorefrozen then
        ask.ignorefrozen = dialog.AddCheckbox("Dont work on frozen",true)
    elseif mode.defignorefrozen then
        ask.nofrozen=dialog.AddLabel("Frozen ones will be auto excluded")
    end
    if mode.askranges then
        ask.R1=dialog.AddLabel("Or put in segmentranges. Above selections also count")
        ask.ranges=dialog.AddTextbox("Ranges","")
    end
    ask.OK = dialog.AddButton("OK",1) ask.Cancel = dialog.AddButton("Cancel",0)
    if dialog.Show(ask) > 0 then
        -- We start with all the segments including ligands
        if mode.askloops then mode.defloops=ask.loops.value end
        if not mode.defloops then
            result=SegmentSetMinus(result,FindAAtype("L"))
        end
        if mode.asksheets then mode.defsheets=ask.sheets.value end
        if not mode.defsheets then
            result=SegmentSetMinus(result,FindAAtype("E"))
        end
        if mode.askhelixes then mode.defhelixes=ask.helixes.value end
        if not mode.defhelixes then
            result=SegmentSetMinus(result,FindAAtype("H"))
        end
        if mode.askligands then mode.defligands=ask.ligands.value end
        if not mode.defligands then
            result=SegmentSetMinus(result,FindAAtype("M"))
        end
        if mode.askignorelocks then mode.defignorelocks=ask.ignorelocks.value end
        if mode.defignorelocks then
            result=SegmentSetMinus(result,FindLocked())
        end
        if mode.askignorefrozen then mode.defignorefrozen=ask.ignorefrozen.value end
        if mode.defignorefrozen then
            result=SegmentSetMinus(result,FindFrozen())
        end
        if mode.askselected then mode.defselected=ask.selected.value end
        if mode.defselected then
            result=SegmentCommSet(result,FindSelected())
        end
        if mode.asknonselected then mode.defnonselected=ask.nonselected.value end
        if mode.defnonselected then
            result=SegmentCommSet(result,SegmentInvertSet(FindSelected()))
        end
        if mode.askranges and ask.ranges.value ~= "" then
            local rangetext=ask.ranges.value
            local function Checknums(nums)
                -- Now checking
                if #nums%2 ~= 0 then
                    print("Not an even number of segments found")
                    return false
                end
                for i=1,#nums do
                    if nums[i]==0 or nums[i]>structure.GetCount() then
                        print("Number "..nums[i].." is not a segment")
                        return false
                    end
                end
                return true
            end

            local function ReadSegmentSet(data)
                local nums = {}
                local NoNegatives='%d+' -- - is not part of a number
                local result={}
                for v in string.gfind(data,NoNegatives) do
                    table.insert(nums, tonumber(v))
                end
                if Checknums(nums) then
                    for i=1,#nums/2 do
                        result[i]={nums[2*i-1],nums[2*i]}
                    end
                    result=SegmentCleanSet(result)
                else Errfound=true result={} end
                return result
            end
            local rangelist=ReadSegmentSet(rangetext)
            if not Errfound then
                result=SegmentCommSet(result,rangelist)
            end
        end
    end
  until not Errfound
    return result
end
-- end of module AskSelections

Title="TvdL MakeContactBands "
ContactVersion="1.3.0"
function askContact(title)
    local ask = dialog.CreateDialog(title)
    local inside = true
    print(Title..ContactVersion)
    ask.doselect=dialog.AddCheckbox("Select which sections to band",false)
    ask.insidesel=dialog.AddCheckbox("Only bands inside selection",inside)
    ask.holdstr=dialog.AddSlider("Keep strength:",KeepBandstrength,0.1,10,1)
    ask.holdlen=dialog.AddSlider("Keep len factor:",KeepFactor,0.5,1.5,2)
    ask.pullstr=dialog.AddSlider("Pull strength:",PullBandstrength,0.1,10,1)
    ask.pulllen=dialog.AddSlider("Pull len factor:",PullFactor,0.5,1,2)
    ask.mindist=dialog.AddSlider("Min seg distance:",MinContactGap,2,8,1)
    ask.info=dialog.AddLabel("If not a contact only band if distance is lower")
    ask.maxdist=dialog.AddSlider("Max no contact dist:",MaxNoContactGap,8,20,0)
    ask.printnon=dialog.AddCheckbox("Print all not satisfied contacts",false)
    ask.OK = dialog.AddButton("OK",1)
    ask.Cancel = dialog.AddButton("Cancel",0)
    local result=dialog.Show(ask)
    if result > 0 then
        KeepBandstrength=ask.holdstr.value
        KeepFactor=ask.holdlen.value
        PullBandstrength=ask.pullstr.value
        PullFactor=ask.pulllen.value
        MinContactGap=ask.mindist.value
        MaxNoContactGap=ask.maxdist.value
        PrintNonContacts=ask.printnon.value
        inside=ask.insidesel.value
        if ask.doselect.value then
            if inside then
                AddContactsSelection(AskForSelections(title),KeepBandstrength,PullBandstrength)
            else AddAllContactsSelection(AskForSelections(title),KeepBandstrength,PullBandstrength)
            end
        else AddContactsInRange(1,segCnt2,KeepBandstrength,PullBandstrength)
        end
    end
    print("Max length contact   : ",MaxLengthContact)
    print("Min length no contact: ",MinLengthNoContact)
end

askContact(Title..ContactVersion)