Drawing a distance matrix

 

This example creates a distance matrix made of the distance between the CAs of two residues. The only input value is selection used to find the CAs. The distances are stored in the 2D array ``dist'', which is indexed by the atom indices of the CAs (remember, atom indices start at 0).

The distance matrix is a new graphics molecule. It consists of two triangles for each element of the matrix (to make up a square). The color for each pair is one of the 32 elements in the color scale and is determined so the range of distance values fits linearly between 34 and 66 (excluding 66 itself). The name of the graphics is ``CA_distance(n)'' where ``n'' is the molecule id of the selection used to create the matrix.

  The first graphics command, ``materials off'' is to keep the lights from reflecting off the matrix (otherwise there is too much glare). At the end, the corners of the matrix are labeled in yellow with the resid values of the CAs.

One extra procedure, ``vecdist'', is created. This computes the difference between two vectors of length 3. It is not as general as the normal method (vecnorm [vecsub $v1 $v2]) but is almost twice as fast, which speeds up the overall subroutine by almost 10%. The script is not very fast; after all, for a 234 residue protein, 27495 distance calculations are made and 54756*2 triangles generated. Nearly all of that is done in Tcl. In terms of times, about 1/3 is spent in the distance calculations, another 1/3 in the math to make the triangles, and another 1/3 in the three ``graphics'' calls. The residue 234 example protein took 70 seconds to do everything on a Crimson.

Example usage:

	mol load pdb $env(VMDDIR)/proteins/alanin.pdb
	set sel [atomselect top all]
	ca_distance $sel

# Input: two vectors of length three (v1 and v2)
# Returns: ||v2-v1||
proc vecdist {v1 v2} {
    lassign $v1 x1 x2 x3
    lassign $v2 y1 y2 y3
    return [expr sqrt(pow($x1-$y1,2) + pow($x2-$y2,2) + pow($x3-$y3, 2))]
}


\index{atom!coordinates}
\index{atomselect!command}
# Input: a selection
# Does: finds the CAs in the selection then computes and draws the
#   CA-CA distance grid with colors based on the color scale
# Returns: the id of the new graphics molecule containing the grid
proc ca_distance {main_sel} {
    # get the CA atoms from the selection
    set mol [$main_sel molindex]
    set sel [atomselect $mol "@$main_sel and name CA"]

    # find distances between each pair
    set coords [$sel get {x y z}]
    set max 0
    set list2 [$sel list]

    foreach atom1 $coords id1 [$sel list] {
        foreach atom2 $coords id2 $list2 {
            set dist($id1,$id2) [vecdist $atom2 $atom1]
            set dist($id2,$id1) $dist($id1,$id2)
            set max [max $max $dist($id1,$id2)]
        }
        lvarpop list2
        lvarpop coords
    }


    # draw the pretty graphic
    puts "Distances calculated, now drawing the distance map ..."
    mol load graphics "CA_distance($mol)"
    set gmol [molinfo top]
    # turn material characteristics off
    graphics $gmol materials off
    # i1 and j1 are i+1 and j+1; this speeds up construction of x{01}{01}
    set i 0
    set i1 1
    # preset the scaling factor (31.95 instead of 32 ensures there will be
    # no color 66 (for the max color), which is transparent)
    set scale [expr 31.95 / ($max + 0.)]
    set list2 [$sel list]
    foreach id1 [$sel list] {
       set j 0
       set j1 1
       set x00 "$i $j 0"
       set x10 "$i1 $j 0"
       set x11 "$i1 $j1 0"
       set x01 "$i $j1 0"
       foreach id2 $list2 {
          set col [expr int($scale * $dist($id1,$id2)) + 34]
          graphics $gmol color $col
          graphics $gmol triangle $x00 $x10 $x11
          graphics $gmol triangle $x00 $x11 $x01
          incr j
          incr j1
          set x00 $x01
          set x10 $x11
          set x11 "$i1 $j1 0"
          set x01 "$i $j1 0"
       }
       incr i
       incr i1
    }
    # put some numbers down to give an idea of what is where
    set resids [$sel get resid]
    set num [llength $resids]
    set start [lindex $resids 0]
    set end [lindex $resids [expr $num - 1]]
    graphics $gmol color yellow
    graphics $gmol text "0 0 0" "$start,$start"
    graphics $gmol text "$num 0 0" "$end,$start"
    graphics $gmol text "0 $num 0" "$start,$end"
    graphics $gmol text "$num $num 0" "$end,$end"
    return $gmol
}