Recipe: Local Mutate 4.0

created by spvincent

Profile

• Name: Local Mutate 4.0
• ID: 104257
• Shared with: Public
• Parent: None
• Children:
  • Local Mutate 4.0 (COPY)
• Created on: December 30, 2020 at 00:31 AM UTC
• Updated on: December 30, 2020 at 00:31 AM UTC
• Description:

  See first comment

Code

min_residue = 1
max_residue =  999
n_residues = 0
n_mutations = 0

delta_residue = 9
kMaxDelta = 20

best_score = 0

clash_importance_during_modify = 0.09
packing_importance_during_modify = 1.0
hiding_importance_during_modify = 1.0
pairwise_importance_during_modify = 1.0
backboneHbond_importance_during_modify = 1.0
sidechainHbond_importance_during_modify = 1.0

kWiggleIterations = 15 

modify_count = 10

score_type = 1

loss_from_modify = {}

original_residues = {}
new_residues = {}

scores = {}
ids = {}

start_time = 0

wiggle_sidechains_first = false

--    Amino acids. Hydrophobics first.

single_letter_codes = { "g","a","v","l","i","m","f","w","p","s","t","c","y","n","q","d","e","k","r","h"}
three_letter_codes  = { "Gly","Ala","Val","Leu","Ile","Met","Phe","Trp","Pro","Ser","Thr","Cys","Tyr","Asn","Gln","Asp","Glu","Lys","Arg","His"}

        -- Use of quicksave slots

kOriginalStructureOrNewBest = 10
base_primary_sequence = {}
n_sequence_changes = 0

function r3 ( x )

-- Round to 3 decimal places

    t  = 10 ^ 3
    return math.floor ( x*t + 0.5 ) /  t 
end

function GetScore ()

    if ( score_type == 1 ) then

        score = current.GetScore () 

    end

    return score

end

function aa1_to_aa3 ( k )

  for  j = 1, 20 do
        if ( k == single_letter_codes [  j ] ) then
            return three_letter_codes [  j  ]
        end
    end

    return "Unk"

end

function get_aa3 ( i )

    k = structure.GetAminoAcid (  i )

    for  j = 1, 20 do
        if ( k == single_letter_codes [  j ] ) then
            return three_letter_codes [  j  ]
        end
    end

    return "Unk"

end

function GetBasePrimarySequence ()

    for i = 1, n_residues do
        base_primary_sequence [ i ] = structure.GetAminoAcid ( i )
    end

end

function GetNumberOfPrimarySequenceChanges ()

    n_changes = 0
    for i = 1, n_residues do
        if ( base_primary_sequence [ i ] ~= structure.GetAminoAcid ( i ) ) then
            n_changes = n_changes + 1
        end
    end

    return n_changes
end


function ShellSort ()

-- Adapted from Numerical Recipes in C

    inc =  1
    
    repeat 
        inc = inc * 3 + 1
    until  inc > n_residues

    repeat

        inc = inc  / 3 
        inc = inc - inc % 1

        for i = inc + 1 ,  n_residues do

            v =  scores [  i  ]
            w = ids [ i ]

            j = i

            flag = false 

            while (  flag == false and scores [  j - inc ] > v ) do
                scores [  j  ] = scores [  j  - inc ]
                ids [ j ] = ids [ j - inc ]

                j = j - inc
                
                if ( j <= inc ) then
                    flag = true
                end

            end

            scores [ j ]  = v
            ids [ j ] = w
        end

    until  inc  <= 1

end


function PostModify ( i )

    behavior.SetClashImportance ( 1.0 )
    behavior.SetPackingImportance ( 1.0 )
    behavior.SetHidingImportance ( 1.0 )
    behavior.SetPairwiseImportance ( 1.0 )
    behavior.SetBackboneHBondImportance ( 1.0 )
    behavior.SetSidechainHBondImportance ( 1.0 )
 
    selection.SelectAll ()
    if ( wiggle_sidechains_first == true ) then
        structure.WiggleSelected ( 1 , false , true )    -- sidechains
    end
 
    score_before_wiggle = GetScore ()
   
    structure.WiggleSelected ( kWiggleIterations )        
    score_after_wiggle = GetScore ()

end


