File: alpha_membrane_center.c

package info (click to toggle)
garlic 1.6-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 4,516 kB
  • sloc: ansic: 52,465; makefile: 2,254
file content (324 lines) | stat: -rw-r--r-- 9,562 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
/* Copyright (C) 2001-2003 Damir Zucic */

/*=============================================================================

			  alpha_membrane_center.c

Purpose:
	Refine the position of the membrane center. This version is suitable
	for helix bundle proteins.

Input:
	(1) Pointer to MolComplexS structure, with the chosen structure.
	(2) Pointer to VectorS structure,  with the vector which is perpend.
	    to the membrane.

Output:
	(1) The position of the membrane center  will be updated on success.
	(2) Return value.

Return value:
	(1) Positive on success.
	(2) Negative on failure.

Notes:
	(1) The method used here is described in alpha_membrane.c.  Read the
	    Note 1 and pay a special attention to the second step.

	(2) On every failure,  remember to free the allocated storage before
	    returning to the caller.

	(3) The CA is used as the representative atom because there are some
	    structures which contain only CA atoms.

	(4) Use calloc to initialize all array elements with zero values.

	(5) The cell width should be similar to the covalent bond length.

========includes:============================================================*/

#include <stdio.h>

#include <stdlib.h>

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <X11/Xatom.h>

#include "defines.h"
#include "typedefs.h"

/*======function prototypes:=================================================*/

double		AbsoluteValue_ (VectorS *);
int		ExtractCA_ (VectorS *, AtomS *, size_t, size_t);
double		ScalarProduct_ (VectorS *, VectorS *);

/*======refine the membrane center (helix bundle version):===================*/

