1 /* @source eprimer3 application
2 ** Picks PCR primers and hybridization oligos
3 **
4 ** @author Copyright (C) Gary Williams (gwilliam@hgmp.mrc.ac.uk)
5 ** 5 Nov 2001 - GWW - written
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 #ifndef WIN32
25 #include <unistd.h>
26 #else
27 #include <io.h>
28 #include <fcntl.h>
29 #include <windows.h>
30 #define fdopen _fdopen
31 #endif
32
33
34 /* estimate size of a sequence's output table */
35 #define TABLEGUESS 251
36
37
38
39
40 static FILE* eprimer3_start_write(int fd);
41 static void eprimer3_write(const AjPStr str, FILE *stream);
42 static void eprimer3_end_write(FILE *stream);
43 static void eprimer3_read(int fd, AjPStr * result);
44 static void eprimer3_send_range(FILE * stream, const char * tag,
45 const AjPRange value,
46 ajint begin);
47 static void eprimer3_send_range2(FILE * stream, const char * tag,
48 const AjPRange value);
49 static void eprimer3_send_int(FILE * stream, const char * tag, ajint value);
50 static void eprimer3_send_float(FILE * stream, const char * tag, float value);
51 static void eprimer3_send_bool(FILE * stream, const char * tag, AjBool value);
52 static void eprimer3_send_string(FILE * stream, const char * tag,
53 const AjPStr value);
54 static void eprimer3_send_stringC(FILE * stream, const char * tag,
55 const char * value);
56 static void eprimer3_send_end(FILE * stream);
57 static void eprimer3_report (AjPFile outfile, const AjPStr output,
58 ajint numreturn, ajint begin);
59 static void eprimer3_output_report(AjPFile outfile, const AjPTable table,
60 ajint numreturn, ajint begin);
61 static const AjPStr eprimer3_tableget(const char * key1, ajint number,
62 const char *key2,
63 const AjPTable table);
64 static void eprimer3_write_primer(AjPFile outfile, const char *tag,
65 const AjPStr pos,
66 const AjPStr tm, const AjPStr gc,
67 const AjPStr seq,
68 AjBool rev, ajint begin);
69
70
71
72
73 /* @prog eprimer3 *************************************************************
74 **
75 ** EMBOSS wrapper for the Whitehead's primer3 program
76 **
77 ******************************************************************************/
78
main(int argc,char ** argv)79 int main(int argc, char **argv)
80 {
81 /* Global details */
82 AjBool explain_flag;
83 AjBool file_flag;
84 AjPStr* task;
85 AjBool do_primer;
86 AjBool do_hybrid;
87 ajint num_return;
88 ajint first_base_index;
89
90 /* "Sequence" Input Tags */
91 AjPSeqall sequence;
92 AjPRange included_region;
93 AjPRange target;
94 AjPRange excluded_region;
95 AjPStr left_input;
96 AjPStr right_input;
97
98 /* Primer details */
99 AjBool pick_anyway;
100 AjPFile mispriming_library;
101 float max_mispriming;
102 float pair_max_mispriming;
103 ajint gc_clamp;
104 ajint opt_size;
105 ajint min_size;
106 ajint max_size;
107 float opt_tm;
108 float min_tm;
109 float max_tm;
110 float max_diff_tm;
111 float opt_gc_percent;
112 float min_gc;
113 float max_gc;
114 float salt_conc;
115 float dna_conc;
116 ajint num_ns_accepted;
117 float self_any;
118 float self_end;
119 ajint max_poly_x;
120
121 /* Sequence Quality. These are not (yet) implemented */
122 /*
123 AjPFile sequence_quality;
124 ajint min_quality;
125 ajint min_end_quality;
126 ajint quality_range_min;
127 ajint quality_range_max;
128 */
129
130 /* Product details */
131 ajint product_opt_size;
132 AjPRange product_size_range;
133 float product_opt_tm;
134 float product_min_tm;
135 float product_max_tm;
136
137 /* Objective Function Penalty Weights for Primers */
138 float max_end_stability;
139
140 /* these are not (yet) implemented */
141 /*
142 float inside_penalty;
143 float outside_penalty;
144 */
145
146 /* Primer penalties */
147 /* these are not (yet) implemented */
148
149 /* Internal Oligo "Sequence" Input Tags */
150 AjPRange internal_oligo_excluded_region;
151
152 /* Internal Oligo "Global" Input Tags */
153 AjPStr internal_oligo_input;
154 ajint internal_oligo_opt_size;
155 ajint internal_oligo_min_size;
156 ajint internal_oligo_max_size;
157 float internal_oligo_opt_tm;
158 float internal_oligo_min_tm;
159 float internal_oligo_max_tm;
160 float internal_oligo_opt_gc_percent;
161 float internal_oligo_min_gc;
162 float internal_oligo_max_gc;
163 float internal_oligo_salt_conc;
164 float internal_oligo_dna_conc;
165 float internal_oligo_self_any;
166 float internal_oligo_self_end;
167 ajint internal_oligo_max_poly_x;
168 AjPFile internal_oligo_mishyb_library;
169 float internal_oligo_max_mishyb;
170
171 /*
172 ajint internal_oligo_min_quality;
173 */
174
175 /* Internal Oligo penalties */
176 /* these are not (yet) implemented */
177
178 /* EMBOSS-wrapper-specific stuff */
179 AjPFile outfile;
180
181 /* other variables */
182 AjPStr result = NULL;
183 AjPStr strand = NULL;
184 AjPStr substr = NULL;
185 AjPSeq seq = NULL;
186 ajint begin = 0;
187 ajint end;
188 FILE* stream;
189 AjPStr taskstr = NULL;
190 const AjPStr program = NULL;
191
192 /* pipe variables */
193
194 int *pipeto; /* pipe to feed the exec'ed program input */
195 int *pipefrom; /* pipe to get the exec'ed program output */
196
197 embInit("eprimer3", argc, argv);
198
199 /* Global details */
200 explain_flag = ajAcdGetBoolean("explainflag");
201 file_flag = ajAcdGetBoolean("fileflag");
202 task = ajAcdGetList("task");
203 do_primer = ajAcdGetToggle("primer");
204 do_hybrid = ajAcdGetToggle("hybridprobe");
205 num_return = ajAcdGetInt("numreturn");
206 first_base_index = ajAcdGetInt("firstbaseindex");
207
208 /* "Sequence" Input Tags */
209 sequence = ajAcdGetSeqall("sequence");
210 included_region = ajAcdGetRange("includedregion");
211 target = ajAcdGetRange("targetregion");
212 excluded_region = ajAcdGetRange("excludedregion");
213 left_input = ajAcdGetString("forwardinput");
214 right_input = ajAcdGetString("reverseinput");
215
216 /* Primer details */
217 pick_anyway = ajAcdGetBoolean("pickanyway");
218 mispriming_library = ajAcdGetInfile("mispriminglibraryfile");
219 max_mispriming = ajAcdGetFloat("maxmispriming");
220 pair_max_mispriming = ajAcdGetFloat("pairmaxmispriming");
221 gc_clamp = ajAcdGetInt("gcclamp");
222 opt_size = ajAcdGetInt("optsize");
223 min_size = ajAcdGetInt("minsize");
224 max_size = ajAcdGetInt("maxsize");
225 opt_tm = ajAcdGetFloat("opttm");
226 min_tm = ajAcdGetFloat("mintm");
227 max_tm = ajAcdGetFloat("maxtm");
228 max_diff_tm = ajAcdGetFloat("maxdifftm");
229 opt_gc_percent = ajAcdGetFloat("ogcpercent");
230 min_gc = ajAcdGetFloat("mingc");
231 max_gc = ajAcdGetFloat("maxgc");
232 salt_conc = ajAcdGetFloat("saltconc");
233 dna_conc = ajAcdGetFloat("dnaconc");
234 num_ns_accepted = ajAcdGetInt("numnsaccepted");
235 self_any = ajAcdGetFloat("selfany");
236 self_end = ajAcdGetFloat("selfend");
237 max_poly_x = ajAcdGetInt("maxpolyx");
238
239 AJCNEW0(pipeto,2);
240 AJCNEW0(pipefrom,2);
241
242 /* Sequence Quality */
243 /* these are not (yet) implemented */
244 /*
245 sequence_quality = ajAcdGetInfile("sequencequality");
246 min_quality = ajAcdGetInt("minquality");
247 min_end_quality = ajAcdGetInt("minendquality");
248 quality_range_min = ajAcdGetInt("qualityrangemin");
249 quality_range_max = ajAcdGetInt("qualityrangemax");
250 */
251
252 /* Product details */
253 product_opt_size = ajAcdGetInt("psizeopt");
254 product_size_range = ajAcdGetRange("prange");
255 product_opt_tm = ajAcdGetFloat("ptmopt");
256 product_min_tm = ajAcdGetFloat("ptmmin");
257 product_max_tm = ajAcdGetFloat("ptmmax");
258
259 /* Objective Function Penalty Weights for Primers */
260 max_end_stability = ajAcdGetFloat("maxendstability");
261 /* these are not (yet) implemented */
262 /*
263 inside_penalty = ajAcdGetFloat("insidepenalty");
264 outside_penalty = ajAcdGetFloat("outsidepenalty");
265 */
266
267 /* Primer penalties */
268 /* these are not (yet) implemented */
269
270 /* Internal Oligo "Sequence" Input Tags */
271 internal_oligo_excluded_region = ajAcdGetRange("oexcludedregion");
272 internal_oligo_input = ajAcdGetString("oligoinput");
273
274 /* Internal Oligo "Global" Input Tags */
275 internal_oligo_opt_size = ajAcdGetInt("osizeopt");
276 internal_oligo_min_size = ajAcdGetInt("ominsize");
277 internal_oligo_max_size = ajAcdGetInt("omaxsize");
278 internal_oligo_opt_tm = ajAcdGetFloat("otmopt");
279 internal_oligo_min_tm = ajAcdGetFloat("otmmin");
280 internal_oligo_max_tm = ajAcdGetFloat("otmmax");
281 internal_oligo_opt_gc_percent = ajAcdGetFloat("ogcopt");
282 internal_oligo_min_gc = ajAcdGetFloat("ogcmin");
283 internal_oligo_max_gc = ajAcdGetFloat("ogcmax");
284 internal_oligo_salt_conc = ajAcdGetFloat("osaltconc");
285 internal_oligo_dna_conc = ajAcdGetFloat("odnaconc");
286 internal_oligo_self_any = ajAcdGetFloat("oanyself");
287 internal_oligo_self_end = ajAcdGetFloat("oendself");
288 internal_oligo_max_poly_x = ajAcdGetInt("opolyxmax");
289 internal_oligo_mishyb_library = ajAcdGetInfile("mishyblibraryfile");
290 internal_oligo_max_mishyb = ajAcdGetFloat("omishybmax");
291 /*
292 internal_oligo_min_quality = ajAcdGetInt("oligominquality");
293 */
294
295 /* Internal Oligo penalties */
296 /* these are not (yet) implemented */
297
298
299 /* EMBOSS-wrapper-specific stuff */
300 outfile = ajAcdGetOutfile("outfile");
301
302
303 ajStrRemoveWhite(&left_input);
304 ajStrRemoveWhite(&right_input);
305
306 /*
307 ** OK - we will now try to do a separate fork-exec for each sequence.
308 */
309
310 result = ajStrNew();
311
312 while(ajSeqallNext(sequence, &seq))
313 {
314 program = ajAcdGetpathC("primer3_core");
315
316 if(!ajSysExecRedirectC(ajStrGetPtr(program),&pipeto,&pipefrom))
317 ajFatal("eprimer3: Could not exec primer3_core");
318
319 stream = eprimer3_start_write(pipeto[1]);
320
321 /* send primer3 Primer "Global" parameters */
322 eprimer3_send_bool(stream, "PRIMER_EXPLAIN_FLAG", explain_flag);
323 eprimer3_send_bool(stream, "PRIMER_FILE_FLAG", file_flag);
324
325 if(do_hybrid)
326 {
327 if(!ajStrCmpC(task[0], "1"))
328 ajStrAssignC(&taskstr, "pick_pcr_primers_and_hyb_probe");
329 else if(!ajStrCmpC(task[0], "2"))
330 ajStrAssignC(&taskstr, "pick_left_only");
331 else if(!ajStrCmpC(task[0], "3"))
332 ajStrAssignC(&taskstr, "pick_right_only");
333 else if(!ajStrCmpC(task[0], "4"))
334 ajStrAssignC(&taskstr, "pick_hyb_probe_only");
335
336 if (!do_primer)
337 ajStrAssignC(&taskstr, "pick_hyb_probe_only");
338 }
339 else
340 {
341 if(!ajStrCmpC(task[0], "1"))
342 ajStrAssignC(&taskstr, "pick_pcr_primers");
343 else if(!ajStrCmpC(task[0], "2"))
344 ajStrAssignC(&taskstr, "pick_left_only");
345 else if(!ajStrCmpC(task[0], "3"))
346 ajStrAssignC(&taskstr, "pick_right_only");
347 else if(!ajStrCmpC(task[0], "4"))
348 ajStrAssignC(&taskstr, "pick_hyb_probe_only");
349 }
350
351 eprimer3_send_string(stream, "PRIMER_TASK", taskstr);
352 eprimer3_send_int(stream, "PRIMER_NUM_RETURN", num_return);
353 eprimer3_send_int(stream, "PRIMER_FIRST_BASE_INDEX",
354 first_base_index);
355 eprimer3_send_bool(stream, "PRIMER_PICK_ANYWAY", pick_anyway);
356
357 /* mispriming library may not have been specified */
358 if(mispriming_library)
359 eprimer3_send_stringC(stream, "PRIMER_MISPRIMING_LIBRARY",
360 ajFileGetPrintnameC(mispriming_library));
361
362 eprimer3_send_float(stream, "PRIMER_MAX_MISPRIMING",
363 max_mispriming);
364 eprimer3_send_float(stream, "PRIMER_PAIR_MAX_MISPRIMING",
365 pair_max_mispriming);
366 eprimer3_send_int(stream, "PRIMER_GC_CLAMP", gc_clamp);
367 eprimer3_send_int(stream, "PRIMER_OPT_SIZE", opt_size);
368 eprimer3_send_int(stream, "PRIMER_MIN_SIZE", min_size);
369 eprimer3_send_int(stream, "PRIMER_MAX_SIZE", max_size);
370 eprimer3_send_float(stream, "PRIMER_OPT_TM", opt_tm);
371 eprimer3_send_float(stream, "PRIMER_MIN_TM", min_tm);
372 eprimer3_send_float(stream, "PRIMER_MAX_TM", max_tm);
373 eprimer3_send_float(stream, "PRIMER_MAX_DIFF_TM", max_diff_tm);
374 eprimer3_send_float(stream, "PRIMER_OPT_GC_PERCENT",
375 opt_gc_percent);
376 eprimer3_send_float(stream, "PRIMER_MIN_GC", min_gc);
377 eprimer3_send_float(stream, "PRIMER_MAX_GC", max_gc);
378 eprimer3_send_float(stream, "PRIMER_SALT_CONC", salt_conc);
379 eprimer3_send_float(stream, "PRIMER_DNA_CONC", dna_conc);
380 eprimer3_send_int(stream, "PRIMER_NUM_NS_ACCEPTED",
381 num_ns_accepted);
382 eprimer3_send_float(stream, "PRIMER_SELF_ANY", self_any);
383 eprimer3_send_float(stream, "PRIMER_SELF_END", self_end);
384 eprimer3_send_int(stream, "PRIMER_MAX_POLY_X", max_poly_x);
385 eprimer3_send_int(stream, "PRIMER_PRODUCT_OPT_SIZE",
386 product_opt_size);
387 eprimer3_send_range2(stream, "PRIMER_PRODUCT_SIZE_RANGE",
388 product_size_range);
389 eprimer3_send_float(stream, "PRIMER_PRODUCT_OPT_TM",
390 product_opt_tm);
391 eprimer3_send_float(stream, "PRIMER_PRODUCT_MIN_TM",
392 product_min_tm);
393 eprimer3_send_float(stream, "PRIMER_PRODUCT_MAX_TM",
394 product_max_tm);
395 eprimer3_send_float(stream, "PRIMER_MAX_END_STABILITY",
396 max_end_stability);
397
398 /* send primer3 Internal Oligo "Global" parameters */
399 eprimer3_send_int(stream, "PRIMER_INTERNAL_OLIGO_OPT_SIZE",
400 internal_oligo_opt_size);
401 eprimer3_send_int(stream, "PRIMER_INTERNAL_OLIGO_MIN_SIZE",
402 internal_oligo_min_size);
403 eprimer3_send_int(stream, "PRIMER_INTERNAL_OLIGO_MAX_SIZE",
404 internal_oligo_max_size);
405 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_OPT_TM",
406 internal_oligo_opt_tm);
407 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_MIN_TM",
408 internal_oligo_min_tm);
409 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_MAX_TM",
410 internal_oligo_max_tm);
411 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_OPT_GC_PERCENT",
412 internal_oligo_opt_gc_percent);
413 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_MIN_GC",
414 internal_oligo_min_gc);
415 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_MAX_GC",
416 internal_oligo_max_gc);
417 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_SALT_CONC",
418 internal_oligo_salt_conc);
419 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_DNA_CONC",
420 internal_oligo_dna_conc);
421 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_SELF_ANY",
422 internal_oligo_self_any);
423 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_SELF_END",
424 internal_oligo_self_end);
425 eprimer3_send_int(stream, "PRIMER_INTERNAL_OLIGO_MAX_POLY_X",
426 internal_oligo_max_poly_x);
427
428 /*
429 ** internal oligo mishybridising library may not have been
430 ** specified
431 */
432 if(internal_oligo_mishyb_library)
433 eprimer3_send_stringC(stream,
434 "PRIMER_INTERNAL_OLIGO_MISHYB_LIBRARY",
435 ajFileGetPrintnameC(internal_oligo_mishyb_library));
436
437 eprimer3_send_float(stream, "PRIMER_INTERNAL_OLIGO_MAX_MISHYB",
438 internal_oligo_max_mishyb);
439
440
441 /*
442 ** Start sequence-specific stuff
443 */
444
445 begin = ajSeqallGetseqBegin(sequence) - 1;
446 end = ajSeqallGetseqEnd(sequence) - 1;
447
448 strand = ajSeqGetSeqCopyS(seq);
449
450 ajStrFmtUpper(&strand);
451 ajStrAssignSubC(&substr,ajStrGetPtr(strand), begin, end);
452
453 /* send flags to turn on using optimal product size */
454 eprimer3_send_float(stream, "PRIMER_PAIR_WT_PRODUCT_SIZE_GT",
455 (float)0.05);
456 eprimer3_send_float(stream, "PRIMER_PAIR_WT_PRODUCT_SIZE_LT",
457 (float)0.05);
458
459 /* send primer3 Primer "Sequence" parameters */
460 eprimer3_send_string(stream, "SEQUENCE", substr);
461
462 /* if no ID name, use the USA */
463 if(ajStrMatchC(ajSeqGetNameS(seq),""))
464 eprimer3_send_string(stream, "PRIMER_SEQUENCE_ID",
465 ajSeqGetUsaS(seq));
466 else
467 eprimer3_send_string(stream, "PRIMER_SEQUENCE_ID",
468 ajSeqGetNameS(seq));
469
470 eprimer3_send_range(stream, "INCLUDED_REGION", included_region,
471 begin);
472 eprimer3_send_range(stream, "TARGET", target, begin);
473 eprimer3_send_range(stream, "EXCLUDED_REGION", excluded_region,
474 begin);
475 eprimer3_send_string(stream, "PRIMER_LEFT_INPUT", left_input);
476 eprimer3_send_string(stream, "PRIMER_RIGHT_INPUT", right_input);
477
478 /* send primer3 Internal Oligo "Sequence" parameters */
479 eprimer3_send_range(stream,
480 "PRIMER_INTERNAL_OLIGO_EXCLUDED_REGION",
481 internal_oligo_excluded_region, begin);
482 eprimer3_send_string(stream, "PRIMER_INTERNAL_OLIGO_INPUT",
483 internal_oligo_input);
484
485
486 /* end the primer3 input sequence record with a '=' */
487 eprimer3_send_end(stream);
488 /* and close the ouput pipe stream */
489 eprimer3_end_write(stream);
490
491 /* read the primer3 output */
492 eprimer3_read(pipefrom[0], &result);
493
494 eprimer3_report(outfile, result, num_return, begin);
495
496 ajStrSetClear(&result);
497
498 #ifndef WIN32
499 close(pipeto[1]);
500 close(pipefrom[0]);
501 #endif
502 } /* end of sequence loop */
503
504
505 ajStrDel(&result);
506 ajSeqDel(&seq);
507 ajStrDel(&strand);
508 ajStrDel(&substr);
509 ajFileClose(&outfile);
510 ajStrDel(&taskstr);
511 ajStrDelarray(&task);
512
513 ajSeqallDel(&sequence);
514 ajSeqDel(&seq);
515
516 ajRangeDel(&included_region);
517 ajRangeDel(&target);
518 ajRangeDel(&excluded_region);
519 ajRangeDel(&product_size_range);
520 ajRangeDel(&internal_oligo_excluded_region);
521
522 ajStrDel(&left_input);
523 ajStrDel(&right_input);
524 ajStrDel(&internal_oligo_input);
525
526 AJFREE(pipeto);
527 AJFREE(pipefrom);
528
529 ajFileClose(&mispriming_library);
530
531 embExit();
532
533 return 0;
534 }
535
536
537
538
539 /* @funcstatic eprimer3_read **************************************************
540 **
541 ** Reads the output from primer3_core into a returned AjPStr until EOF
542 **
543 ** @param [r] fd [int] file descriptor
544 ** @param [u] result [AjPStr*] Returned string
545 ** @return [void]
546 **
547 ******************************************************************************/
548
eprimer3_read(int fd,AjPStr * result)549 static void eprimer3_read(int fd, AjPStr* result)
550 {
551 FILE *stream;
552 int ch;
553
554 ajDebug("reading primer3_core output\n");
555
556 if((stream = fdopen(fd, "r")) == NULL)
557 ajFatal("fdopen() read error");
558
559 while((ch = getc( stream )) != EOF)
560 ajStrAppendK(result, ch);
561
562 ajDebug("primer3_core output:\n%S\n", *result);
563
564
565 fclose(stream);
566 ajDebug("reading done\n");
567
568 return;
569 }
570
571
572
573
574 /* @funcstatic eprimer3_send_end **********************************************
575 **
576 ** Writes the end-of-input flag '=' to the input stream of primer3_core
577 **
578 ** @param [u] stream [FILE *] File handle
579 ** @return [void]
580 **
581 ******************************************************************************/
582
eprimer3_send_end(FILE * stream)583 static void eprimer3_send_end(FILE * stream)
584 {
585 fputs( "=\n", stream );
586
587 return;
588 }
589
590
591
592
593 /* @funcstatic eprimer3_send_range ********************************************
594 **
595 ** Write range data to primer3_core as 'start,length start2,length2,etc'
596 **
597 ** @param [u] stream [FILE *] File handle
598 ** @param [r] tag [const char *] Tag of primer3 data type
599 ** @param [r] value [const AjPRange] Ranges to write
600 ** @param [r] begin [ajint] Start position of subsequence (-sbegin)
601 ** @return [void]
602 **
603 ******************************************************************************/
604
eprimer3_send_range(FILE * stream,const char * tag,const AjPRange value,ajint begin)605 static void eprimer3_send_range(FILE * stream, const char * tag,
606 const AjPRange value,
607 ajint begin)
608 {
609 AjPStr str;
610 ajuint n;
611 ajuint start;
612 ajuint end;
613
614 str = ajStrNew();
615
616 if(ajRangeGetSize(value))
617 {
618 ajFmtPrintS(&str, "%s=", tag);
619 eprimer3_write(str, stream);
620 ajStrSetClear(&str);
621
622 for(n=0; n < ajRangeGetSize(value); n++)
623 {
624 ajRangeElementGetValues(value, n, &start, &end);
625 start -= begin;
626 end -= begin;
627 ajFmtPrintS(&str, "%d,%d ", start, end-start+1);
628 eprimer3_write(str, stream);
629 ajStrSetClear(&str);
630 }
631
632 ajFmtPrintS(&str, "\n");
633 eprimer3_write(str, stream);
634 }
635
636 ajStrDel(&str);
637
638 return;
639 }
640
641
642
643
644 /* @funcstatic eprimer3_send_range2 *******************************************
645 **
646 ** Write alternate display of ranges as 'a-b c-d' to primer3_core
647 **
648 ** @param [u] stream [FILE *] File handle
649 ** @param [r] tag [const char *] Tag of primer3 data type
650 ** @param [r] value [const AjPRange] Ranges to write
651 ** @return [void]
652 **
653 ******************************************************************************/
654
eprimer3_send_range2(FILE * stream,const char * tag,const AjPRange value)655 static void eprimer3_send_range2(FILE * stream, const char * tag,
656 const AjPRange value)
657 {
658 AjPStr str;
659 ajuint n;
660 ajuint start;
661 ajuint end;
662
663 str=ajStrNew();
664
665 if(ajRangeGetSize(value))
666 {
667 ajFmtPrintS(&str, "%s=", tag);
668 eprimer3_write(str, stream);
669 ajStrSetClear(&str);
670
671 for(n=0; n < ajRangeGetSize(value); n++)
672 {
673 ajRangeElementGetValues(value, n, &start, &end);
674 ajFmtPrintS(&str, "%d-%d ", start, end);
675 eprimer3_write(str, stream);
676 ajStrSetClear(&str);
677 }
678
679 ajFmtPrintS(&str, "\n");
680 eprimer3_write(str, stream);
681 }
682
683 ajStrDel(&str);
684
685 return;
686 }
687
688
689
690
691 /* @funcstatic eprimer3_send_int **********************************************
692 **
693 ** Write integer to primer3_core
694 **
695 ** @param [u] stream [FILE *] File handle
696 ** @param [r] tag [const char *] Tag of primer3 data type
697 ** @param [r] value [ajint] Integer value to write
698 ** @return [void]
699 **
700 ******************************************************************************/
701
eprimer3_send_int(FILE * stream,const char * tag,ajint value)702 static void eprimer3_send_int(FILE * stream, const char * tag, ajint value)
703 {
704 AjPStr str;
705
706 str = ajStrNew();
707
708 ajFmtPrintS(&str, "%s=%d\n", tag, value);
709 eprimer3_write(str, stream);
710
711 ajStrDel(&str);
712
713 return;
714 }
715
716
717
718
719 /* @funcstatic eprimer3_send_float ********************************************
720 **
721 ** Write float to primer3_core
722 **
723 ** @param [u] stream [FILE *] File handle
724 ** @param [r] tag [const char *] Tag of primer3 data type
725 ** @param [r] value [float] Float value to write
726 ** @return [void]
727 **
728 ******************************************************************************/
729
eprimer3_send_float(FILE * stream,const char * tag,float value)730 static void eprimer3_send_float(FILE * stream, const char * tag, float value)
731 {
732 AjPStr str;
733
734 str = ajStrNew();
735
736 ajFmtPrintS(&str, "%s=%f\n", tag, value);
737 eprimer3_write(str, stream);
738
739 ajStrDel(&str);
740
741 return;
742 }
743
744
745
746
747 /* @funcstatic eprimer3_send_bool *********************************************
748 **
749 ** Write boolean to primer3_core
750 **
751 ** @param [u] stream [FILE *] File handle
752 ** @param [r] tag [const char *] Tag of primer3 data type
753 ** @param [r] value [AjBool] Boolean value to write
754 ** @return [void]
755 **
756 ******************************************************************************/
757
eprimer3_send_bool(FILE * stream,const char * tag,AjBool value)758 static void eprimer3_send_bool(FILE * stream, const char * tag, AjBool value)
759 {
760 AjPStr str;
761
762 str = ajStrNew();
763
764 ajFmtPrintS(&str, "%s=%d\n", tag, value?1:0);
765 eprimer3_write(str, stream);
766
767 ajStrDel(&str);
768
769 return;
770 }
771
772
773
774
775 /* @funcstatic eprimer3_send_string *******************************************
776 **
777 ** Write string to primer3_core
778 **
779 **
780 ** @param [u] stream [FILE *] File handle
781 ** @param [r] tag [const char *] Tag of primer3 data type
782 ** @param [r] value [const AjPStr] String value to write
783 ** @return [void]
784 **
785 ******************************************************************************/
786
eprimer3_send_string(FILE * stream,const char * tag,const AjPStr value)787 static void eprimer3_send_string(FILE * stream, const char * tag,
788 const AjPStr value)
789 {
790 AjPStr str;
791
792 str = ajStrNew();
793
794 if(ajStrGetLen(value))
795 {
796 ajFmtPrintS(&str, "%s=%S\n", tag, value);
797 eprimer3_write(str, stream);
798 }
799
800 ajStrDel(&str);
801
802 return;
803 }
804
805
806
807
808 /* @funcstatic eprimer3_send_stringC ******************************************
809 **
810 ** Write char * to primer3_core
811 **
812 ** @param [u] stream [FILE*] File handle
813 ** @param [r] tag [const char*] Tag of primer3 data type
814 ** @param [r] value [const char*] Char * value to write
815 ** @return [void]
816 **
817 ******************************************************************************/
818
eprimer3_send_stringC(FILE * stream,const char * tag,const char * value)819 static void eprimer3_send_stringC(FILE *stream, const char *tag,
820 const char *value)
821 {
822 AjPStr str;
823
824 str = ajStrNew();
825
826 if(strlen(value))
827 {
828 ajFmtPrintS(&str, "%s=%s\n", tag, value);
829 eprimer3_write(str, stream);
830 }
831
832 ajStrDel(&str);
833
834 return;
835 }
836
837
838
839
840 /* @funcstatic eprimer3_start_write *******************************************
841 **
842 ** Open a file descriptor as a stream to pipe to primer3_core
843 **
844 ** @param [r] fd [int] File descriptor
845 ** @return [FILE*] File stream
846 **
847 ******************************************************************************/
848
eprimer3_start_write(int fd)849 static FILE* eprimer3_start_write(int fd)
850 {
851 FILE *stream;
852
853 ajDebug( "start writing\n");
854 if((stream = fdopen( fd, "w" )) == NULL)
855 ajFatal("Couldn't open pipe with fdopen()");
856
857 return stream;
858 }
859
860
861
862
863 /* @funcstatic eprimer3_write *************************************************
864 **
865 ** Write a tag=value AjPStr to the primer3_core input stream
866 **
867 ** @param [r] str [const AjPStr] Input string
868 ** @param [u] stream [FILE*] Stream piped to primer3_core
869 ** @return [void]
870 **
871 ******************************************************************************/
872
eprimer3_write(const AjPStr str,FILE * stream)873 static void eprimer3_write(const AjPStr str, FILE *stream)
874 {
875 fputs(ajStrGetPtr(str), stream);
876 ajDebug("eprimer3_write '%S'\n", str);
877
878 return;
879 }
880
881
882
883
884 /* @funcstatic eprimer3_end_write *********************************************
885 **
886 ** Close the stream piping in to primer3_core
887 **
888 ** @param [u] stream [FILE *] Stream
889 ** @return [void]
890 **
891 ******************************************************************************/
892
eprimer3_end_write(FILE * stream)893 static void eprimer3_end_write(FILE *stream)
894 {
895 ajDebug("Closing output pipe stream\n");
896 fclose(stream);
897
898 return;
899 }
900
901
902
903
904 /* @funcstatic eprimer3_report ************************************************
905 **
906 ** Read output of primer3_core into a temporary table of tag/value results
907 **
908 ** @param [u] outfile [AjPFile] Report outfile
909 ** @param [r] output [const AjPStr] Output from primer3_core
910 ** @param [r] numreturn [ajint] Number of results to return for each sequence
911 ** @param [r] begin [ajint] Start position of subsequence (-sbegin)
912 ** @return [void]
913 **
914 ******************************************************************************/
915
eprimer3_report(AjPFile outfile,const AjPStr output,ajint numreturn,ajint begin)916 static void eprimer3_report(AjPFile outfile, const AjPStr output,
917 ajint numreturn, ajint begin)
918 {
919 AjPStr line = NULL;
920 AjPStrTok linetokenhandle;
921 char eol[] = "\n\r";
922 AjPStrTok keytokenhandle;
923 char equals[] = "=";
924 AjPStr key = NULL;
925 AjPStr value = NULL;
926 AjBool gotsequenceid = AJFALSE;
927 AjPTable table;
928
929 linetokenhandle = ajStrTokenNewC(output, eol);
930
931 /* get next line of relevant results */
932 while(ajStrTokenNextParseC(linetokenhandle, eol, &line))
933 {
934 if(!gotsequenceid)
935 {
936 /*
937 ** Att the start of another sequence's results?
938 ** Start storing the results in the table.
939 */
940
941 if(ajStrCmpLenC(line, "PRIMER_SEQUENCE_ID=", 19) == 0)
942 {
943 gotsequenceid = AJTRUE;
944 table = ajTablestrNew(TABLEGUESS);
945
946 }
947 else
948 continue;
949 }
950 else
951 {
952 /*
953 ** At the end of this sequence? - marked by a '=' in the primer3
954 ** output - then output the results.
955 */
956 if(ajStrCmpC(line, "=") == 0)
957 {
958 gotsequenceid = AJFALSE;
959 eprimer3_output_report(outfile, table, numreturn, begin);
960 ajTablestrFree(&table);
961 continue;
962 }
963 }
964
965 /*
966 ** store key and value in table and parse values
967 ** when have all of the sequence
968 ** results in the table because the LEFT, RIGHT
969 ** and INTERNAL results for each
970 ** resulting primer are interleaved
971 */
972
973 keytokenhandle = ajStrTokenNewC(line, equals);
974
975 key = ajStrNew();
976 ajStrTokenNextParse(keytokenhandle, &key);
977
978 value = ajStrNew();
979 ajStrTokenNextParse(keytokenhandle, &value);
980
981 ajDebug("key=%S\tvalue=%S\n", key, value);
982
983 ajTablePut(table,(void *)key, (void *)value);
984
985 ajStrTokenDel(&keytokenhandle);
986 }
987
988 ajStrDel(&line);
989 ajStrTokenDel(&linetokenhandle);
990 ajTablestrFree(&table);
991
992 return;
993 }
994
995
996
997
998 /* @funcstatic eprimer3_output_report *****************************************
999 **
1000 ** Read the results out of the tag/value table and write to report
1001 **
1002 ** @param [u] outfile [AjPFile] Report outfile
1003 ** @param [r] table [const AjPTable] Table of tag/value result pairs
1004 ** @param [r] numreturn [ajint] Number of results to return for each sequence
1005 ** @param [r] begin [ajint] Start position of subsequence (-sbegin)
1006 ** @return [void]
1007 **
1008 ******************************************************************************/
1009
eprimer3_output_report(AjPFile outfile,const AjPTable table,ajint numreturn,ajint begin)1010 static void eprimer3_output_report(AjPFile outfile, const AjPTable table,
1011 ajint numreturn, ajint begin)
1012 {
1013 AjPStr key = NULL;
1014 AjPStr explain = NULL;
1015 AjPStr explainstr = NULL;
1016 const AjPStr error = NULL; /* value from table - do not delete */
1017 const AjPStr seqid = NULL;
1018 ajint i;
1019 const AjPStr size = NULL;
1020 const AjPStr seq = NULL;
1021 const AjPStr pos = NULL;
1022 const AjPStr tm = NULL;
1023 const AjPStr gc = NULL;
1024
1025 /*
1026 ** Typical primer3 output looks like:
1027 ** PRIMER_SEQUENCE_ID=em-id:HSFAU
1028 ** PRIMER_PAIR_PENALTY=0.1974
1029 ** PRIMER_LEFT_PENALTY=0.007073
1030 ** PRIMER_RIGHT_PENALTY=0.190341
1031 ** PRIMER_INTERNAL_OLIGO_PENALTY=0.132570
1032 ** PRIMER_LEFT_SEQUENCE=AGTCGCCAATATGCAGCTCT
1033 ** PRIMER_RIGHT_SEQUENCE=GGAGCACGACTTGATCTTCC
1034 ** PRIMER_INTERNAL_OLIGO_SEQUENCE=GGAGCTACACACCTTCGAGG
1035 ** PRIMER_LEFT=47,20
1036 ** PRIMER_RIGHT=183,20
1037 ** PRIMER_INTERNAL_OLIGO=80,20
1038 ** PRIMER_LEFT_TM=60.007
1039 ** PRIMER_RIGHT_TM=59.810
1040 ** PRIMER_INTERNAL_OLIGO_TM=59.867
1041 ** PRIMER_LEFT_GC_PERCENT=50.000
1042 ** PRIMER_RIGHT_GC_PERCENT=55.000
1043 ** PRIMER_INTERNAL_OLIGO_GC_PERCENT=60.000
1044 ** PRIMER_LEFT_SELF_ANY=4.00
1045 ** PRIMER_RIGHT_SELF_ANY=5.00
1046 ** PRIMER_INTERNAL_OLIGO_SELF_ANY=4.00
1047 ** PRIMER_LEFT_SELF_END=2.00
1048 ** PRIMER_RIGHT_SELF_END=1.00
1049 ** PRIMER_INTERNAL_OLIGO_SELF_END=4.00
1050 ** PRIMER_LEFT_END_STABILITY=7.9000
1051 ** PRIMER_RIGHT_END_STABILITY=8.2000
1052 ** PRIMER_PAIR_COMPL_ANY=5.00
1053 ** PRIMER_PAIR_COMPL_END=1.00
1054 ** PRIMER_PRODUCT_SIZE=137
1055 */
1056
1057 /* check for errors */
1058 ajStrAssignC(&key, "PRIMER_ERROR");
1059 error = ajTableFetchS(table, key);
1060 if(error != NULL)
1061 ajErr("%S", error);
1062
1063 ajStrAssignC(&key, "PRIMER_WARNING");
1064 error = ajTableFetchS(table, key);
1065 if(error != NULL)
1066 ajWarn("%S", error);
1067
1068
1069 /* get the sequence id */
1070 ajStrAssignC(&key, "PRIMER_SEQUENCE_ID");
1071 seqid = ajTableFetchS(table, key);
1072 ajFmtPrintF(outfile, "\n# EPRIMER3 RESULTS FOR %S\n\n", seqid);
1073
1074 /* get information on the analysis */
1075 ajStrAssignC(&key, "PRIMER_LEFT_EXPLAIN");
1076 explain = ajTableFetchmodS(table, key);
1077
1078 if(explain != NULL)
1079 {
1080 ajStrAssignS(&explainstr, explain);
1081 ajStrExchangeCC(&explainstr, ",", "\n#");
1082 ajFmtPrintF(outfile, "# FORWARD PRIMER STATISTICS:\n# %S\n\n",
1083 explainstr);
1084 }
1085 ajStrAssignC(&key, "PRIMER_RIGHT_EXPLAIN");
1086 explain = ajTableFetchmodS(table, key);
1087
1088 if(explain != NULL)
1089 {
1090 ajStrAssignS(&explainstr, explain);
1091 ajStrExchangeCC(&explainstr, ",", "\n#");
1092 ajFmtPrintF(outfile, "# REVERSE PRIMER STATISTICS:\n# %S\n\n",
1093 explainstr);
1094 }
1095 ajStrAssignC(&key, "PRIMER_PAIR_EXPLAIN");
1096 explain = ajTableFetchmodS(table, key);
1097
1098 if(explain != NULL)
1099 {
1100 ajStrAssignS(&explainstr, explain);
1101 ajStrExchangeCC(&explainstr, ",", "\n#");
1102 ajFmtPrintF(outfile, "# PRIMER PAIR STATISTICS:\n# %S\n\n",
1103 explainstr);
1104 }
1105 ajStrAssignC(&key, "PRIMER_INTERNAL_OLIGO_EXPLAIN");
1106 explain = ajTableFetchmodS(table, key);
1107
1108 if(explain != NULL)
1109 {
1110 ajStrAssignS(&explainstr, explain);
1111 ajStrExchangeCC(&explainstr, ",", "\n#");
1112 ajFmtPrintF(outfile, "# INTERNAL OLIGO STATISTICS:\n# %S\n\n",
1113 explainstr);
1114 }
1115
1116 /* table header */
1117 ajFmtPrintF(outfile,"# Start Len Tm "
1118 "GC%% Sequence\n\n");
1119
1120 /* get the results */
1121 for(i=0; i <= numreturn; i++)
1122 {
1123 /* product data */
1124 size = eprimer3_tableget("PRIMER_PRODUCT_SIZE", i, "", table);
1125 if(size != NULL)
1126 ajFmtPrintF(outfile, "%4d PRODUCT SIZE: %S\n",
1127 i+1, size);
1128
1129 /* left primer data */
1130 tm = eprimer3_tableget("PRIMER_LEFT", i, "TM", table);
1131 gc = eprimer3_tableget("PRIMER_LEFT", i, "GC_PERCENT", table);
1132 pos = eprimer3_tableget("PRIMER_LEFT", i, "", table);
1133 seq = eprimer3_tableget("PRIMER_LEFT", i, "SEQUENCE", table);
1134 eprimer3_write_primer(outfile, "FORWARD PRIMER",
1135 pos, tm, gc, seq, AJFALSE, begin);
1136
1137
1138 /* right primer data */
1139 tm = eprimer3_tableget("PRIMER_RIGHT", i, "TM", table);
1140 gc = eprimer3_tableget("PRIMER_RIGHT", i, "GC_PERCENT", table);
1141 pos = eprimer3_tableget("PRIMER_RIGHT", i, "", table);
1142 seq = eprimer3_tableget("PRIMER_RIGHT", i, "SEQUENCE", table);
1143 eprimer3_write_primer(outfile, "REVERSE PRIMER",
1144 pos, tm, gc, seq, AJTRUE, begin);
1145
1146
1147 /* internal oligo data */
1148
1149 tm = eprimer3_tableget("PRIMER_INTERNAL_OLIGO", i, "TM", table);
1150 gc = eprimer3_tableget("PRIMER_INTERNAL_OLIGO", i, "GC_PERCENT",
1151 table);
1152 pos = eprimer3_tableget("PRIMER_INTERNAL_OLIGO", i, "", table);
1153 seq = eprimer3_tableget("PRIMER_INTERNAL_OLIGO", i, "SEQUENCE", table);
1154 eprimer3_write_primer(outfile, "INTERNAL OLIGO",
1155 pos, tm, gc, seq, AJFALSE, begin);
1156
1157 ajFmtPrintF(outfile, "\n");
1158
1159 }
1160
1161
1162 ajStrDel(&key);
1163 ajStrDel(&explainstr);
1164
1165 return;
1166 }
1167
1168
1169
1170
1171 /* @funcstatic eprimer3_tableget **********************************************
1172 **
1173 ** Read the results out of the tag/value table
1174 **
1175 ** @param [r] key1 [const char *] First half of table key base string
1176 ** @param [r] number [ajint] Table key numeric part
1177 ** @param [r] key2 [const char *] Second half of table key base string
1178 ** (minus '_')
1179 ** @param [r] table [const AjPTable] Table of tag/value result pairs
1180 ** @return [const AjPStr] Table value
1181 **
1182 ******************************************************************************/
1183
eprimer3_tableget(const char * key1,ajint number,const char * key2,const AjPTable table)1184 static const AjPStr eprimer3_tableget(const char *key1, ajint number,
1185 const char *key2,
1186 const AjPTable table)
1187 {
1188 AjPStr fullkey = NULL;
1189 AjPStr keynum = NULL;
1190 const AjPStr value = NULL;
1191
1192 ajStrAssignC(&fullkey, key1);
1193
1194 if(number > 0)
1195 {
1196 ajStrAppendC(&fullkey, "_");
1197 ajStrFromInt(&keynum, number);
1198 ajStrAppendS(&fullkey, keynum);
1199 }
1200
1201 if(strcmp(key2, ""))
1202 ajStrAppendC(&fullkey, "_");
1203
1204 ajStrAppendC(&fullkey, key2);
1205 ajDebug("Constructed key=%S\n", fullkey);
1206 value = ajTableFetchS(table, fullkey);
1207
1208
1209 ajStrDel(&fullkey);
1210 ajStrDel(&keynum);
1211
1212 return value;
1213 }
1214
1215
1216
1217
1218 /* @funcstatic eprimer3_write_primer ******************************************
1219 **
1220 ** Write out one primer or oligo line to the output file
1221 **
1222 ** @param [u] outfile [AjPFile] Report outfile
1223 ** @param [r] tag [const char *] Tag on output line
1224 ** @param [r] pos [const AjPStr] Start and length string
1225 ** @param [r] tm [const AjPStr] Tm of primer
1226 ** @param [r] gc [const AjPStr] GC% of primer
1227 ** @param [r] seq [const AjPStr] Sequence of primer
1228 ** @param [r] rev [AjBool] Sequence is the reverse-complement primer
1229 ** @param [r] begin [ajint] Start position of subsequence (-sbegin)
1230 ** @return [void]
1231 **
1232 ******************************************************************************/
1233
eprimer3_write_primer(AjPFile outfile,const char * tag,const AjPStr pos,const AjPStr tm,const AjPStr gc,const AjPStr seq,AjBool rev,ajint begin)1234 static void eprimer3_write_primer(AjPFile outfile, const char *tag,
1235 const AjPStr pos,
1236 const AjPStr tm, const AjPStr gc,
1237 const AjPStr seq,
1238 AjBool rev, ajint begin)
1239 {
1240 ajint startint;
1241 ajint lenint;
1242 float tmfloat;
1243 float gcfloat;
1244 AjPStr start = NULL;
1245 AjPStr lenstr = NULL;
1246 ajlong comma;
1247
1248 if(ajStrGetLen(pos))
1249 {
1250 ajStrToFloat(tm, &tmfloat);
1251 ajStrToFloat(gc, &gcfloat);
1252 comma = ajStrFindC(pos, ",");
1253 ajStrAssignS(&start, pos);
1254 ajStrCutRange(&start, comma, ajStrGetLen(start)-1);
1255 ajStrToInt(start, &startint);
1256 startint += begin;
1257 ajStrAssignS(&lenstr, pos);
1258 ajStrCutRange(&lenstr, 0, comma);
1259 ajStrToInt(lenstr, &lenint);
1260 if(rev)
1261 startint = startint - lenint + 1;
1262
1263 ajFmtPrintF(outfile, " %s %6d %4d %2.2f %2.2f %S\n\n",
1264 tag, startint, lenint, tmfloat, gcfloat, seq);
1265 }
1266
1267
1268 ajStrDel(&start);
1269 ajStrDel(&lenstr);
1270
1271 return;
1272 }
1273