function MutateResidue ( i  )

    local	j
    local 	k

    save.Quickload ( kOriginalStructureOrNewBest )

     for  j = 1 , n_residues do
        scores [ j ] = structure.GetDistance ( i , j )
        ids [ j ] = j
    end 

    ShellSort ()

    selection.DeselectAll ()
    for j = 1 , modify_count do
        selection.Select ( ids [ j ] )
    end

    behavior.SetClashImportance ( clash_importance_during_modify )
    behavior.SetPackingImportance ( packing_importance_during_modify )
    behavior.SetHidingImportance ( hiding_importance_during_modify )
    behavior.SetPairwiseImportance ( pairwise_importance_during_modify )
    behavior.SetBackboneHBondImportance ( backboneHbond_importance_during_modify )
    behavior.SetSidechainHBondImportance ( sidechainHbond_importance_during_modify )
 
    structure.MutateSidechainsSelected ( 1 )  
 
    score = GetScore ()
    table.insert ( loss_from_modify , best_score - score )

    n_mutations = GetNumberOfPrimarySequenceChanges ()

 --   print ( "Hid " .. r3 ( scores [ modify_count ] ) .. " lim " .. r3 ( best_score - score ) .. "  nc " .. n_mutations )
    print (  "      lom " .. r3 ( best_score - score ) .. "  nc " .. n_mutations )


    for k = 1, n_residues do
        if ( base_primary_sequence [ k ] ~= structure.GetAminoAcid ( k ) ) then
            print ( "         " .. k .. "  "  .. aa1_to_aa3 (  base_primary_sequence [ k ] ) .. " -> " .. aa1_to_aa3 ( structure.GetAminoAcid ( k ) ) )
        end
    end


    if ( math.abs ( score - best_score ) > 0.5 ) then
        PostModify ( i )

        score = GetScore ()
        if ( score > best_score ) then

            save.Quicksave ( kOriginalStructureOrNewBest )
            GetBasePrimarySequence ()
            best_score = score
           print ( "Improvement to " .. r3 ( best_score ) )
       end

    end

end

function MutateAll ()

    n = 1

    for start_idx = 1 , delta_residue do

        i = min_residue - 1 + start_idx
        
        while ( i <= max_residue ) do

            print ( get_aa3 (i) .. " " .. i .. " (" .. n .. "/" .. max_residue - min_residue + 1 .. ")" )

            MutateResidue ( i )
            i = i + delta_residue
            n = n + 1
        end

    end

end

function GetOptionalParameters ()

    local dlog = dialog.CreateDialog ( "Options" )

    dlog.pidm = dialog.AddSlider ( "Packing importance" , packing_importance_during_modify , 0 , 3.0 , 2 )  
    dlog.hidm = dialog.AddSlider ( "Hiding importance" , hiding_importance_during_modify , 0 , 3.0 , 2 )  
    dlog.widm = dialog.AddSlider ( "Pairwise importance" , pairwise_importance_during_modify , 0 , 3.0 , 2 )  
    dlog.bidm = dialog.AddSlider ( "Backbone hbond importance" , backboneHbond_importance_during_modify , 0 , 3.0 , 2 )  
    dlog.sidm = dialog.AddSlider ( "Sidechain hbond importance" , sidechainHbond_importance_during_modify , 0 , 3.0 , 2 )  

    dlog.ok = dialog.AddButton ( "OK" , 1 )
    dlog.cancel = dialog.AddButton ( "Cancel" , 0 )

    if ( dialog.Show ( dlog ) > 0 ) then
    
        packing_importance_during_modify = dlog.pidm.value
        hiding_importance_during_modify = dlog.hidm.value
        pairwise_importance_during_modify = dlog.widm.value
        backboneHbond_importance_during_modify = dlog.bidm.value
        sidechainHbond_importance_during_modify = dlog.sidm.value
 
        return true
      else
        return false
    end

end

function GetParameters ()

    local dlog = dialog.CreateDialog ( "Local Mutate 4.0" )

    dlog.min_residue = dialog.AddSlider ( "Min residue" , min_residue , 1 , n_residues , 0 )  
    dlog.max_residue = dialog.AddSlider ( "Max residue" , max_residue , 1 , n_residues , 0 )  
    dlog.delta_residue = dialog.AddSlider ( "Delta residue" , delta_residue  , 1 , kMaxDelta , 0 )  
    dlog.cidm = dialog.AddSlider ( "Clash importance" , clash_importance_during_modify , 0 , 1.0 , 2 )  
    dlog.modify_count = dialog.AddSlider ( "Modify count" , modify_count , 1 , 40 , 0 )
   
    dlog.ok = dialog.AddButton ( "OK" , 1 )
    dlog.cancel = dialog.AddButton ( "Cancel" , 0 )
    dlog.other_options = dialog.AddButton ( "Other options" , 2 )

    return_code = dialog.Show ( dlog )

    if ( return_code > 0 ) then
    
        min_residue = dlog.min_residue.value
        max_residue = dlog.max_residue.value
        delta_residue = dlog.delta_residue.value
        clash_importance_during_modify = dlog.cidm.value
        modify_count = dlog.modify_count.value
 
    end
    
    return return_code

end