int AlphaMembraneCenter_ (MolComplexS *mol_complexSP,
			  VectorS *normal_vectorSP)
{
int			residuesN, residueI;
size_t			double_size, int_size;
int			cellsN = 10000;
static double		*hydrophobicityP;
static int		*projectionsNP;
static int		*cells_usedNP;
double			abs_value, reciprocal_denominator;
VectorS			unit_vectorS;
VectorS			center_vectorS;
int			centralI;
double			cell_width = 1.0, half_cell_width;
double			reciprocal_cell_width;
ResidueS		*curr_residueSP;
int			first_atomI, last_atomI;
int			n;
AtomS			*first_atomSP;
double			hydrophobicity;
VectorS			CA_vectorS;
VectorS			radius_vectorS;
double			projection;
int			cellI;
int			used_cellsN;
double			average;
int			first_cellI, last_cellI;
int			half_window_width, windowI, combinedI;
double			max_average;
int			peak_cellI;
double			shift;

/* Copy the number of residues: */
residuesN = mol_complexSP->residuesN;

/* Double and int size: */
double_size = sizeof (double);
int_size = sizeof (int);

/* Allocate memory for the array of hydrophobicity values: */
hydrophobicityP = (double *) calloc (cellsN, double_size);
if (!hydrophobicityP)
	{
	return -1;
	}

/* Allocate the memory required to store the number of projections per cell: */
projectionsNP = (int *) calloc (cellsN, int_size);
if (!projectionsNP)
	{
	free (hydrophobicityP);
	return -2;
	}

/* Allocate the memory for the total number of used cells: */
cells_usedNP = (int *) calloc (cellsN, int_size);
if (!cells_usedNP)
	{
	free (hydrophobicityP);
	free (projectionsNP);
	return -3;
	}

/* Prepare the unit vector required to project CA atoms to the axis: */
abs_value = AbsoluteValue_ (normal_vectorSP);
if (abs_value == 0.0)
	{
	free (hydrophobicityP);
	free (projectionsNP);
	free (cells_usedNP);
	return -4;
	}
reciprocal_denominator = 1.0 / abs_value;
unit_vectorS.x = reciprocal_denominator * normal_vectorSP->x;
unit_vectorS.y = reciprocal_denominator * normal_vectorSP->y;
unit_vectorS.z = reciprocal_denominator * normal_vectorSP->z;

/* Prepare the vector which defines the membrane center: */
center_vectorS.x = mol_complexSP->membraneS.center_x;
center_vectorS.y = mol_complexSP->membraneS.center_y;
center_vectorS.z = mol_complexSP->membraneS.center_z;

/* Prepare  the central  index,  half of the */
/* cell width and the reciprocal cell width: */
centralI = cellsN / 2;
if (cell_width == 0.0)
	{
	free (hydrophobicityP);
	free (projectionsNP);
	free (cells_usedNP);
	return -4;
	}
half_cell_width = 0.5 * cell_width;
reciprocal_cell_width = 1.0 / cell_width;

/* Now project the hydrophobicities of all residues to the reference axis */
/* defined by the current membrane center  and by the normal vector.  The */
/* axis is divided into cells.  The width of a single cell is cell_width. */

/* Scan all residues and project CA atoms to the reference axis: */
for (residueI = 0; residueI < residuesN; residueI++)
	{
	/* Pointer to the current residue: */
	curr_residueSP = mol_complexSP->residueSP + residueI;

	/* Copy the indices of the first and of the last atom: */
	first_atomI = curr_residueSP->residue_startI;
	last_atomI  = curr_residueSP->residue_endI;

	/* Extract the position of CA atom. */
	n = ExtractCA_ (&CA_vectorS, mol_complexSP->atomSP,
			first_atomI, last_atomI);
	if (n != 1) continue;

	/* Copy the hydrophobicity associated with this residue: */
	first_atomSP = (mol_complexSP->atomSP + first_atomI); 
	hydrophobicity = first_atomSP->raw_atomS.hydrophobicity;

	/* Prepare the vector from the membrane center to the CA atom: */
	radius_vectorS.x = CA_vectorS.x - center_vectorS.x;
	radius_vectorS.y = CA_vectorS.y - center_vectorS.y;
	radius_vectorS.z = CA_vectorS.z - center_vectorS.z;

	/* Project this vector: */
	projection = ScalarProduct_ (&radius_vectorS, &unit_vectorS);

	/* Prepare and check the cell index: */
	cellI = centralI +
		reciprocal_cell_width * (projection - half_cell_width);
	if (cellI < 0) continue;
	if (cellI >= cellsN) continue;

	/* Update the total hydrophobicity stored */
	/* to the cell with the given cell index. */

	/* Update the number of CA atoms projected to the given cell: */
	(*(projectionsNP + cellI))++;

	/* Update the total hydrophobicity: */
	*(hydrophobicityP + cellI) += hydrophobicity;
	}

/* Find the range of cells which covers the entire macromolecular complex. */

/* Find the first cell which is not empty (search from bottom to top): */
first_cellI = 0;
for (cellI = 0; cellI < cellsN; cellI++)
	{
	if (*(projectionsNP + cellI) != 0)
		{
		first_cellI = cellI;
		break;
		}
	}

/* Find the last cell which is not empty (search from top to bottom): */
last_cellI = cellsN - 1;
for (cellI = cellsN - 1; cellI >= 0; cellI--)
	{
	if (*(projectionsNP + cellI) != 0)
		{
		last_cellI = cellI;
		break;
		}
	}

/* Prepare  the sliding window  which is required for  the averaging of */
/* total hydrophobicities. In this step the total hydrophobicities will */
/* be averaged over the specified number  of adjacent cells.  The width */
/* of the sliding window  should be  slightly larger  than the membrane */
/* thickness. The minimal half-width is 1 and the min. full width is 3. */

half_window_width = (int) (mol_complexSP->membraneS.thickness / cell_width);
half_window_width /= 2;
half_window_width += 1;

/* Average the total hydrophobicities over */
/* the specified number of adjacent cells: */
for (cellI = first_cellI; cellI <= last_cellI; cellI++)
	{
	/* Reset the number of cells in a window and the average value: */
	used_cellsN = 0;
	average = 0.0;

	/* Scan the sliding window: */
	for (windowI = -half_window_width;
	     windowI <= half_window_width; windowI++)
		{
		/* Prepare and check the combined index: */
		combinedI = cellI + windowI;
		if (combinedI < 0) continue;
		else if (combinedI >= cellsN) continue;

		/* Skip empty cells: */
		if (*(projectionsNP + combinedI) == 0) continue;

		/* Update the number of used cells and */
		/* the sum of  total hydrophobicities: */
		used_cellsN++;
		average += *(hydrophobicityP + combinedI);
		}

	/* Calculate the average value and store it. In addition, store the */
	/* number of cells which were used  to calculate the average value. */
	if (used_cellsN != 0)
		{
		*(hydrophobicityP + cellI) = average / (double) used_cellsN;
		*(cells_usedNP + cellI) = used_cellsN;
		}
	}

/* After averaging the total projected hydrophobicities,  find the */
/* index of the cell with the highest peak. In this loop, use only */
/* the average values  which were  calculated  from  a significant */
/* number of cells  (at least half of  the entire sliding window). */
/* Be sure to skip the cells  close to the edge of  the structure. */
max_average = -999999.0;
peak_cellI = -1;
for (cellI = first_cellI + half_window_width;
     cellI <= last_cellI - half_window_width; cellI++)
	{
	/* Skip  the cells for which  the average value */
	/* was calculated from a small number of cells: */
	if (*(cells_usedNP + cellI) < half_window_width) continue;

	/* Copy the average value stored to the current cell: */
	average = *(hydrophobicityP + cellI);

	/* Check is this value the maximum: */
	if (average > max_average)
		{
		peak_cellI = cellI;
		max_average = average;
		}
	}

/* Check the index of the cell with the highest peak: */
if (peak_cellI == -1)
	{
	free (hydrophobicityP);
	free (projectionsNP);
	free (cells_usedNP);
	return -5;
	}

/* Calculate the shift along the axis defined by the normal vector: */
shift = (double) (peak_cellI - centralI) * cell_width;

/* Calculate the new position of the membrane center: */
mol_complexSP->membraneS.center_x += shift * unit_vectorS.x;
mol_complexSP->membraneS.center_y += shift * unit_vectorS.y;
mol_complexSP->membraneS.center_z += shift * unit_vectorS.z;

/* Free storage: */
free (hydrophobicityP);
free (projectionsNP);
free (cells_usedNP);

/* Return positive value on success: */
return 1;
}

/*===========================================================================*/