1 #include <objres.h>
2 #include <objseq.h>
3 #include <seqport.h>
4 #include <aacomp.h>
5
6 /* defines */
7 #define F7 7
8 #define F24 24
9 #define F42 42
10 #define PI 3.14159
11
12 /* functions - internal */
13
AAC(SeqPortPtr spp,Int4 start,Int4 end,Int4 numval,Int4 window,FloatHiPtr scr,CharPtr res)14 static FloatHiPtr AAC (SeqPortPtr spp, Int4 start, Int4 end,
15 Int4 numval, Int4 window,
16 FloatHiPtr scr, CharPtr res)
17 {
18 Int4 i, n, col;
19 Int4 stop;
20 Char aa;
21 Boolean flagMatch;
22 FloatHi temp, fmin = 0.0, fmax = 0.0, frange;
23 FloatHiPtr frealheadptr, fheadptr;
24 FloatHiPtr fptr;
25
26 fptr = (FloatHiPtr) MemNew ((sizeof (FloatHi)) * numval);
27 frealheadptr = fptr;
28 for (i = start; i < end+1; i++)
29 *fptr++ = 0.0;
30 fptr = frealheadptr;
31
32 for (i = 0; i < window/2; i++)
33 fptr++;
34 fheadptr = fptr;
35 stop = end + 1 - window + 1;
36 for (i = start; i < stop; i++)
37 {
38 temp = 0;
39 SeqPortSeek (spp, i, SEEK_SET);
40 for (col = 0; col < window; col++)
41 {
42 aa = SeqPortGetResidue (spp);
43 flagMatch = FALSE;
44 n = 0;
45 while (res[n] != '\0')
46 {
47 if (aa == res[n])
48 {
49 flagMatch = TRUE;
50 break;
51 }
52 n++;
53 }
54 if (flagMatch)
55 temp += scr[n];
56 }
57 temp /= window;
58 *fptr = temp;
59 fptr++;
60 }
61 fptr = fheadptr;
62 for (i = start; i < stop; i++)
63 {
64 if (*fptr > fmax)
65 fmax = *fptr;
66 fptr++;
67 }
68 fptr = fheadptr;
69 for (i = start; i < stop; i++)
70 {
71 if (*fptr < fmin)
72 fmin = *fptr;
73 fptr++;
74 }
75 frange = fmax-fmin;
76 fptr = fheadptr;
77 for (i = start; i < stop; i++)
78 {
79 *fptr = (*fptr - fmin)/frange*100;
80 fptr++;
81 }
82
83 return frealheadptr;
84 }
85
AACGraph(SeqPortPtr spp,Int4 start,Int4 end,Int4 window,FloatHiPtr scr,CharPtr res)86 static SeqGraphPtr AACGraph (SeqPortPtr spp, Int4 start, Int4 end, Int4 window,
87 FloatHiPtr scr, CharPtr res)
88 {
89 SeqGraphPtr sgp;
90 Int4 gwidth = 500;
91
92 if ((sgp = SeqGraphNew ()) != NULL)
93 {
94 /* type and number of values and compression */
95 sgp->numval = end + 1;
96 sgp->compr = (Int4) (sgp->numval / gwidth);
97 if ((sgp->numval%gwidth) != 0)
98 sgp->compr += 1;
99 /* graph type */
100 sgp->flags[2] = 1;
101 sgp->values = (Pointer) MemNew ((sizeof (FloatHi)) * sgp->numval);
102 /* min/max */
103 sgp->max.realvalue = 100.0;
104 sgp->min.realvalue = 0.0;
105 sgp->axis.realvalue = 0.0;
106 /* scaling */
107 sgp->flags[1] = 1;
108 sgp->a = 2;
109 sgp->b = 0;
110 /* do it */
111 AAC (spp, start, end, end+1, window, scr, res);
112 }
113 return sgp;
114 }
115
116 /*
117 this stores a score for each of 24 resides
118 20 amino acids
119 2 ambiguous amino acids (ASX GLX)
120 1 unknown amino acid (the average score of the 20 generally)
121 1 stop codon (currently also the average score of the 20)
122 the scoring table should have the scores in a column
123 in order of the residues stored in static Char c[]
124
125 the caller must pass pointers to arrays that can hold
126 the 24 residues (res) and scores (scr)
127
128 a successful, tho not necessarily correct, read will return 24
129 */
130
ReadAAC(CharPtr filename,FloatHiPtr scr,CharPtr res)131 extern Int4 ReadAAC (CharPtr filename, FloatHiPtr scr, CharPtr res)
132 {
133 Int2 val, count;
134 Int4 type;
135 long int ltype;
136 FloatHi score;
137 Char buf[PATH_MAX], buff[PATH_MAX];
138 CharPtr bptr;
139 FILE *fio;
140
141 static Char c[] = {'A','B','C','D','E','F','G','H','I','K','L','M',
142 'N','P','Q','R','S','T','V','W','X','Y','Z','*'};
143
144 if (!FindPath ("ncbi", "ncbi", "data", buf, sizeof (buf)))
145 return 0;
146 FileBuildPath (buf, NULL, filename);
147 if ((fio = FileOpen (buf, "r")) == NULL)
148 {
149 if ((fio = FileOpen (filename, "r")) == NULL)
150 {
151 return 0;
152 }
153 }
154
155 val = 0;
156 type = -1;
157 while ((FileGets (buff, sizeof (buff), fio)) != NULL)
158 {
159 if (buff[0] == ';')
160 continue;
161 if (type == -1)
162 {
163 sscanf (buff, "%ld", <ype);
164 type = (Int4) ltype;
165 if (type != 1)
166 {
167 FileClose (fio);
168 return 0;
169 }
170 FileGets (buff, sizeof (buff), fio);
171 continue;
172 }
173
174 res[val] = c[val];
175 bptr = buff;
176
177 sscanf (bptr, "%lf", &score);
178 scr[val] = score;
179
180 val++;
181 if (val == F24)
182 break;
183 }
184 count = val;
185
186 FileClose (fio);
187 return (Int4) count;
188 }
189
190 /* functions - external */
191
AAComposition(SeqPortPtr spp,Int4 start,Int4 end,Int4 numval,FloatHiPtr scr,CharPtr res,Int4 window)192 extern FloatHiPtr AAComposition (SeqPortPtr spp,
193 Int4 start, Int4 end, Int4 numval,
194 FloatHiPtr scr, CharPtr res, Int4 window)
195 {
196 FloatHiPtr fptr;
197
198 fptr = AAC (spp, start, end, numval, window, scr, res);
199 return fptr;
200 }
201
202