VMD-L Mailing List

From: Tristan Croll (tristan.croll_at_qut.edu.au)
Date: Wed Mar 05 2014 - 05:48:46 CST

Hi all,

The current chirality plugin is only able to handle chiral centers for which all bonded atoms fall within the same residue. This doesn't work for glycans, where the glycosidic bond is both chiral and involves an oxygen atom from an adjacent residue. Furthermore, that oxygen can have multiple different names depending on the particular link.

I've put together a new version of ::chirality::find_chiral_centers (copied below) using a more generalised algorithm. Rather than looking for all the atom names in the chiral center within the one residue, it first selects only the chiral carbon, uses getbonds to get the indices of bonded atoms, and then matches these back to the listed atom names. Multiple atom names are allowable for all but the chiral center itself, and it's entirely agnostic about which residue or segment each atom lies on. It's working for protein and for the four different glycan residue types (as N-linked glycans) I'm working with, and I don't expect that DNA/RNA should be problematic.

Cheers,

Tristan

New entries for glycan residues in the chiralCenters variable. Note that BGLN is renamed from the standard residue name in the CHARMM-36 forcefield, which was over four characters):

  {{BGLN NAG} {{H1 C1 {O1 O2 O3 O4 O6 ND2} O5 C2} {H2 C2 C1 N C3} {H3 C3 O3 C2 C4} {H4 C4 C5 C3 O4} {H5 C5 C4 O5 C6}}}
  {{BMAN BMA} {{H1 C1 {O1 O2 O3 O4 O6} O5 C2} {H2 C2 O2 C1 C3} {H3 C3 O3 C2 C4} {H4 C4 C5 C3 O4} {H5 C5 C4 O5 C6}}}
  {{AMAN MAN} {{H1 C1 {O1 O2 O3 O4 O6} C2 O5} {H2 C2 O2 C1 C3} {H3 C3 O3 C2 C4} {H4 C4 C5 C3 O4} {H5 C5 C4 O5 C6}}}
  {{AFUC FUC} {{H1 C1 {O1 O2 O3 O4 O6} O5 C2} {H2 C2 O2 C1 C3} {H3 C3 C2 O3 C4} {H4 C4 C5 C3 O4} {H5 C5 O5 C4 C6}}}

New ::chirality::find_chiral_centers:

