Plot contact map between aminoacids (alpha carbons) of a protein using AWK and gnuplot

plot-contact-map.awk

#!/usr/bin/awk -f

BEGIN {
  num_atoms = 0;

  if (ARGC < 3) {
    print "Usage: plot-contact-map.awk PDBFILE THRESHOLD" > "/dev/stderr";
    error_exit = 1;
    exit 1;
  }

  threshold = ARGV[2];
  delete ARGV[2];
}

$1 == "ATOM" && $3 == "CA" {
  atoms[num_atoms,type] = $3;
  for (i = 0; i < 3; i++)
    atoms[num_atoms,i] = $(7+i);
  ++num_atoms;
}

function distance(i, j)
{
  s = 0.0;
  for (k = 0; k < 3; k++)
    s += (atoms[i, k] - atoms[j, k])^2;
  return sqrt(s);
}

END {
  if (error_exit)
    exit 1;

  if (num_atoms == 0) {
    print "No coordinates were found in:", FILENAME > "/dev/stderr";
    exit 1;
  }

  gnuplot = "gnuplot -p";

  print "set palette gray" | gnuplot;
  print "unset colorbox" | gnuplot;
  print "set tics out nooffset" | gnuplot;
  print "set xtics 0, 5,", num_atoms | gnuplot;
  print "set ytics 0, 5,", num_atoms | gnuplot;
  print "plot '-' matrix notitle with image" | gnuplot;
  for (i = 0; i < num_atoms; i++) {
    for (j = 0; j < num_atoms; j++) {
      printf "%s ", (distance(i, j) < threshold) ? "0" : "1" | gnuplot;
    }
    print "" | gnuplot;
  }
  print "e\n" | gnuplot;

  close(gnuplot);
}