method for finding axis of alpha helix from coordinates?

[Djamal Bouzida]

In response to <2511km$bnj at tamsun.tamu.edu>:

I've done some work on helices some time ago, and this is a C routine 
that I've written to find each axis corresponding to any number of 
helices (Nh which is read and returned from routine) given the C-alpha
coordinated.  The axis is found through the eigenvalues and eigenvectors 
of the least square matrix.

Please don't ask me any programming details. Looking for a job right now 
is in itself a job.  Just hire me ... (PhD in computational biophysics 
specializing in Monte Carlo simulations of biological molecules).

Djamal Bouzida.

--

int fnd_axes()
{
  void jacobi();

  int i,j,ih,Nh;
  int ivec;
  double lmd_max; /* largest eigenvalue */
  double t,sum,*r2;
  char cc;
  int ii;

  int nrot;
  double *lmbda,**V,**A;

  int count = 1;

  printf("Least Square Fit:\n");
  
  printf("\nHow many helices?\n");
  fscanf(in,"%d",&Nh);
  printf("%d\n",Nh);  

  /* Allocate Memory */
  helix = (struct helix_type *) malloc( (Nh+1) * sizeof (struct helix_type));
  if(!helix) nrerror("Can't allocate memory to structure helix");
  lmbda = dvector(1,3);
  V = dmatrix(1,3,1,3);
  A = dmatrix(1,3,1,3);

  for (ih=1; ih<=Nh; ih++) { /* Loop through all helices */

    printf(">> Helix %d:  ",ih);
    fscanf(in,"%d", &helix[ih].N);
    printf("%d points\n", helix[ih].N);

    helix[ih].xg = dvector(1,3);
    helix[ih].vg = dvector(1,3);
    helix[ih].d2 = dvector(1,helix[ih].N);
    helix[ih].x = dmatrix(1,helix[ih].N,1,3);
    helix[ih].xc = dmatrix(1,helix[ih].N,1,3);

    for (i=1; i<=helix[ih].N; i++){
      printf(">> Point %d:",i);
      for (ii=0; ii<20; ii++)
	fgetc(in,"%c", &cc);
      fscanf(in,"%lf %lf %lf", &helix[ih].x[i][1], &helix[ih].x[i][2], 
	&helix[ih].x[i][3]);

    }
    printf("\n");

    /* Find centroid */
    for (j=1; j<=3; j++) {
      helix[ih].xg[j] = 0.0;
      for (i=1; i<=helix[ih].N; i++) 
	helix[ih].xg[j] += helix[ih].x[i][j];
      helix[ih].xg[j] /= helix[ih].N;
      /* Coordinates in centroid system frame */
      for (i=1; i<=helix[ih].N; i++) 
	helix[ih].xc[i][j] = helix[ih].x[i][j] - helix[ih].xg[j];
    }

    /* Compute least square matrix */
    
    A[1][1] = 0.0;
    for (i=1; i<=helix[ih].N; i++)
      A[1][1] += helix[ih].xc[i][1]*helix[ih].xc[i][1];
    
    A[2][2] = 0.0;
    for (i=1; i<=helix[ih].N; i++)
      A[2][2] += helix[ih].xc[i][2]*helix[ih].xc[i][2];
    
    A[3][3] = 0.0;
    for (i=1; i<=helix[ih].N; i++)
      A[3][3] += helix[ih].xc[i][3]*helix[ih].xc[i][3];
    
    A[1][2] = 0.0;
    for (i=1; i<=helix[ih].N; i++)
      A[1][2] += helix[ih].xc[i][1]*helix[ih].xc[i][2];
    
    A[1][3] = 0.0;
    for (i=1; i<=helix[ih].N; i++)
      A[1][3] += helix[ih].xc[i][1]*helix[ih].xc[i][3];
    
    A[2][3] = 0.0;
    for (i=1; i<=helix[ih].N; i++)
      A[2][3] += helix[ih].xc[i][2]*helix[ih].xc[i][3];
    
    A[2][1]=A[1][2];
    A[3][1]=A[1][3];
    A[3][2]=A[2][3];
    
    printf("\nThe least square matrix is:\n\n");
    printf("    [ %8.3lf   %8.3lf   %8.3lf ]\n",A[1][1],A[1][2],A[1][3]);
    printf("    [ %8.3lf   %8.3lf   %8.3lf ]\n",A[2][1],A[2][2],A[2][3]);
    printf("    [ %8.3lf   %8.3lf   %8.3lf ]\n",A[3][1],A[3][2],A[3][3]);
    
    /* Compute eigenvectors and eigenvalues */
    jacobi(A,3,lmbda,V,&nrot);
    
    /* The eigenvalues are all positive! */
    lmd_max = lmbda[1];
    ivec = 1;
    if (lmbda[2] > lmd_max) {
      lmd_max = lmbda[2];
      ivec = 2;
    }
    if (lmbda[3] > lmd_max) {
      lmd_max = lmbda[3];
      ivec = 3;
    }
    
    /* Choose the direction of the eigenvector to be along the
       helix direction. To do that, the first point of the helix
       should correspond to a negative t. */
    t = 0.0;
    for (j=1; j<=3; j++) 
      t += helix[ih].xc[1][j]*V[j][ivec];
    if (t < 0.0) {
      helix[ih].vg[1] = V[1][ivec];
      helix[ih].vg[2] = V[2][ivec];
      helix[ih].vg[3] = V[3][ivec];
    }
    else {
      helix[ih].vg[1] = -V[1][ivec];
      helix[ih].vg[2] = -V[2][ivec];
      helix[ih].vg[3] = -V[3][ivec];
    }

    /* Store the t corresponding to the first point of the helix */
    helix[ih].ti = t;

    printf("\nLeast square line fit:\n");
    printf("            Direction:  %8.3lf %8.3lf %8.3lf\n",
	   helix[ih].vg[1],helix[ih].vg[2],helix[ih].vg[3]);
    printf("passing through point:  %8.3lf %8.3lf %8.3lf\n\n",
	   helix[ih].xg[1],helix[ih].xg[2],helix[ih].xg[3]);
    
    /* For each point i, find minimum distance^2 to line: d2[i] */
    printf(" #     t        d         r         l\n");
    r2 = dvector(1,helix[ih].N);
    sum = 0.0;
    for (i=1; i<=helix[ih].N; i++) {
      t = 0.0;
      for (j=1; j<=3; j++) 
	t += helix[ih].xc[i][j]*helix[ih].vg[j];
      printf(" %d %8.3lf ",i,t);
      helix[ih].d2[i] = 0.0;
      r2[i] = 0.0;
      for (j=1; j<=3; j++) {
	helix[ih].d2[i] += (helix[ih].xc[i][j] - helix[ih].vg[j]*t)*
	                   (helix[ih].xc[i][j] - helix[ih].vg[j]*t);
	r2[i] += helix[ih].xc[i][j]*helix[ih].xc[i][j];
      }
      sum += sqrt(helix[ih].d2[i]);
      printf("%8.3lf  ",sqrt(helix[ih].d2[i]));
      printf("%8.3lf  ",sqrt(r2[i]));
      printf("%8.3lf\n",sqrt(r2[i]-helix[ih].d2[i]));

/*      fprintf(out,"%d %8.3lf\n", count++, sqrt(helix[ih].d2[i]));  */
    }
/*    fprintf(out,"\n"); */

    /* Store the t corresponding to the last point of the helix */
    helix[ih].tf = t;

    printf("Sum of the distances: %8.3lf\n", sum);
    printf("Estimated Radius Rho: %8.3lf\n\n", sum/helix[ih].N);

  } /* Loop through all helices */

  free_dvector(lmbda,1,3);
  free_dmatrix(V,1,3,1,3);
  free_dmatrix(A,1,3,1,3);

  return(Nh);
}




--
--
Djamal Bouzida                            E-mail: djb at engcd.bu.edu