VMD-L Mailing List

From: Mgr. Lubos Vrbka (lubos.vrbka_at_uochb.cas.cz)
Date: Thu Jul 28 2005 - 05:00:35 CDT

hi,

you might find following script useful. usage is described in the .txt file.

best regards,
lubos 

-- 
.....................................................
Mgr. Lubos Vrbka
Center for Biomolecules and Complex Molecular Systems
Institute of Organic Chemistry and Biochemistry
Academy of Sciences of the Czech Republic
Prague, Czech Republic
.....................................................

SASA All Residue Calculator 1.2
-------------------------------

REQUIREMENTS: VMD Version 1.8.3 or greater

CHANGELOG:
        1.1 fixed some wrong variable names
        1.2 fixed bug in internal representation of SASA values
                worked correctly only for selection 0 - some value, otherwise
                crashed
                script is now using arrays instead of lists
                
DESCRIPTION:
        Script for detailed SASA analysis. Calculates SASA for every residue
        of a given selection (e.g. protein) ad prints it out. It can be used
        to identify residues more/less exposed to the solvent. Works also
        for whole trajectories, were it calculates average SASA value.

HOW IT WORKS:
        Takes one frame after another. Splits given selection to single residues
        and calculates SASA for each of them by using
        measure sasa _radius_ _whole_sel_ -restrict _residue_
        construct.
        This procedure can be very time consuming (360 residue protein, 2000 frames
        trajectory = 7 hours on 3GHz Athlon).
        See the script itself for more information.

PROCEDURES:
        getAllResSASA - main procedure

EXAMPLE USAGE AND OUTPUT:
        usage:
        getAllResSASA "selection" probe_radius <startframe <endframe>>
        Optional startframe/endframe parameters define range of frames from
        the trajectory, that should be analyzed. If omitted, whole trajectory
        is parsed.
        
        example:
        source _sasa_res.tcl
        getAllResSASA "protein" 1.5 1 100
        This performs the analysis on frames 1 to 100 of the trajectory and
        prints average SASA values for every residue of the protein to file
        res_sasa.dat using following format
            residue resname SASA
        i.e.,
            0 GLU 195.54
            1 THR 140.93
            2 THR 66.38
            3 ALA 7.58
            4 LEU 0.51
            5 VAL 0.59
            ...
        
DOWNLOAD FILE:
        _sasa_res.tcl

AUTHOR:
        Lubos Vrbka (shnek (at) tiscali (dot) cz)


# version 1.2
# calculate SASA for every residue of a selection
# usage: getAllResSASA "selection" probe_radius <startframe <endframe>>

proc getAllResSASA { selection radius args } {
  puts "SASA value calculator for $selection residues, probe radius $radius"

  # get number of frames
  set numframes [molinfo top get numframes]
  puts "total number of frames in trajectory: $numframes"
  if {[llength $args] == 0} {
    set startframe 1
    set endframe $numframes
  }
  if {[llength $args] == 1} {
    set startframe [lindex $args 0]
    set endframe $numframes
    set stepframe 1
    if {$startframe <= 0 || $startframe > $numframes} {
      puts "illegal value of startframe, changing to first frame"
      set startframe 1
    }
  }
  if {[llength $args] >= 2} {
    set startframe [lindex $args 0]
    set endframe [lindex $args 1]
    if {$startframe <= 0 || $startframe > $numframes} {
      puts "illegal value of startframe, changing to first frame"
      set startframe 1
    }
    if {$endframe < $startframe || $endframe > $numframes} {
      puts "illegal value of endframe, changing to last frame"
      set endframe $numframes
    }
  }

  set totframes [expr ($endframe - $startframe + 1)]
  puts "analysis will be performed on $totframes frame(s) ($startframe to $endframe)"

  # get list of residues
  set allsel [atomselect top $selection]
  set residlist [lsort -unique -integer [ $allsel get residue ]]

  # resid map for every atom
  set allResid [$allsel get residue]
  # resname map for every atom
  set allResname [$allsel get resname]
  # create resid->resname map
  foreach resID $allResid resNAME $allResname {
    set mapResidResname($resID) $resNAME
  }

  # set sasa for every residue to 0.0
  foreach r $residlist {
    set resSASA($r) 0.0
    # puts "setting 0.0 for residue $r"
}
  
  # now subtract 1 from all frame indexes - numbering starts with 0
  set startframe [expr $startframe - 1]
  set endframe [expr $endframe - 1]

  puts "go and get coffee..."
  for {set i $startframe; set d 1} {$i <= $endframe} {incr i; incr d} {
    # show activity - one dot for every frame
    puts -nonewline "."
    if { [expr $d % 50] == 0 } {
      puts " "
    }
    flush stdout
    
    $allsel frame $i
    $allsel update
  
    foreach r $residlist {
      set sel [atomselect top "residue $r" frame $i]
      set rsasa [measure sasa $radius $allsel -restrict $sel]
      # set user value for this frame
      $sel set user $rsasa
      $sel delete
      set valuesasa 0.0
      foreach {tmp valuesasa} [split [array get resSASA $r] ] break
      # puts "residue $r, sasa: $rsasa, old: $valuesasa"
      # if sasa is above threshold, store it
      #if {$rsasa >= $sasa_limit} {
      set resSASA($r) [ expr { $valuesasa + $rsasa/$totframes} ]
      # }
    }
  }
  
  set fw [open "res_sasa.dat" w]
  foreach r $residlist {
    foreach {tmp resName} [split [array get mapResidResname $r] ] break
    foreach {tmp valuesasa} [split [array get resSASA $r] ] break
    puts $fw "$r $resName $valuesasa"
  }
  close $fw
  puts "done"

  # uncomment following code to change coloring method to "user based"
  #mol modcolor 0 [molinfo top] User
  #mol colupdate 0 [molinfo top] 1
  #mol scaleminmax [molinfo top] 0 auto
}