function main ()

    band.DisableAll ()
    n_residues = structure.GetCount ()
 
    -- Trim off locked terminal sequences as a UI convenience
    
    min_residue = 1
    while ( ( structure.IsLocked ( min_residue ) == true ) and ( min_residue <= n_residues ) ) do
      min_residue = min_residue + 1
    end
    
    max_residue = n_residues
    while ( ( structure.IsLocked ( max_residue ) == true ) and ( max_residue >= 1 ) ) do
      max_residue = max_residue - 1
    end  

    save.Quicksave ( kOriginalStructureOrNewBest )
    GetBasePrimarySequence ()

    repeat
       dialog_code = GetParameters ()
       if ( dialog_code == 2 ) then
           GetOptionalParameters ()
       end
    until ( dialog_code < 2 )

    if ( dialog_code == 0 ) then
        return
    end

    print ( "Local Mutate 4.0" )
 
    max_residue = math.min ( n_residues , max_residue )
    min_residue = math.max ( 1 , min_residue )

    print  ( "Mutate range " .. min_residue .. " to " .. max_residue )
    if ( max_residue < min_residue ) then
        print ( "Mutate range error." )
        return
    end
 
    print  ( "Delta " .. delta_residue )
    
    print ( "Clashing importance " .. r3 ( clash_importance_during_modify ) )
    print ( "Packing importance " .. r3 ( packing_importance_during_modify ) )
    print ( "Hiding importance " .. r3 ( hiding_importance_during_modify ) )
    print ( "Pairwise importance " .. r3 ( pairwise_importance_during_modify ) )
    print ( "Backbone Hbond importance " .. r3 ( backboneHbond_importance_during_modify ) )
    print ( "Sidechain Hbond importance " .. r3 ( sidechainHbond_importance_during_modify ) )
    
    print  ( "CI on modify " .. r3 ( clash_importance_during_modify ) )
    print  ( "Modify count " .. modify_count )
    print  ( "" )

    best_score = GetScore ()
    print (  "Start  score : " .. r3 ( best_score ) )

    for i = 1 , n_residues do
        original_residues [ i ] = structure.GetAminoAcid ( i )
    end

    if ( ( math.abs ( packing_importance_during_modify ) - 1 > 1e-3 ) or 
         ( math.abs ( hiding_importance_during_modify ) - 1 > 1e-3 ) or 
         ( math.abs ( pairwise_importance_during_modify ) - 1 > 1e-3 ) or 
         ( math.abs ( backboneHbond_importance_during_modify ) - 1 > 1e-3 ) or 
         ( math.abs ( sidechainHbond_importance_during_modify ) - 1 > 1e-3 ) ) then 
       print ( "" )
       print ( "MAKE SURE Recipe Score Modding IS CHECKED" )
       print ( "" )
    end
    
    start_time = os.time ()

    MutateAll ()

    cleanup ()

end

function cleanup ()
   
    print ( "Cleaning up" )
   
    behavior.SetClashImportance ( 1.0 )
    behavior.SetPackingImportance ( 1.0 )
    behavior.SetHidingImportance ( 1.0 )
    behavior.SetPairwiseImportance ( 1.0 )
    behavior.SetBackboneHBondImportance ( 1.0 )
    behavior.SetSidechainHBondImportance ( 1.0 )

    save.Quickload ( kOriginalStructureOrNewBest )

    selection.SelectAll ()
    band.EnableAll ()

    if ( #loss_from_modify > 1 ) then
        table.sort ( loss_from_modify )
        mid_pt = math.ceil ( #loss_from_modify / 2 )
        print ( "Median loss from modify " .. r3 ( loss_from_modify [ mid_pt ] ) )
     end

    for i = 1 , n_residues do
        new_residues [ i ] = structure.GetAminoAcid ( i )
    end

    for i = 1 , n_residues do
        if ( new_residues [ i ] ~= original_residues [ i ] ) then
            print ( i .. " : " .. aa1_to_aa3 ( original_residues  [ i ]  ) .. " -> " .. aa1_to_aa3 ( new_residues  [ i ] ) )
        end
    end

    end_time = os.time ()
    print ( "Elapsed time : " .. os.difftime ( end_time , start_time ) .. " secs" )

end

--main ()
xpcall ( main , cleanup )

Comments

[spvincent]

Not that different from 3.0: the main change is to allow control over various part score importances like hiding.

Algorithm outline

For each amino acid
Set the Clash Importance to something low
Select everything in the immediate neighbourhood of that residue
Mutate
Set the CI to 1 and select the whole protein
Wiggle
Record any improvement
end

Dialog settings

Min residue, Max residue

The range to be mutated : by default the whole protein with locked ends trimmed off.

Delta residue

Processing each residue in turn is inefficient: the selection in each case will be largely duplicated. This setting indicates the difference between successive residues. For example, if the min and max residues above are 1 and 100 respectively and the delta value is 5, then residues will be processed in the order 1,5,10,15,etc…..2,6,11 and so on.

Clash importance

The value during the initial mutate. Something that reduces the energy score by about 200 seems about right (the median reduction is printed out at the end of the script)

Modify count

Determines how many amino acids will be mutated at each step. They will be the closest (as measured by the distance between alpha carbons) to the residue being processed. Fewer means faster processing and will enable finer control over the mutating process: for example it may be desirable to minimize the chance of a mutation disrupting any existing Hydrogen Bond Networks.

Packing importance
Hiding importance
Pairwise importance
Backbone Hbond importance
Sidechain Hbond importance

Sliders added to control these values (in Other Options). If you change them make sure Behavior->Recipe Score Modding is checked. Have you to get a handle on how useful changing these will be.