From: alexandra.marques_at_fc.up.pt
Date: Thu Jan 11 2007 - 13:42:07 CST

Hi again,

I used the pdb file with and without hydrogens and yes, the servers donīt use
hydrogens.
I think I found the problem. If I use only selected residues of the protein to
do the calculation, the sasa values are different and much higher than
if I use
the whole protein residues. For instance, for the first 5 aa:

   Sasa (using the whole protein) Sasa (using only the 5 aa in the selection)
0 124.216300964 148.595306396
1 2.31623363495 127.548957825
2 174.497711182 204.40007019
3 55.33360672 87.0589981079
4 49.2410736084 174.806564331

And the sasa values computed using the whole protein are similar to the ones
reported by the servers. Also, now I can obtain zero values for the buried
residues.
Do you know the reason for this difference? Do I have to use the whole protein
residues to get correct sasa values?

Thanks
Alexandra

Quoting Peter Freddolino <petefred_at_ks.uiuc.edu>:

> 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
>> }
>>
>

-------------------------------------------------------------
A FCUP utiliza o sistema de webmail Horde/IMP (www.horde.org)

Visite: http://www.fc.up.pt/