VMD-L Mailing List

From: Ritu Arora (ritu.arora_at_mail.concordia.ca)
Date: Mon Oct 26 2020 - 15:52:33 CDT

Hi Giacomo,

This is the full contactFreq.tcl script:

proc contactFreq { {sel1} {sel2} {percent 0} {outFile stdout} {mol top} } {
  if { $outFile != "stdout" } {
     set outFile [open $outFile w]
  }
  puts $outFile "[clock format [clock scan now]] Search started."

  set allAtoms {}
  set allCount {}
  set numberOfFrames [molinfo $mol get numframes]
  for { set i 0 } { $i < $numberOfFrames } { incr i } {
      molinfo $mol set frame $i

      set frameCount {}
      set frameAtoms [atomselect $mol "$sel1 and noh and within 4.5 of ($sel2 and noh)"]
      foreach {segid} [$frameAtoms get segid] {resname} [$frameAtoms get resname] {resid} [$frameAtoms get resid] {name} [$frameAtoms get name] {
 set atom [list $resid $resname $segid]
 if {[lsearch $frameCount $atom] != -1} continue
 lappend frameCount $atom
 set loc [lsearch $allAtoms $atom]
          if { $loc == -1 } {
             lappend allAtoms $atom
             lappend allCount 1
          } else {
             lset allCount $loc [expr { [lindex $allCount $loc] + 1 } ]
          }
     }
     $frameAtoms delete
  }

  puts $outFile "[clock format [clock scan now]] Search finished."

  puts $outFile "Find interactions:"
  puts $outFile "Residue \t\tfraction"
  #print count after sorting
  set outData {}
  foreach { a } $allAtoms { c } $allCount {
      lappend outData [concat $c $a]
  }
  foreach { data } [lsort -integer -index 1 $outData] {
      set c [lindex $data 0]
      set fraction [expr { 100*$c/($numberOfFrames+0.0) }]
      if { $fraction >= $percent } {
puts $outFile [format "%s-%s-%s \t\t %.2f%%" [lindex $data 3] [lindex $data 2] [lindex $data 1] $fraction]
      }
      #set beta according to the fraction, this is optional
      set atom [atomselect $mol "segid [lindex $data 3] and resname [lindex $data 2] and resid [lindex $data 1]"]
      $atom set beta $fraction
      $atom delete
  }

  if { $outFile != "stdout" } {
      close $outFile
  }

}

Thanks,
Ritu

________________________________
From: Giacomo Fiorin <giacomo.fiorin_at_gmail.com>
Sent: Monday, October 26, 2020 4:49 PM
To: Ritu Arora <ritu.arora_at_mail.concordia.ca>
Cc: vmd-l_at_ks.uiuc.edu <vmd-l_at_ks.uiuc.edu>
Subject: Re: how does contactFreq.tcl work?

Hello Ritu, I am not familiar with the script in question. I'm CC'ing the VMD mailing list below to see if anybody has an idea.

Giacomo

On Mon, Oct 26, 2020 at 4:44 PM Ritu Arora <ritu.arora_at_mail.concordia.ca<mailto:ritu.arora_at_mail.concordia.ca>> wrote:
Hi Giacomo,

Can I please ask you a small query related to the contactFreq.tcl script? How the following highlighted portions of the script is working:

 lappend frameCount $atom
 set loc [lsearch $allAtoms $atom]
          if { $loc == -1 } {
             lappend allAtoms $atom
             lappend allCount 1
          } else {
             lset allCount $loc [expr { [lindex $allCount $loc] + 1 } ]
          }
     }
     $frameAtoms delete
  }

  puts $outFile "[clock format [clock scan now]] Search finished."

  puts $outFile "Find interactions:"
  puts $outFile "Residue \t\tfraction"
  #print count after sorting
  set outData {}
  foreach { a } $allAtoms { c } $allCount {
      lappend outData [concat $c $a]
  }
  foreach { data } [lsort -integer -index 1 $outData] {
      set c [lindex $data 0]
      set fraction [expr { 100*$c/($numberOfFrames+0.0) }]

Sorry for my ignorance, but I fail to decode the script, how does the script count sustained contacts? The script counts the number of atoms of a selection within 4.5 of another selection. So, if there are more than one atom of a residue in a frame, it keeps on adding those atoms..not sure?

I look forward to your reply.

Thank you.
Best,
Ritu