Subject: Summary PDB superposition and RMS

From: "l.r.cooper" <[email-masked]>
Date: Fri, 14 May 1999 11:36:20 +0100

Dear Lee,

This question had also been posted several years
ago, and I enclose the response received in that
time, together with a useful piece of software.
Regards.

Jose I. Garcia

From: 
       "Rino Ragno (uniroma1.it)" <[email-masked]>
    To: 
       "l.r.cooper" <[email-masked]>

Hi,

PROFIT 17 does what you are looking for, have a look at:

www.biochem.ucl.ac.uk/~martin/programs/index.html

Ciao

From: 
       [email-masked] (Paul Davis)
    To: 
       [email-masked]

Dear Lee,

Unless you are looking at different conformations of the
same molecule, this is not a trivial problem because you
have to decide which atoms in each molecule have to be
equivalenced.

Here are two suggestions.

1) ProFit by Andrew Martin is an excellent, flexible program.

http://www.biochem.ucl.ac.uk/~martin/programs/index.html#profit

2) The 'Motif' routine in WhatIf (by Gert Vriend) is
especially worth looking at; it offers automatic alignment
of molecules in three dimensions. The program suggests which
atoms should be equivalenced. WhatIf is a large and complex
modelling and structural analysis package, and so may not be
exactly what you were asking for; however you can write scripts
to routinely perform a given set of tasks.

http://www.sander.embl-heidelberg.de/whatif/

ProFit probably fits the bill.

    http://www.biochem.ucl.ac.uk/~martin/programs/index.html

cheers,

Malcolm
--
Malcolm Gillies <[email-masked]>
PhD student, computational medicinal chemistry
Dept Medicinal Chemistry, Faculty of Pharmacy,
Utrecht University, The Netherlands

There are many programs try

 http://www.biochem.ucl.ac.uk/~martin/swreg.html = FroFit v1.8

Miri

From: 
       "Mauro A. Cremonini" <[email-masked]>
    To: 
       "l.r.cooper" <[email-masked]>

Hi!
I did it once in fortran. I used the orthogonal procrustes algorithm by
Golub and Pereira.
Unfortunately I haven't the code anymore, but it was some 30 lines of
fortran code and one call to the lapack library.
If you want I can fax you 4 pages with the original reference.
Hope it helps,
Mauro

Try SwissPDB viewer

Claude Aflalo ######################## Phones: Office 972-7-6472118 ####
Dept. of Life Sciences                          Lab   972-7-6472119
Ben Gurion University of the Negev              Fax   972-7-6472890
P.O.Box 653                             email: [email-masked]
Beer Sheva 84105 Israel       URL: http://www.bgu.ac.il/life/aflalo.html

Hi Lee,

I suggest you use my ProFit program:
http://www.biochem.ucl.ac.uk/~martin/programs/

Hi, I recommand you use:
(1) the superpose module in Sybly
(2) TINKER
(3) Modeller

both tinker and modeller are public available without any
charge for academic coummunity, and Tinker also offers the
F77 (?) source codes, and the c codes translated by f2c.

Good luck.
Dayong

    I don't know of any c code that will do that 
but I know that you can do it with MSI's Cerius.

Nathalie Godbout                           |
| University of Illinois at Urbana-Champaign |
| 600 S. Goodwin Av. Box 23.6                |
| Urbana, IL 61801     

