VMD-L Mailing List

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Tue Jan 09 2007 - 10:02:00 CST

Hi Alexandra,
as a first check, could you verify that you're using input that includes
explicit hydrogens? The servers you mentioned ignore them (I believe),
but VMD doesn't. As described in the VMD user's manual
(http://www.ks.uiuc.edu/Research/vmd/current/ug/node124.html), the sasa
is calculated by extending the radius of each atom in the target
selection by the solvent radius and calculating the area of the
resultant spheres which is not overlapping another atom. By this
definition even buried residues are likely to have some amount of
exposed surface (since they do not fit together perfectly); the question
of whether or not a given residue is actually in a solvent-accessible
crevice or not is not addressed here (this is part of the purpose of the
-restrict option).

Peter

alexandra.marques_at_fc.up.pt wrote:
> Hi
>
> I am using a script (_sasa_res.tcl) that I found in the mailing list archives to
> calculate the sasa by residue. However the values obtained are much higher than
> expected and buried residues do not have zero values. For instance, I compared
> the values obtained for a pdb file with those reported by the PQS or EBI PISA
> servers and they are very different.
> I am sending the script in attachment. Could someone please take a look to it
> and/or suggest me another way to calculate sasa by residue? What is the method
> used by vmd to calculate SASA?
>
> Thanks in advance
> Alexandra
>
> -------------------------------------------------------------
> A FCUP utiliza o sistema de webmail Horde/IMP (www.horde.org)
>
> Visite: http://www.fc.up.pt/
>
> ------------------------------------------------------------------------
>
> # 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
> }
>