VMD-L Mailing List

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Tue Dec 11 2007 - 12:05:04 CST

Hi Michel,
this shouldn't be difficult using the approach I described. If I
understand you correctly, you should just be able to keep a counter for
each residue that is updated every frame, eg,

set sel [atomselect top alpha]

set helcounts [list]

foreach ind [$sel get index] {
    lappend helcounts 0
}

for {set i 0} {$i < [molinfo top get numframes]} {incr i} {
    animate goto $i
    display update ui
    mol reanalyze top
    set structs [$sel get structure]

    for {set j 0} {$j < [llength $structs]} {incr j} {
       if {[lindex $structs $j] == "H"} {
          set helcounts [lreplace $helcounts $j $j [expr [lindex
$helcounts $j] + 1] ]
       }
    }
}

(note that this is off the top of my head, so there may be bugs)

Which gives you a list containing the number of frames in which each
residue of your selection is tagged as helical. Is this what you're after?

Best,
Peter

   

animate goto $frame
    display update ui

L. Michel Espinoza-Fonseca wrote:
> Hi Peter,
>
> Thank you for the prompt suggestion. Perhaps I wasn't clear enough,
> but basically what I want to get is an average % of helical content
> per residue over time (i.e., a graph of % helical content vs residue).
> The suggestion you gave me solves the problem if you want to get the
> total helical content of the protein. So I wonder if it'll be
> possible to do that with simple tcl scripting in VMD (I hope so!).
>
> Thanks!
> Michel
>
> On Dec 11, 2007 5:39 PM, Peter Freddolino <petefred_at_ks.uiuc.edu> wrote:
>
>> Hi Michel,
>> I'd recommend instead making a selection with the CA atoms of all the
>> residues you're interested in, and then after going to each frame do
>>
>> mol reanalyze $mol
>> set mystruct [$sel get structure]
>>
>> You can then parse the returned list for the number of "H" entries.
>> Best,
>> Peter
>>
>>
>> L. Michel Espinoza-Fonseca wrote:
>>
>>> Hi all,
>>>
>>> I'm calculating the percentage of helical content of a protein (using
>>> a NAMD-generated trajectory). I've been using the following script:
>>>
>>> set outfile [open protein_SS_16.txt w]
>>> set numframes [molinfo top get numframes]
>>> set sel [atomselect top "resid 16 and name CA"]
>>> for {set frame 0} {$frame < $numframes} {incr frame} {
>>> animate goto $frame
>>> display update ui
>>> vmd_calculate_structure top
>>> $sel frame $frame
>>> $sel update
>>> set helixlist [$sel get alpha_helix]
>>> set helixcount 0
>>> puts $outfile " Alpha helix count: $helixcount "
>>> }
>>> close $outfile
>>>
>>> However, the script calculates the content *ONLY* for one residue at a
>>> time, so I was wondering if you have any ideas on how to modify it so
>>> I can get a table of residue vs count of alpha helix, or something
>>> similar. This is basically to avoid calculating it for every residue
>>> (which will take forever!). If you have alternative ideas, I'd be more
>>> than happy to read them.
>>>
>>> Thanks!
>>> Michel
>>>
>>>