1 #include <objres.h>
2 #include <objseq.h>
3 #include <seqport.h>
4 #include <seqpcc.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
PPCC(FloatHi score)14 static FloatHi PPCC (FloatHi score)
15 {
16 FloatHi meang = 0.77;
17 FloatHi stdg = 0.20;
18 FloatHi meanc = 1.63;
19 FloatHi stdc = 0.24;
20 FloatHi constg, constc;
21 FloatHi evalg, evalc;
22 FloatHi Gg, Gc;
23 FloatHi prob;
24
25 constg = 1 / (stdg * 2 * PI);
26 constc = 1 / (stdc * 2 * PI);
27 evalg = (((score - meang) / stdg) * ((score - meang) / stdg)) / 2;
28 evalc = (((score - meanc) / stdc) * ((score - meanc) / stdc)) / 2;
29 evalg = exp (evalg);
30 evalc = exp (evalc);
31 Gg = constg / evalg;
32 Gc = constc / evalc;
33 prob = Gc / ((30 * Gg) + Gc);
34 return prob;
35 }
36
nroot(FloatHi fv,Int4 wv)37 static FloatHi nroot (FloatHi fv, Int4 wv)
38 {
39 FloatHi fi = 1.0;
40
41 fi = fi / (FloatHi) wv;
42 fv = pow (fv, fi);
43 return fv;
44 }
45
PCC(SeqPortPtr spp,Int4 start,Int4 end,Int4 window,FloatHi scr[F24][F42],Char res[F24],FloatHiPtr fptr)46 static void PCC (SeqPortPtr spp, Int4 start, Int4 end, Int4 window,
47 FloatHi scr[F24][F42], Char res[F24], FloatHiPtr fptr)
48 {
49 Int4 i, n, col;
50 Int4 stop;
51 Char aa;
52 Boolean flagMatch;
53 FloatHi temp;
54 FloatHiPtr fheadptr;
55
56 fheadptr = fptr;
57 for (i = start; i < end; i++)
58 *fptr++ = 0.0;
59 fptr = fheadptr;
60 stop = end - window;
61 for (i = start; i < stop; i++)
62 {
63 SeqPortSeek (spp, i, SEEK_SET);
64 temp = 1.0;
65 for (col = 0; col < window; col++)
66 {
67 aa = SeqPortGetResidue (spp);
68 flagMatch = FALSE;
69 n = 0;
70 while (res[n] != '\0')
71 {
72 if (aa == res[n])
73 {
74 flagMatch = TRUE;
75 break;
76 }
77 n++;
78 }
79 if (flagMatch)
80 temp *= scr[n][col];
81 }
82 temp = nroot (temp, window);
83 temp = PPCC (temp);
84 if (temp < 0.0)
85 temp = 0.0;
86 temp *= 100.0;
87 fheadptr = fptr;
88 for (n = 0; n < window; n++)
89 {
90 if (temp > *fptr)
91 *fptr = temp;
92 fptr++;
93 }
94 fptr = fheadptr;
95 fptr++;
96 }
97 return;
98 }
99
PCCGraph(SeqPortPtr spp,Int4 start,Int4 end,Int4 window,FloatHi scr[F24][F42],Char res[F24])100 static SeqGraphPtr PCCGraph (SeqPortPtr spp, Int4 start, Int4 end, Int4 window,
101 FloatHi scr[F24][F42], Char res[F24])
102 {
103 SeqGraphPtr sgp;
104 Int4 gwidth = 500;
105 FloatHiPtr fptr;
106
107 if ((sgp = SeqGraphNew ()) != NULL)
108 {
109 /* type and number of values and compression */
110 sgp->numval = end + 1;
111 sgp->compr = (Int4) (sgp->numval / gwidth);
112 if ((sgp->numval%gwidth) != 0)
113 sgp->compr += 1;
114 /* graph type */
115 sgp->flags[2] = 1;
116 sgp->values = (Pointer) MemNew ((sizeof (FloatHi)) * sgp->numval);
117 /* min/max */
118 sgp->max.realvalue = 100.0;
119 sgp->min.realvalue = 0.0;
120 sgp->axis.realvalue = 0.0;
121 /* scaling */
122 sgp->flags[1] = 1;
123 sgp->a = 2;
124 sgp->b = 0;
125 /* do it */
126 fptr = (FloatHiPtr) sgp->values;
127 PCC (spp, start, end, window, scr, res, fptr);
128 }
129 return sgp;
130 }
131
ReadPCC(CharPtr res,FloatHi scr[F24][F42],Int4 window)132 static Int4 ReadPCC (CharPtr res, FloatHi scr[F24][F42], Int4 window)
133 {
134 FILE *fin;
135 Int2 i, n, val, count;
136 Int2 numcol;
137 FloatHi score;
138 Char buf[PATH_MAX], buff[PATH_MAX];
139 CharPtr bptr;
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 numcol = F7;
145
146 if (!FindPath ("ncbi", "ncbi", "data", buf, sizeof (buf)))
147 return 0;
148 FileBuildPath (buf, NULL, "KSpcc.dat");
149 if ((fin = FileOpen (buf, "r")) == NULL)
150 {
151 if ((fin = FileOpen ("KSpcc.dat", "r")) == NULL)
152 {
153 return 0;
154 }
155 }
156
157 val = 0;
158 while ((FileGets (buff, sizeof (buff), fin)) != NULL)
159 {
160 if (buff[0] == ';')
161 continue;
162 res[val] = c[val];
163 bptr = buff;
164 for (n = 0; n < numcol; n++)
165 {
166 sscanf (bptr, "%lf", &score);
167 scr[val][n] = score;
168 bptr += 7;
169 }
170 val++;
171 if (val == F24)
172 break;
173 }
174 count = val;
175
176 for (n = numcol; n < window; n++)
177 {
178 if (n == F42)
179 break;
180 val = (Int2) (n % numcol);
181 for (i = 0; i < F24; i++)
182 {
183 scr[i][n] = scr[i][val];
184 }
185 }
186
187 FileClose (fin);
188 return (Int4) count;
189 }
190
191 /* functions - external */
192
PCCProc(BioseqPtr bsp,SeqFeatPtr sfp,Int4 window)193 extern SeqGraphPtr PCCProc (BioseqPtr bsp, SeqFeatPtr sfp, Int4 window)
194 {
195 SeqIdPtr psip;
196 BioseqPtr pbsp;
197 Int4 start, end;
198 FloatHi scr[F24][F42];
199 Char res[F24];
200 SeqPortPtr spp;
201 SeqGraphPtr sgp = NULL;
202
203 if (bsp != NULL)
204 {
205 pbsp = bsp;
206
207 if (!ISA_aa (pbsp->mol))
208 return sgp;
209 }
210 else if (sfp != NULL)
211 {
212 if (sfp->data.choice != SEQFEAT_CDREGION)
213 return sgp;
214
215 psip = SeqLocId (sfp->product);
216 pbsp = BioseqFind (psip);
217 if (pbsp == NULL)
218 return sgp;
219
220 if (!ISA_aa (pbsp->mol))
221 return sgp;
222 }
223
224 if (ReadPCC (res, scr, window) != F24)
225 return sgp;
226
227 start = 0;
228 end = pbsp->length-1;
229 if ((spp = SeqPortNew (pbsp, start, end, 0, Seq_code_iupacaa)) != NULL)
230 {
231 end -= window;
232 sgp = PCCGraph (spp, start, end, window, scr, res);
233 SeqPortFree (spp);
234 }
235
236 return sgp;
237 }
238