1 /* @source stssearch application
2 **
3 ** Gribskov statistical plot of synonymous codon usage
4 **
5 ** @author Unknown
6 ** @@
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 "emboss.h"
24
25
26
27
28 /* @datastatic Primer *********************************************************
29 **
30 ** stssearch internals
31 **
32 ** @alias primers
33 **
34 ** @attr Name [AjPStr] Undocumented
35 ** @attr Prima [AjPRegexp] Undocumented
36 ** @attr Primb [AjPRegexp] Undocumented
37 ** @attr Oligoa [AjPStr] Undocumented
38 ** @attr Oligob [AjPStr] Undocumented
39 ******************************************************************************/
40
41 typedef struct primers
42 {
43 AjPStr Name;
44 AjPRegexp Prima;
45 AjPRegexp Primb;
46 AjPStr Oligoa;
47 AjPStr Oligob;
48 } *Primer;
49
50
51
52
53 AjPFile out = NULL;
54 AjPSeq seq = NULL;
55 AjPStr seqstr = NULL;
56 AjPStr revstr = NULL;
57 ajint nprimers = 0;
58 ajint ntests = 0;
59
60
61
62
63 static void stssearch_primDel(void **x, void *cl);
64 static void stssearch_primTest(void **x,void *cl);
65
66
67
68
69 /* @prog stssearch ************************************************************
70 **
71 ** Searches a DNA database for matches with a set of STS primers
72 **
73 ******************************************************************************/
74
main(int argc,char ** argv)75 int main(int argc, char **argv)
76 {
77
78 AjPSeqall seqall;
79 AjPFile primfile;
80 AjPStr rdline = NULL;
81
82 Primer primdata;
83 AjPStrTok handle = NULL;
84
85 AjPList primList = NULL;
86
87 embInit("stssearch", argc, argv);
88
89 primfile = ajAcdGetInfile("infile");
90 out = ajAcdGetOutfile("outfile");
91 seqall = ajAcdGetSeqall("seqall");
92
93 while(ajReadlineTrim(primfile, &rdline))
94 {
95 if(ajStrGetCharFirst(rdline) == '#')
96 continue;
97 if(ajStrSuffixC(rdline, ".."))
98 continue;
99
100 AJNEW(primdata);
101 primdata->Name = NULL;
102 primdata->Oligoa = NULL;
103 primdata->Oligob = NULL;
104
105 handle = ajStrTokenNewC(rdline, " \t");
106 ajStrTokenNextParse(handle, &primdata->Name);
107
108 if(!(nprimers % 1000))
109 ajDebug("Name [%d]: '%S'\n", nprimers, primdata->Name);
110
111 ajStrTokenNextParse(handle, &primdata->Oligoa);
112 ajStrFmtUpper(&primdata->Oligoa);
113 primdata->Prima = ajRegComp(primdata->Oligoa);
114
115 ajStrTokenNextParse(handle, &primdata->Oligob);
116 ajStrFmtUpper(&primdata->Oligob);
117 primdata->Primb = ajRegComp(primdata->Oligob);
118 ajStrTokenDel(&handle);
119
120 if(!nprimers)
121 primList = ajListNew();
122
123 ajListPushAppend(primList, primdata);
124 nprimers++;
125 }
126
127 if(!nprimers)
128 ajFatal("No primers read\n");
129
130 ajDebug("%d primers read\n", nprimers);
131
132 while(ajSeqallNext(seqall, &seq))
133 {
134 ajSeqFmtUpper(seq);
135 ajStrAssignS(&seqstr, ajSeqGetSeqS(seq));
136 ajStrAssignS(&revstr, ajSeqGetSeqS(seq));
137 ajSeqstrReverse(&revstr);
138 ajDebug("Testing: %s\n", ajSeqGetNameC(seq));
139 ntests = 0;
140 ajListMap(primList, &stssearch_primTest, NULL);
141 }
142
143 ajFileClose(&out);
144
145 ajSeqallDel(&seqall);
146 ajSeqDel(&seq);
147 ajFileClose(&out);
148 ajStrDel(&revstr);
149 ajStrDel(&seqstr);
150 ajFileClose(&primfile);
151 ajListMap(primList, &stssearch_primDel, NULL);
152 ajListFree(&primList);
153 ajStrDel(&rdline);
154
155
156 embExit();
157
158 return 0;
159 }
160
161
162
163
164 /* @funcstatic stssearch_primDel **********************************************
165 **
166 ** Undocumented.
167 **
168 ** @param [r] x [void**] Undocumented
169 ** @param [r] cl [void*] Undocumented
170 ** @return [void]
171 ** @@
172 ******************************************************************************/
173
174
stssearch_primDel(void ** x,void * cl)175 static void stssearch_primDel(void **x, void *cl)
176 {
177 Primer* p;
178 Primer primdata;
179
180 (void) cl; /* make it used */
181
182 p = (Primer*) x;
183 primdata = *p;
184
185 ajStrDel(&primdata->Name);
186 ajRegFree(&primdata->Prima);
187 ajRegFree(&primdata->Primb);
188 ajStrDel(&primdata->Oligoa);
189 ajStrDel(&primdata->Oligob);
190 AJFREE(*p);
191 }
192
193
194
195
196 /* @funcstatic stssearch_primTest *********************************************
197 **
198 ** Undocumented.
199 **
200 ** @param [r] x [void**] Undocumented
201 ** @param [r] cl [void*] Undocumented
202 ** @return [void]
203 ** @@
204 ******************************************************************************/
205
206
stssearch_primTest(void ** x,void * cl)207 static void stssearch_primTest(void **x,void *cl)
208 {
209 Primer* p;
210 Primer primdata;
211
212 AjBool testa;
213 AjBool testb;
214 AjBool testc;
215 AjBool testd;
216 ajint ioff;
217
218 (void) cl; /* make it used */
219
220 p = (Primer*) x;
221 primdata = *p;
222
223 ntests++;
224
225 if(!(ntests % 1000))
226 ajDebug("completed tests: %d\n", ntests);
227
228 testa = ajRegExec(primdata->Prima, seqstr);
229
230 if(testa)
231 {
232 ioff = ajRegOffset(primdata->Prima);
233 ajDebug("%s: %S PrimerA matched at %d\n",
234 ajSeqGetNameC(seq), primdata->Name, ioff);
235 ajFmtPrintF(out, "%s: %S PrimerA matched at %d\n",
236 ajSeqGetNameC(seq), primdata->Name, ioff);
237 ajRegTrace(primdata->Prima);
238 }
239
240 testb = ajRegExec(primdata->Primb, seqstr);
241 if(testb)
242 {
243 ioff = ajRegOffset(primdata->Primb);
244 ajDebug("%s: %S PrimerB matched at %d\n",
245 ajSeqGetNameC(seq), primdata->Name, ioff);
246 ajFmtPrintF(out, "%s: %S PrimerB matched at %d\n",
247 ajSeqGetNameC(seq), primdata->Name, ioff);
248 ajRegTrace(primdata->Primb);
249 }
250
251 testc = ajRegExec(primdata->Prima, revstr);
252 if(testc)
253 {
254 ioff = ajStrGetLen(seqstr) - ajRegOffset(primdata->Prima);
255 ajDebug("%s: (rev) %S PrimerA matched at %d\n",
256 ajSeqGetNameC(seq), primdata->Name, ioff);
257 ajFmtPrintF(out, "%s: (rev) %S PrimerA matched at %d\n",
258 ajSeqGetNameC(seq), primdata->Name, ioff);
259 ajRegTrace(primdata->Prima);
260 }
261
262 testd = ajRegExec(primdata->Primb, revstr);
263 if(testd)
264 {
265 ioff = ajStrGetLen(seqstr) - ajRegOffset(primdata->Primb);
266 ajDebug("%s: (rev) %S PrimerB matched at %d\n",
267 ajSeqGetNameC(seq), primdata->Name, ioff);
268 ajFmtPrintF(out, "%s: (rev) %S PrimerB matched at %d\n",
269 ajSeqGetNameC(seq), primdata->Name, ioff);
270 ajRegTrace(primdata->Primb);
271 }
272
273 return;
274 }
275