1 /* @source findkm.c
2 **
3 ** @author Copyright (C) Sinead O'Leary (soleary@hgmp.mrc.ac.uk),
4 ** David Martin (david.martin@biotek.uio.no)
5 **
6 ** Application to calculate the Michaelis Menton Constants (Km) of different
7 ** enzymes and their substrates.
8 **
9 ** This program is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU General Public License
11 ** as published by the Free Software Foundation; either version 2
12 ** of the License, or (at your option) any later version.
13 **
14 ** This program is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 ** GNU General Public License for more details.
18 **
19 ** You should have received a copy of the GNU General Public License
20 ** along with this program; if not, write to the Free Software
21 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22 ******************************************************************************/
23
24 #include "emboss.h"
25 #include <math.h>
26 #include <stdlib.h>
27 #include <limits.h>
28
29
30
31
32 static float findkm_summation(const float *arr, ajint number);
33 static float findkm_multisum (const float *arr1, const float *arr2,
34 ajint number);
35 static float findkm_findmax(const float *arr1, ajint number);
36 static float findkm_findmin(const float *arr1, ajint number);
37
38
39
40
41 /* @prog findkm ***************************************************************
42 **
43 ** Find Km and Vmax for an enzyme reaction by a Hanes/Woolf plot
44 **
45 ******************************************************************************/
46
main(int argc,char ** argv)47 int main(int argc, char **argv)
48 {
49 AjPFile infile = NULL;
50 AjPFile outfile = NULL;
51 AjPStr line;
52 AjPGraph graphLB = NULL;
53 AjPGraphdata xygraph = NULL;
54 AjPGraphdata xygraph2 = NULL;
55 AjBool doplot;
56
57 ajint N=0;
58
59 float *xdata = NULL;
60 float *ydata = NULL;
61 float *V = NULL;
62 float *S = NULL;
63
64 float a;
65 float b;
66 float upperXlimit;
67 float upperYlimit;
68
69 float A;
70 float B;
71 float C;
72 float D;
73 float xmin;
74 float xmax;
75 float ymin;
76 float ymax;
77 float xmin2;
78 float xmax2;
79 float ymin2;
80 float ymax2;
81
82 float Vmax;
83 float Km;
84 float cutx;
85 float cuty;
86
87 float amin = 0.;
88 float amax = 0.;
89 float bmin = 0.;
90 float bmax = 0.;
91
92
93 embInit("findkm", argc, argv);
94
95 infile = ajAcdGetInfile("infile");
96 outfile = ajAcdGetOutfile ("outfile");
97 doplot = ajAcdGetBoolean("plot");
98 graphLB = ajAcdGetGraphxy("graphLB");
99 line = ajStrNew();
100
101
102 /* Determine N by reading infile */
103
104 while(ajReadlineTrim(infile, &line))
105 if(ajStrGetLen(line) >0)
106 N++;
107
108
109 /* only allocate memory to the arrays */
110
111 AJCNEW(xdata, N);
112 AJCNEW(ydata, N);
113 AJCNEW(S, N);
114 AJCNEW(V, N);
115
116 ajFileSeek(infile, 0L, 0);
117
118 N=0;
119 while(ajReadlineTrim(infile, &line))
120 {
121 if(ajStrGetLen(line) > 0)
122 {
123 sscanf(ajStrGetPtr(line),"%f %f",&S[N],&V[N]);
124 if(S[N] > 0.0 && V[N] > 0.0)
125 {
126 xdata[N] = S[N];
127 ydata[N] = S[N]/V[N];
128 N++;
129 }
130 }
131 }
132
133
134 /* find the max and min values for the graph parameters*/
135 xmin = (float)0.5*findkm_findmin(xdata, N);
136 xmax = (float)1.5*findkm_findmax(xdata, N);
137 ymin = (float)0.5*findkm_findmin(ydata, N);
138 ymax = (float)1.5*findkm_findmax(ydata, N);
139
140 xmin2 = (float)0.5*findkm_findmin(S, N);
141 xmax2 = (float)1.5*findkm_findmax(S, N);
142 ymin2 = (float)0.5*findkm_findmin(V, N);
143 ymax2 = (float)1.5*findkm_findmax(V, N);
144
145
146
147 /*
148 ** In case the casted ints turn out to be same number on the axis,
149 ** make the max number larger than the min so graph can be seen.
150 */
151
152 if((ajint)xmax == (ajint)xmin)
153 ++xmax;
154 if((ajint)ymax == (ajint)ymin)
155 ++ymax;
156
157
158 if((ajint)xmax2 == (ajint)xmin2)
159 ++xmax2;
160 if((ajint)ymax2 == (ajint)ymin2)
161 ++ymax2;
162
163
164
165 /*
166 ** Gaussian Elimination for Best-fit curve plotting and
167 ** calculating Km and Vmax
168 */
169
170 A = findkm_summation(xdata, N);
171 B = findkm_summation(ydata, N);
172
173 C = findkm_multisum(xdata, ydata, N);
174 D = findkm_multisum(xdata, xdata, N);
175
176
177 /*
178 ** To find the best fit line, Least Squares Fit: y =ax +b;
179 ** Two Simultaneous equations, REARRANGE FOR b
180 **
181 ** findkm_summation(ydata, N) - findkm_summation(xdata,N)*a - N*b =0;
182 ** b = (findkm_summation(ydata,N) - findkm_summation(xdata,N)*a) / N;
183 ** b = (B - A*a)/ N;
184 **
185 ** C - D*a - A*((B - A*a)/ N) =0;
186 ** C - D*a - A*B/N + A*A*a/N =0;
187 ** C - A*B/N = D*a - A*A*a/N;
188 */
189
190 /* REARRANGE FOR a */
191
192 a = (N*C - A*B)/ (N*D - A*A);
193 b = (B - A*a)/ N;
194
195 /*
196 ** Equation of Line - Lineweaver burk eqn
197 ** 1/V = (Km/Vmax)*(1/S) + 1/Vmax;
198 */
199
200
201 Vmax = 1/a;
202 Km = b/a;
203
204 cutx = -1/Km;
205 cuty = Km/Vmax;
206
207 /* set limits for last point on graph */
208
209 upperXlimit = findkm_findmax(xdata,N)+3;
210 upperYlimit = (upperXlimit)*a + b;
211
212 ajFmtPrintF(outfile, "---Hanes Woolf Plot Calculations---\n");
213 ajFmtPrintF(outfile, "Slope of best fit line is a = %.2f\n", a);
214 ajFmtPrintF(outfile,"Coefficient in Eqn of line y = ma +b is b "
215 "= %.2f\n", b);
216
217 ajFmtPrintF(outfile, "Where line cuts x axis = (%.2f, 0)\n", cutx);
218 ajFmtPrintF(outfile, "Where line cuts y axis = (0, %.2f)\n", cuty);
219 ajFmtPrintF(outfile, "Limit-point of graph for plot = (%.2f, %.2f)\n\n",
220 upperXlimit, upperYlimit);
221 ajFmtPrintF(outfile, "Vmax = %.2f, Km = %f\n",Vmax, Km);
222
223 /* draw graphs */
224
225 if(doplot)
226 {
227 xygraph = ajGraphdataNewI(N);
228 ajGraphdataAddXY(xygraph, S, V);
229 ajGraphDataAdd(graphLB, xygraph);
230 ajGraphdataSetTitleC(xygraph, "Michaelis Menten Plot");
231 ajGraphdataSetXlabelC(xygraph, "[S]");
232 ajGraphdataSetYlabelC(xygraph, "V");
233
234 ajGraphxySetXstartF(graphLB, 0.0);
235 ajGraphxySetXendF(graphLB, xmax2);
236 ajGraphxySetYstartF(graphLB, 0.0);
237 ajGraphxySetYendF(graphLB, ymax2);
238 ajGraphxySetXrangeII(graphLB, (ajint)0.0, (ajint)xmax2);
239 ajGraphxySetYrangeII(graphLB, (ajint)0.0, (ajint)ymax2);
240 ajGraphdataAddposLine(xygraph, 0.0, 0.0, S[0], V[0], (ajint)BLACK);
241 ajGraphxyShowPointsCircle(graphLB, ajTrue);
242 ajGraphdataSetMinmax(xygraph,0.0,xmax2,0.0,ymax2);
243
244
245 ajGraphicsCalcRange(S,N,&amin,&amax);
246 ajGraphicsCalcRange(V,N,&bmin,&bmax);
247 ajGraphdataSetTruescale(xygraph,amin,amax,bmin,bmax);
248 ajGraphdataSetTypeC(xygraph,"2D Plot Float");
249
250 xygraph2 = ajGraphdataNewI(N);
251 ajGraphdataAddXY(xygraph2, xdata, ydata);
252 ajGraphDataAdd(graphLB, xygraph2);
253
254 ajGraphdataSetTitleC(xygraph2, "Hanes Woolf Plot");
255 ajGraphdataSetXlabelC(xygraph2, "[S]");
256 ajGraphdataSetYlabelC(xygraph2, "[S]/V");
257
258 ajGraphxySetXstartF(graphLB, cutx);
259 ajGraphxySetXendF(graphLB, upperXlimit);
260 ajGraphxySetYstartF(graphLB, 0.0);
261 ajGraphxySetYendF(graphLB, upperYlimit);
262 ajGraphxySetXrangeII(graphLB, (ajint)cutx, (ajint)upperXlimit);
263 ajGraphxySetYrangeII(graphLB, (ajint)0.0, (ajint)upperYlimit);
264
265 ajGraphxyShowPointsCircle(graphLB, ajTrue);
266 ajGraphdataSetMinmax(xygraph2, cutx,upperXlimit,0.0,upperYlimit);
267 ajGraphicsCalcRange(xdata,N,&amin,&amax);
268 ajGraphicsCalcRange(ydata,N,&bmin,&bmax);
269 ajGraphdataSetTruescale(xygraph2,amin,amax,bmin,bmax);
270 ajGraphdataSetTypeC(xygraph2,"2D Plot");
271
272
273
274 ajGraphSetTitleC(graphLB,"FindKm");
275 ajGraphxySetflagOverlay(graphLB,ajFalse);
276 ajGraphxyDisplay(graphLB, ajTrue);
277 }
278
279 AJFREE(xdata);
280 AJFREE(ydata);
281
282 AJFREE(S);
283 AJFREE(V);
284
285 ajFileClose(&infile);
286 ajFileClose(&outfile);
287
288 ajGraphxyDel(&graphLB);
289 ajStrDel(&line);
290
291 embExit();
292
293 return 0;
294 }
295
296
297
298
299 /* @funcstatic findkm_summation ***********************************************
300 **
301 ** Undocumented.
302 **
303 ** @param [r] arr [const float*] Undocumented
304 ** @param [r] number [ajint] Array size
305 ** @return [float] Undocumented
306 ** @@
307 ******************************************************************************/
308
309
findkm_summation(const float * arr,ajint number)310 static float findkm_summation(const float *arr, ajint number)
311 {
312 ajint i;
313 float sum=0;
314
315 for(i = 0; i < number; i++)
316 sum += arr[i];
317
318 return sum;
319 }
320
321
322
323
324 /* @funcstatic findkm_multisum ************************************************
325 **
326 ** Undocumented.
327 **
328 ** @param [r] arr1 [const float*] Undocumented
329 ** @param [r] arr2 [const float*] Undocumented
330 ** @param [r] number [ajint] Array size
331 ** @return [float] Undocumented
332 ** @@
333 ******************************************************************************/
334
findkm_multisum(const float * arr1,const float * arr2,ajint number)335 static float findkm_multisum(const float *arr1, const float *arr2,
336 ajint number)
337 {
338 ajint i;
339 float sum = 0;
340
341 for(i = 0; i < number; i++)
342 sum += arr1[i]*arr2[i];
343
344 return sum;
345 }
346
347
348
349
350 /* @funcstatic findkm_findmax *************************************************
351 **
352 ** Undocumented.
353 **
354 ** @param [r] arr [const float*] Undocumented
355 ** @param [r] number [ajint] Array size
356 ** @return [float] Undocumented
357 ** @@
358 ******************************************************************************/
359
findkm_findmax(const float * arr,ajint number)360 static float findkm_findmax(const float *arr, ajint number)
361 {
362 ajint i;
363 float max;
364
365 max = arr[0];
366
367 for(i=1; i<number; ++i)
368 if(arr[i] > max)
369 max = arr[i];
370
371 return max;
372 }
373
374
375
376
377 /* @funcstatic findkm_findmin *************************************************
378 **
379 ** Undocumented.
380 **
381 ** @param [r] arr [const float*] Undocumented
382 ** @param [r] number [ajint] Array size
383 ** @return [float] Undocumented
384 ** @@
385 ******************************************************************************/
386
findkm_findmin(const float * arr,ajint number)387 static float findkm_findmin(const float *arr, ajint number)
388 {
389 ajint i;
390 float min;
391
392 min = arr[0];
393
394 for(i=1; i<number; ++i)
395 if(arr[i] < min)
396 min = arr[i];
397
398 return min;
399 }
400