1 /* @source pepwindow application **********************************************
2 **
3 ** Displays protein hydropathy.
4 ** @author Copyright (C) Ian Longden (il@sanger.ac.uk)
5 ** @@
6 ** Original program by Jack Kyte and Russell F. Doolittle.
7 **
8 ** This program is free software; you can redistribute it and/or
9 ** modify it under the terms of the GNU General Public License
10 ** as published by the Free Software Foundation; either version 2
11 ** of the License, or (at your option) any later version.
12 **
13 ** This program is distributed in the hope that it will be useful,
14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 ** GNU General Public License for more details.
17 **
18 ** You should have received a copy of the GNU General Public License
19 ** along with this program; if not, write to the Free Software
20 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
21 ******************************************************************************/
22
23 #include <limits.h>
24 #include <float.h>
25 #include "emboss.h"
26
27 #define AZ 28
28
29
30
31
32 static AjBool pepwindow_getnakaidata(AjPFile file, float matrix[],
33 AjBool normal);
34
35
36
37
38 /* @prog pepwindow ************************************************************
39 **
40 ** Displays protein hydropathy
41 **
42 ******************************************************************************/
43
main(int argc,char ** argv)44 int main(int argc, char **argv)
45 {
46 AjPFile datafile = NULL;
47 AjPStr aa0str = NULL;
48 const char *s1 = NULL;
49 AjPSeq seq = NULL;
50 ajuint llen = 0U;
51 AjBool normal = AJFALSE;
52 float matrix[AZ];
53 ajuint i = 0U;
54 ajuint midpoint = 0U;
55 ajuint j = 0U;
56 ajuint k = 0U;
57 AjPGraphdata graphdata = NULL;
58 AjPGraph mult = NULL;
59 float totalmin = FLT_MAX;
60 float totalmax = -FLT_MAX;
61 float total = 0.0F;
62 float fstart = 0.0F;
63 float fend = 0.0F;
64 ajuint ilen = 0U;
65 ajuint istart = 0U;
66 ajuint iend = 0U;
67
68 embInit("pepwindow", argc, argv);
69
70 seq = ajAcdGetSeq("sequence");
71
72 mult = ajAcdGetGraphxy("graph");
73 datafile = ajAcdGetDatafile("datafile");
74 llen = ajAcdGetInt("window");
75 normal = ajAcdGetBoolean("normalize");
76
77 s1 = ajStrGetPtr(ajSeqGetSeqS(seq));
78
79 midpoint = (ajint)((llen+1)/2);
80
81 istart = ajSeqGetBegin(seq) - 1;
82 iend = ajSeqGetEnd(seq);
83
84 aa0str = ajStrNewRes(iend+1);
85
86 if((iend-istart) > llen)
87 ilen = iend-istart+1-llen;
88 else
89 ilen = 1;
90
91 graphdata = ajGraphdataNewI(ilen);
92
93 ajGraphdataSetTypeC(graphdata,"2D Plot");
94
95 ajGraphDataAdd(mult,graphdata);
96
97 for(i = 0U; i < iend; i++)
98 ajStrAppendK(&aa0str,(char)ajBasecodeToInt(*s1++));
99
100
101 if(!pepwindow_getnakaidata(datafile,&matrix[0], normal))
102 ajExitBad();
103
104 s1 = ajStrGetPtr(aa0str) + istart;
105
106 k = 0U;
107 for(i = istart; iend >= llen && i <= (iend - llen); i++)
108 {
109 total = 0.0F;
110 for(j=0;j<llen;j++)
111 total += matrix[(ajint)s1[j]];
112
113 total /= (float)llen;
114 graphdata->x[k] = (float) (i + midpoint);
115 graphdata->y[k] = total;
116
117 totalmax = AJMAX(total, totalmax);
118 totalmin = AJMIN(total, totalmin);
119
120 k++;
121 s1++;
122 }
123 fstart = (float) istart;
124 fend = (float) iend;
125
126 ajGraphdataSetTruescale(graphdata,fstart,fend,totalmin,totalmax);
127
128 if(totalmin <0)
129 totalmin = totalmin*(float)1.1;
130 else
131 totalmin = totalmin/(float)1.1;
132
133 if(totalmax > 0)
134 totalmax = totalmax*(float)1.1;
135 else
136 totalmax = totalmax/(float)1.1;
137
138 ajGraphdataSetMinmax(graphdata,fstart,fend,totalmin,totalmax);
139 ajGraphxySetMinmax(mult,fstart,fend,totalmin,totalmax);
140
141 ajGraphxyDisplay(mult,AJTRUE);
142 ajGraphxyDel(&mult);
143
144 ajFileClose(&datafile);
145 ajSeqDel(&seq);
146
147 ajStrDel(&aa0str);
148
149 embExit();
150
151 return EXIT_SUCCESS;
152 }
153
154
155
156
157 /* @funcstatic pepwindow_getnakaidata *****************************************
158 **
159 ** Read the NAKAI (AAINDEX) data file
160 **
161 ** @param [u] file [AjPFile] Input file
162 ** @param [w] matrix [float[]] Data values for each amino acid
163 ** @param [r] normal [AjBool] If true, normalize data to mean 0.0 and SD 1.0
164 ** @return [AjBool] ajTrue on success
165 ** @@
166 ******************************************************************************/
167
168
pepwindow_getnakaidata(AjPFile file,float matrix[],AjBool normal)169 static AjBool pepwindow_getnakaidata(AjPFile file, float matrix[],
170 AjBool normal)
171 {
172 AjPStr buffer = NULL;
173 AjPStr buf2 = NULL;
174 AjPStr delim = NULL;
175 AjPStr description = NULL;
176 AjPStrTok token = NULL;
177 ajint line = 0;
178 const char *ptr = NULL;
179 ajint i = 0;
180
181 if(!file)
182 return ajFalse;
183
184 delim = ajStrNewC(" :\t\n");
185 buffer = ajStrNew();
186 buf2 = ajStrNew();
187 description = ajStrNew();
188
189 for (i = 0; i < 26; i++) {
190 matrix[i] = FLT_MIN;
191 }
192
193 while(ajReadline(file,&buffer))
194 {
195 ptr = ajStrGetPtr(buffer);
196 if(*ptr == 'D') /* save description */
197 ajStrAssignS(&description, buffer);
198 else if(*ptr == 'I')
199 line = 1;
200 else if(line == 1)
201 {
202 line++;
203 ajStrRemoveWhiteExcess(&buffer);
204
205 token = ajStrTokenNewS(buffer,delim);
206
207 ajStrTokenNextParseS(token,delim,&buf2);
208 ajStrToFloat(buf2,&matrix[0]);
209
210 ajStrTokenNextParseS(token,delim,&buf2);
211 ajStrToFloat(buf2,&matrix[17]);
212
213 ajStrTokenNextParseS(token,delim,&buf2);
214 ajStrToFloat(buf2,&matrix[13]);
215
216 ajStrTokenNextParseS(token,delim,&buf2);
217 ajStrToFloat(buf2,&matrix[3]);
218
219 ajStrTokenNextParseS(token,delim,&buf2);
220 ajStrToFloat(buf2,&matrix[2]);
221
222 ajStrTokenNextParseS(token,delim,&buf2);
223 ajStrToFloat(buf2,&matrix[16]);
224
225 ajStrTokenNextParseS(token,delim,&buf2);
226 ajStrToFloat(buf2,&matrix[4]);
227
228 ajStrTokenNextParseS(token,delim,&buf2);
229 ajStrToFloat(buf2,&matrix[6]);
230
231 ajStrTokenNextParseS(token,delim,&buf2);
232 ajStrToFloat(buf2,&matrix[7]);
233
234 ajStrTokenNextParseS(token,delim,&buf2);
235 ajStrToFloat(buf2,&matrix[8]);
236
237 ajStrTokenDel(&token);
238 }
239 else if(line == 2)
240 {
241 line++;
242
243 ajStrRemoveWhiteExcess(&buffer);
244 token = ajStrTokenNewS(buffer,delim);
245
246 ajStrTokenNextParseS(token,delim,&buf2);
247 ajStrToFloat(buf2,&matrix[11]);
248
249 ajStrTokenNextParseS(token,delim,&buf2);
250 ajStrToFloat(buf2,&matrix[10]);
251
252 ajStrTokenNextParseS(token,delim,&buf2);
253 ajStrToFloat(buf2,&matrix[12]);
254
255 ajStrTokenNextParseS(token,delim,&buf2);
256 ajStrToFloat(buf2,&matrix[5]);
257
258 ajStrTokenNextParseS(token,delim,&buf2);
259 ajStrToFloat(buf2,&matrix[15]);
260
261 ajStrTokenNextParseS(token,delim,&buf2);
262 ajStrToFloat(buf2,&matrix[18]);
263
264 ajStrTokenNextParseS(token,delim,&buf2);
265 ajStrToFloat(buf2,&matrix[19]);
266
267 ajStrTokenNextParseS(token,delim,&buf2);
268 ajStrToFloat(buf2,&matrix[22]);
269
270 ajStrTokenNextParseS(token,delim,&buf2);
271 ajStrToFloat(buf2,&matrix[24]);
272
273 ajStrTokenNextParseS(token,delim,&buf2);
274 ajStrToFloat(buf2,&matrix[21]);
275
276 ajStrTokenDel(&token);
277 }
278 }
279
280 if(normal)
281 embPropNormalF(matrix, FLT_MIN);
282
283 embPropFixF(matrix, FLT_MIN);
284
285 ajStrDel(&buffer);
286 ajStrDel(&description);
287 ajStrDel(&buf2);
288 ajStrDel(&delim);
289
290 return ajTrue;
291 }
292