Biosequences .. Software .. Molbio soft .. Network News .. FTP

# method for finding axis of alpha helix from coordinates?

Djamal Bouzida djb at engcd.bu.edu
Fri Sep 3 14:25:48 EST 1993

```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);

} /* 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

```