-- 
°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`
//=\Mr. Lee R. Cooper Bsc(Hons) MSc, Parallel, Emergent//=\evolution,
\=//& Distributed Architechture Laboratory, Dept. of   \=//creationists  
 //=\Comp. Science, University of Reading, Whiteknights //=\have
argued,  
 \=//PoBox 225, Reading, Berkshire. RG6 6AY (UK)        \=/was as likely
  //=\Tel +44 (0)118 987 5123 x7645 (0915-1700)|Fax	 //=to have made 
  \=//Tel +44 (0)118 931 8188 x7817 (1900-2300)|0118 975 \=//man as a 
   //=\http://www.rdg.ac.uk/~ssr97lrc/		   1994   //=\monkey...,
   °¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`´°¤ø,¸¸,ø¤°`
randomly typing the bible without error. with natural selection acting
as an eraser of mistakes; assuming thirteen mistakes for every correct
character & a rate of 1 keystroke per second, the 6 million characters
of the holy script could be typed in 13x6e6 s, or <3 yrs.
Dear net people,

        With the recent interest in molecular super-positioning I
thought I'd post a copy of my program that I sent to Mike Smith (who
originally wrote to the list wanting a program to overlay two
molecules).

	The program takes two PDB files and two lists of atom numbers
and produces a translated/rotated version of the second PDB file so the
specified atoms are overlayed as closely as possible.

        I'm sending the source with this message, hopefully it should
get through uncorrupted. If there are any problems due to more than 80
characters per line (my text editor is set up with about 150 columns)
I could ftp it to you or something.

        There are 3 source files, "overlay.c", "utils.c" and "utils.h".
After you've chopped them out of this file, make the program by doing

cc -o overlay overlay.c utils.c

	You'll also need a "-lm" flag on SGI's as for some reason the
maths part of the standard C library isn't automatically linked in.

        Then to overlay the two PDB files you do

overlay <pdb1> <atom_list1> <pdb2> <atom_list2>

        This will produce a file called "<pdb2>.super". And tell you
the RMS distance between the atoms you specified to be overlayed.

        The atom list files should have one atom id per line. The
atom id can be just its index number, or you can specify the chain
too. The format of each line is...

<atom index>[<white space><chain id>]

        where things in brackets [] are optional.

eg.
118
4560
etc

or

118     A
118     B
4560    B
etc

        Obviously there must be the same number of atoms in each list.

        I hope this is useful. If there are any problems or queries send me a
message.

        yours,

	Simon Kilvington ([email-masked])

The source...

=================================================================
=================================================================

/*
   overlay.c

   Simon Kilvington, 1994
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "utils.h"

#define MAXFILENAMELEN		256
#define PDBLINELEN		82		/* length of a line in a pdb file, inc space for "\n\0" terminator */
#define INDEXCOLUMN		6		/* position of atom number in ATOM records */
#define CHAINCOLUMN		21		/* position of chain id in ATOM records */
#define XPOSCOLUMN		30		/* position of x coord in ATOM records */
#define YPOSCOLUMN		38
#define ZPOSCOLUMN		46
#define COORDLEN		8		/* width of coord fields */
#define COORDDPS		3		/* number of decimal places for coords */
#define SUPERPOSTOLERANCE	0.000000175	/* when superpositioning give up when tan(rotation needed) is less than this, ~tan(0.00001degs) */
#define MAXSUPERPOSITS		100000		/* or we need more than this number of iterations */

typedef struct
{
   char   pdb[PDBLINELEN];
   int    index;				/* pdb atom number */
   char   chain;				/* chain id, resolves atom number clashes */
   vector pos;
} atomdetails;

typedef struct
{
   int         natoms;
   atomdetails atom[1];				/* array extends to atom[natoms-1] */
} fragdetails;

#define FRAGHDRSIZE	(sizeof(fragdetails) - sizeof(atomdetails))

fragdetails *readpdbdata(char *);
void writepdbdata(char *, fragdetails *);

int *readatomlist(char *, fragdetails *, int *);
int findatomindex(fragdetails *, int, char);

void superpositionfrags(int, fragdetails *, int *, fragdetails *, int *);
void calcsuperposmatrix(int, fragdetails *, int *, fragdetails *, int *, matrix *, vector *, vector *);

/*
   main
   params should be <pdb file1> <atom file1> <pdb file2> <atom file2>
   the <atom files> should be lists of atoms (and, optionally, chain id's) to be overlapped, there must be the same number in each file
   produces a pdb file called "<pdb file2>.super"
*/

int
main(int argc, char **argv)
{
   fragdetails *pdb1, *pdb2;
   int *atom1, *atom2, natoms1, natoms2, a;
   float rms;
   vector dr;
   char filename[MAXFILENAMELEN];

   if(argc != 5)
   {
      printf("Useage: overlay <pdb1> <atom list1> <pdb2> <atom list2>\n");
      return 1;
   }

/* read in the 	files */
   if((pdb1 = readpdbdata(argv[1])) == NULL)
      return 1;
   if((atom1 = readatomlist(argv[2], pdb1, &natoms1)) == NULL)
      return 1;
   if((pdb2 = readpdbdata(argv[3])) == NULL)
      return 1;
   if((atom2 = readatomlist(argv[4], pdb2, &natoms2)) == NULL)
      return 1;
   if(natoms1 != natoms2)
   {
      printf("Both atom lists must have the same number of atoms\n");
      return 1;
   }

   printf("Moving \"%s\"\n", argv[3]);
   superpositionfrags(natoms1, pdb1, atom1, pdb2, atom2);

   strcpy(filename, argv[3]);
   strcat(filename, ".super");
   writepdbdata(filename, pdb2);

/* report the rms distance between atom1[] and atom2[] */
   rms = 0.0;
   for(a=0; a<natoms1; a++)
      rms += vector_len2(vector_sub(&pdb1->atom[atom1[a]].pos, &pdb2->atom[atom2[a]].pos, &dr));
   rms = sqrt(rms / natoms1);
   printf("RMS distance between atoms to overlap: %.3f\n", rms);
   printf("Number of atom pairs: %d\n", natoms1);

   free(pdb1);
   free(atom1);
   free(pdb2);
   free(atom2);

   return 0;
}

fragdetails *
readpdbdata(char *filename)
{
   fragdetails *frag, *newptr;
   atomdetails *aptr;
   FILE *file;
   char buffer[PDBLINELEN];
   BOOL err;

   if((file = fopen(filename, "r")) == NULL)
   {
      printf("Unable to read \"%s\"\n", filename);
      return NULL;
   }

   printf("Reading \"%s\"\n", filename);

   frag = malloc(FRAGHDRSIZE);
   frag->natoms = 0;

   err = FALSE;
   while(!feof(file) && !err)
   {
   /* read in the line */
      fgets(buffer, PDBLINELEN, file);
   /* if this is an ATOM record add it to the growing fragment */
      if(strncmp(buffer, "ATOM  ", 6) != 0)
	 continue;
      frag->natoms++;
      newptr = realloc(frag, FRAGHDRSIZE + (frag->natoms * sizeof(atomdetails)));
      if(err = (newptr == NULL))
	 continue;
      frag = newptr;
      aptr = &frag->atom[frag->natoms-1];
      strncpy(aptr->pdb, buffer, PDBLINELEN);
      aptr->index = atoi(&buffer[INDEXCOLUMN]);
      aptr->chain = buffer[CHAINCOLUMN];
      aptr->pos.x = atof(&buffer[XPOSCOLUMN]);
      aptr->pos.y = atof(&buffer[YPOSCOLUMN]);
      aptr->pos.z = atof(&buffer[ZPOSCOLUMN]);
   }

   fclose(file);

   if(err)
   {
      printf("No memory to read \"%s\"\n", filename);
      free(frag);
   }

   return (err) ? NULL : frag;
}

void
writepdbdata(char *filename, fragdetails *frag)
{
   FILE *file;
   int a;
   atomdetails *aptr;

   if((file = fopen(filename, "w")) == NULL)
   {
      printf("Unable to write to \"%s\"\n", filename);
      return;
   }

   printf("Writing \"%s\"\n", filename);

   for(a=0; a<frag->natoms; a++)
   {
      aptr = &frag->atom[a];
      sprintf(&aptr->pdb[XPOSCOLUMN], "%*.*f", COORDLEN, COORDDPS, aptr->pos.x);
      sprintf(&aptr->pdb[YPOSCOLUMN], "%*.*f", COORDLEN, COORDDPS, aptr->pos.y);
      sprintf(&aptr->pdb[ZPOSCOLUMN], "%*.*f", COORDLEN, COORDDPS, aptr->pos.z);
      fprintf(file, "%s\n", aptr->pdb);
   }

   fclose(file);

   return;
}

int *
readatomlist(char *filename, fragdetails *frag, int *natoms)
{
   int *list, *newptr, index, bpos;
   char chain, buffer[PDBLINELEN];
   FILE *file;
   BOOL err;

   *natoms = 0;

   if((file = fopen(filename, "r")) == NULL)
   {
      printf("Unable to read \"%s\"\n", filename);
      return NULL;
   }

   printf("Reading \"%s\"\n", filename);

   list = malloc(sizeof(int));	/* so we can use realloc from the start */

   err = FALSE;
   while(!feof(file) && !err)
   {
      (*natoms)++;
      newptr = realloc(list, (*natoms) * sizeof(int));
      if(err = (newptr == NULL))
      {
	 printf("No memory read \"%s\"\n", filename);
	 continue;
      }
      list = newptr;
      fgets(buffer, PDBLINELEN, file);
      sscanf(buffer, "%d%n", &index, &bpos);
      while(buffer[bpos] != '\0' && isspace(buffer[bpos]))
	  bpos++;
      chain = (buffer[bpos] == '\0') ? ' ' : buffer[bpos];
      if(err = ((list[(*natoms)-1] = findatomindex(frag, index, chain)) == -1))
	 printf("Error while reading \"%s\"; PDB atom with index %d%c does not exist\n", filename, index, chain);
   }

   fclose(file);

   if(err)
      free(list);

   return (err) ? NULL : list;
}

int
findatomindex(fragdetails *frag, int index, char chain)
{
   int atom;
   BOOL match;

   atom = 0;
   match = FALSE;
   while(!match && atom < frag->natoms)
   {
      match = (frag->atom[atom].index == index && frag->atom[atom].chain == chain);
      if(!match)
         atom++;
   }

   return (atom == frag->natoms) ? -1 : atom;
}

/*
   superpositionfrags
   this moves yfrag so the specified atoms in it and xfrag overlap as closely as possible
   the xatom and yatom arrays should be natoms long and hold the numbers of the atoms that are required to be overlappped
   uses the Ferro-Hermans algorithm, Acta Cryst. 33 (1977) p345-347
*/

void
superpositionfrags(int natoms, fragdetails *xfrag, int *xatom, fragdetails *yfrag, int *yatom)
{
   matrix spos, ypos, result;
   vector cgx, cgy;
   int a;
   atomdetails *aptr;

   matrix_alloc(&spos, 3, 3, TRUE);
   matrix_alloc(&ypos, 3, 1, TRUE);

   calcsuperposmatrix(natoms, xfrag, xatom, yfrag, yatom, &spos, &cgx, &cgy);

   for(a=0; a<yfrag->natoms; a++)
   {
      aptr = &yfrag->atom[a];
      matrix_setelement(&ypos, 0, 0, aptr->pos.x - cgy.x);
      matrix_setelement(&ypos, 1, 0, aptr->pos.y - cgy.y);
      matrix_setelement(&ypos, 2, 0, aptr->pos.z - cgy.z);
      matrix_mult(&spos, &ypos, &result);
      aptr->pos.x = cgx.x + matrix_getelement(&result, 0, 0);
      aptr->pos.y = cgx.y + matrix_getelement(&result, 1, 0);
      aptr->pos.z = cgx.z + matrix_getelement(&result, 2, 0);
      matrix_forget(&result);
   }

   matrix_forget(&ypos);
   matrix_forget(&spos);

   return;
}

/*
   calcsuperposmatrix
   as superpositionfrags, but yfrag is unmoved by the experience
   on entry spos should be a matrix_alloc'd 3x3 matrix
   on exit spos is the superposition matrix, cgx is the centre of gravity of xatom, and cgy is the centre of gravity of yatom
   to superposition the two frags yfrag should be moved to "cgx + spos(y - cgy)" where y are the original coords of yfrag
   !!! this routine depends on the fact that a "vector" is composed of 3 consecutive floats that represent x, y and z coords !!!
*/

#define veccoord(VEC, AXIS)	(*(((float *) (VEC)) + AXIS))
#define xcoord(ATOM, AXIS)	veccoord(&(xfrag->atom[xatom[ATOM]].pos), AXIS)
#define ycoord(ATOM, AXIS)	veccoord(&(yfrag->atom[yatom[ATOM]].pos), AXIS)

enum {XAXIS, YAXIS, ZAXIS};

void
calcsuperposmatrix(int natoms, fragdetails *xfrag, int *xatom, fragdetails *yfrag, int *yatom, matrix *spos, vector *cgx, vector *cgy)
{
   matrix corr;
   int i, j, k, p, q, r, itno;
   float oneovern, sum, sigma, gamma, dist, qk, rk;
   BOOL rotated;

   if(natoms < 3)
   {
      printf("Not enough atoms (%d) to superposition molecules\n", natoms);
      return;
   }

   matrix_alloc(&corr, 3, 3, TRUE);

/* calc cgx and cgy */
   cgx->x = cgx->y = cgx->z = 0.0;
   cgy->x = cgy->y = cgy->z = 0.0;
   for(i=0; i<natoms; i++)
   {
      vector_add(cgx, &xfrag->atom[xatom[i]].pos, cgx);
      vector_add(cgy, &yfrag->atom[yatom[i]].pos, cgy);
   }
   oneovern = 1.0 / (float) natoms;
   vector_scale(cgx, oneovern);
   vector_scale(cgy, oneovern);

/* initialise the correlation matrix, c'' in the paper */
   for(j=XAXIS; j<=ZAXIS; j++)
   {
      for(k=XAXIS; k<=ZAXIS; k++)
      {
	 sum = 0.0;
	 for(i=0; i<natoms; i++)
	    sum += (xcoord(i, j) - veccoord(cgx, j)) * (ycoord(i, k) - veccoord(cgy, k));
	 matrix_setelement(&corr, k, j, sum);
      }
   }
/* spos is a bit easier to do, this is M in the paper */
   matrix_makeunit(spos);

   itno = 0;
   do
   {
   /* calc the rotation needed about each axis */
      rotated = FALSE;
      for(p=XAXIS; p<=ZAXIS; p++)
      {
	 q = (p+1 <= ZAXIS) ? p+1 : XAXIS;
	 r = (q+1 <= ZAXIS) ? q+1 : XAXIS;
      /* calc the rotation needed, angle = arctan(sigma/gamma) */
	 sigma = matrix_getelement(&corr, r, q) - matrix_getelement(&corr, q, r);
	 gamma = matrix_getelement(&corr, q, q) + matrix_getelement(&corr, r, r);
      /* if the angle is not too small update the two matrices */
	 if(fabs(sigma / gamma) > SUPERPOSTOLERANCE)
	 {
	    dist = (float) sqrt(gamma*gamma + sigma*sigma);
	    if(dist != 0.0)
	    {
	       for(k=XAXIS; k<=ZAXIS; k++)
	       {
	       /* the superposition matrix */
		  qk = matrix_getelement(spos, q, k);
		  rk = matrix_getelement(spos, r, k);
		  matrix_setelement(spos, q, k, (gamma*qk + sigma*rk) / dist);
		  matrix_setelement(spos, r, k, (-sigma*qk + gamma*rk) / dist);
	       /* and the correlation matrix */
		  qk = matrix_getelement(&corr, q, k);
		  rk = matrix_getelement(&corr, r, k);
		  matrix_setelement(&corr, q, k, (gamma*qk + sigma*rk) / dist);
		  matrix_setelement(&corr, r, k, (-sigma*qk + gamma*rk) / dist);
	       }
	    }
	    rotated = TRUE;			/* make a note of the fact that a rotation was still needed */
	 }
      }
      itno++;
   }
   while(rotated && itno < MAXSUPERPOSITS);	/* until we didnt rotate about any axis, or we give up */

   matrix_forget(&corr);

   return;
}

#undef veccoord
#undef xcoord
#undef ycoord

=================================================================
=================================================================

/*
   utils.c

   Simon Kilvington, 1994
*/

#include "utils.h"

float
vector_len2(vector *vp)
{
   return (vp->x * vp->x) + (vp->y * vp->y) + (vp->z * vp->z);
}

vector *
vector_add(vector *vp1, vector *vp2, vector *vp3)
{
   vp3->x = vp1->x + vp2->x;
   vp3->y = vp1->y + vp2->y;
   vp3->z = vp1->z + vp2->z;

   return vp3;
}

vector *
vector_sub(vector *vp1, vector *vp2, vector *vp3)
{
   vp3->x =vp1->x - vp2->x;
   vp3->y =vp1->y - vp2->y;
   vp3->z =vp1->z - vp2->z;

   return vp3;
}

vector *
vector_scale(vector *v, float scalar)
{
   v->x *= scalar;
   v->y *= scalar;
   v->z *= scalar;

   return v;
}

BOOL
matrix_alloc(matrix *m, int rows, int cols, BOOL quit)
{
   if(!(m->element = (matrix_element_t *) malloc(rows * cols * sizeof(matrix_element_t))))
   {
      printf("Out of memory in matrix_alloc\n");
      if(quit)
	 exit(1);
      return FALSE;
   }

   m->rows = rows;
   m->cols = cols;

   return TRUE;
}

matrix *
matrix_makeunit(matrix *m)
{
   int i, j;

   for(i=0; i<m->rows; i++)
      for(j=0; j<m->cols; j++)
	 matrix_setelement(m, i, j, (i==j) ? 1.0 : 0.0);

   return m;
}

matrix *
matrix_mult(matrix *m1, matrix *m2, matrix *r)
{
   int i, j, k;
   matrix_element_t e;

   if(m1->cols != m2->rows)
   {
      printf("Matrices are of incompatible sizes in matrix_mult\n");
      exit(1);
   }

   matrix_alloc(r, m1->rows, m2->cols, TRUE);

   for(i=0; i<m1->rows; i++)
   {
      for(j=0; j<m2->cols; j++)
      {
         e=0;
         for(k=0; k<m1->cols; k++)
            e += matrix_getelement(m1, i, k) * matrix_getelement(m2, k, j);
         matrix_setelement(r, i, j, e);
      }
   }

   return r;
}

=================================================================
=================================================================

/*
   utils.h

   Simon Kilvington, 1994
*/

#include <math.h>
#include <stdio.h>

typedef enum {FALSE, TRUE} BOOL;

typedef struct
{
   float x, y, z;
} vector;

float vector_len2(vector *);
#define vector_len(VP)		sqrt(vector_len2(VP))

vector *vector_add(vector *, vector *, vector *);
vector *vector_sub(vector *, vector *, vector *);

vector *vector_scale(vector *, float);

typedef float matrix_element_t;

/* a matrix is just a dope vector and a ptr to an array of elements */
typedef struct
{
   int rows, cols;
   matrix_element_t *element;
} matrix;

#define matrix_elementptr(mat, i, j)		(((mat)->element) + ((j) * (mat)->rows) + (i))
#define matrix_setelement(mat, i, j, val)	(*(matrix_elementptr(mat, i, j)) = (val))
#define matrix_getelement(mat, i, j)		(*(matrix_elementptr(mat, i, j)))

BOOL matrix_alloc(matrix *, int, int, BOOL);
#define matrix_forget(mat)			(free((mat)->element), (mat)->rows = (mat)->cols = 0)

matrix *matrix_makeunit(matrix *);

matrix *matrix_mult(matrix *, matrix *, matrix *);

=================================================================
=================================================================

Here endeth todays source.