proc ::chirality::find_chiral_centers {seltext molid} {
  variable chiralCenters

  # List returns by this proc. Each element corresponds to an
  # entry in chiralCenters, which contains five lists of atoms,
  # plus a flag if the chiral centers contain hydrogens or not,
  # plus the selected atom names for each chiral center.
  set returnList {}

  # Check each chiral center for correctness
  for {set i 0} {$i < [llength $chiralCenters]} {incr i} {

    set resnamesList [lindex [lindex $chiralCenters $i] 0]
    set atomnamesList [lindex [lindex $chiralCenters $i] 1]
    set sel [atomselect $molid "($seltext) and resname $resnamesList"]

    if {[$sel num] > 0} {

      # If there are alternative conformations, pick the first one in
      # the atom selections below
      set altlocSelector ""
      set altlocList [lsort -unique [$sel get altloc]]
      if {[llength $altlocList] > 1} {
        if {[lindex $altlocList 0] == ""} {
          set altlocSelector "and altloc \"\" [lindex $altlocList 1]"
        } else {
          set altlocSelector "and altloc [lindex $altlocList 0]"
        }
      }

      foreach atomnames $atomnamesList {
        set selatom0 [atomselect $molid "($seltext) and resname $resnamesList and name [lindex $atomnames 1] $altlocSelector"]
        set selatomH [atomselect $molid "($seltext) and resname $resnamesList and name [lindex $atomnames 0] $altlocSelector"]
if {[$selatomH num] > 0} {
set hasH 1
} else {
set hasH 0
}
set atomidH {}
set atomid0 [$selatom0 get index]
set atomid1 {}
set atomid2 {}
set atomid3 {}
for {set sel0index 0} {$sel0index < [llength $atomid0]} {incr sel0index} {
set currentCenterIndex [lindex $atomid0 $sel0index]
set currentCenterSel [atomselect top "index $currentCenterIndex"]
set bondedList [lindex [$currentCenterSel getbonds] 0]
set bondedsel [atomselect $molid "index $bondedList"]
set bondedIndexList [$bondedsel get index]
set bondedNameList [$bondedsel get name]
$bondedsel delete
$currentCenterSel delete

        # Work through the $atomnames set to match to $bondedNameList and thereby assign an index. This is a general
        # approach which allows multiple names for all but the chiral carbon, and doesn't require all atoms to be on the
        # same residue
if {$hasH} {
foreach name [lindex $atomnames 0] {
set bondedNameListPos [lsearch -exact $bondedNameList $name]
if {$bondedNameListPos != -1} {
lappend atomidH [lindex $bondedIndexList $bondedNameListPos]
}
}
}
foreach name [lindex $atomnames 2] {
set bondedNameListPos [lsearch -exact $bondedNameList $name]
if {$bondedNameListPos != -1} {
lappend atomid1 [lindex $bondedIndexList $bondedNameListPos]
}
}

foreach name [lindex $atomnames 3] {
set bondedNameListPos [lsearch -exact $bondedNameList $name]
if {$bondedNameListPos != -1} {
lappend atomid2 [lindex $bondedIndexList $bondedNameListPos]
}
}
foreach name [lindex $atomnames 4] {
set bondedNameListPos [lsearch -exact $bondedNameList $name]
if {$bondedNameListPos != -1} {
lappend atomid3 [lindex $bondedIndexList $bondedNameListPos]
}
}
}

        # If any of the lists is empty, i.e., all chiral centers contain
        # missing atoms, move on to the next chiral center.
        if {([llength $atomid0] == 0 || [llength $atomid1] == 0 ||
             [llength $atomid2] == 0 || [llength $atomid3] == 0) ||
             ($hasH && [llength $atomidH] == 0)} {
          puts "WARNING: Ignoring residue(s) with missing atoms..."
          continue
        }
set selatom1 [atomselect $molid "index $atomid1"]
        set selatom2 [atomselect $molid "index $atomid2"]
        set selatom3 [atomselect $molid "index $atomid3"]

        # Check if index lists have the same length. Otherwise, prune
        # residues with missing atoms from the lists.
        set listLengthOK 1
        set listLength [llength $atomid0]
        if {[llength $atomid1] != $listLength ||
            [llength $atomid2] != $listLength ||
            [llength $atomid3] != $listLength} {
   set listLengthOK 0
        }
        if {$hasH && [llength $atomidH] != $listLength} {
          set listLengthOK 0
        }
        if {!$listLengthOK} {
          puts "WARNING: Ignoring residue(s) with missing atoms..."
        }

        # This check is no longer compatible with the new algorithm
        if {0} {
          set residuesList {}
          lappend residuesList [$selatom0 get residue]
          lappend residuesList [$selatom1 get residue]
          lappend residuesList [$selatom2 get residue]
          lappend residuesList [$selatom3 get residue]

          if {$hasH} {
            lappend residuesList [$selatomH get residue]
          }
          set residuesIntersection [::chirality::list_intersect $residuesList]
          $selatom0 delete
          $selatom1 delete
          $selatom2 delete
          $selatom3 delete

          # If after pruning indices due to missing atoms we end up with an empty
          # index list, silently move on to the next kind of chiral center.
          if {([llength $atomid0] == 0 || [llength $atomid1] == 0 ||
               [llength $atomid3] == 0) || ($hasH && [llength $atomidH] == 0)} {

            continue
          }

          # An empty intersection of residues is VERY likely an error
          if {[llength $residuesIntersection] == 0} {
            error "Empty residuesList, which indicates an error with the dictionary describing chiral centers. Please report this error to the plugin authors."
          }

          set selatom0 [atomselect $molid "(index $atomid0) and (residue $residuesIntersection)"]
          set selatom1 [atomselect $molid "(index $atomid1) and (residue $residuesIntersection)"]
          set selatom2 [atomselect $molid "(index $atomid2) and (residue $residuesIntersection)"]
          set selatom3 [atomselect $molid "(index $atomid3) and (residue $residuesIntersection)"]

          set atomid0 [$selatom0 get index]
          set atomid1 [$selatom1 get index]
          set atomid2 [$selatom2 get index]
          set atomid3 [$selatom3 get index]
          if {$hasH} {
            $selatomH delete
            set selatomH [atomselect $molid "(index $atomidH) and (residue $residuesIntersection)"]
            set atomidH [$selatomH get index]
          }
        }

        $selatom0 delete
        $selatom1 delete
        $selatom2 delete
        $selatom3 delete
        $selatomH delete

        lappend returnList [list $atomidH $atomid0 $atomid1 $atomid2 $atomid3 $hasH $atomnames]

      }
    }

    $sel delete
  }

  return $returnList
}