1 #include <U2Core/disable-warnings.h>
2 U2_DISABLE_WARNINGS
3
4 /*
5 Copyright (c) 1996,1997,1998,1999,2000,2001,2004,2006
6 Whitehead Institute for Biomedical Research, Steve Rozen
7 (http://jura.wi.mit.edu/rozen), and Helen Skaletsky
8 All rights reserved.
9
10 Redistribution and use in source and binary forms, with or without
11 modification, are permitted provided that the following conditions are
12 met:
13
14 * Redistributions of source code must retain the above copyright
15 notice, this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above
17 copyright notice, this list of conditions and the following disclaimer
18 in the documentation and/or other materials provided with the
19 distribution.
20 * Neither the names of the copyright holders nor contributors may
21 be used to endorse or promote products derived from this software
22 without specific prior written permission.
23
24 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
28 OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
29 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
30 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
31 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
32 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
33 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35 */
36
37 #include <errno.h>
38 #include <stdio.h>
39 #include <math.h>
40 #include <signal.h>
41 //#include <unistd.h>
42 #include <float.h>
43 #include "primer3.h"
44 #include "primer3_main.h"
45 #include "boulder_input.h" //for free_seq_lib only
46 #include "oligotm.h"
47
48 /* #define's */
49
50 /*
51 * Panic messages for when the program runs out of memory. pr_program_name and
52 * pr_program_name_len must be set at the beginning of main.
53 */
54 #define OOM_MESSAGE ": out of memory\n"
55 #define OOM_MESSAGE_LEN 16
56 //#define OOM_STMT1 write(2, pr_program_name, pr_program_name_len)
57 #define OOM_STMT2 write(2, OOM_MESSAGE, OOM_MESSAGE_LEN), exit(-2)
58 //#define OOM_ERROR OOM_STMT1, OOM_STMT2
59 #define OOM_ERROR {int * p = NULL; *p = 1;}
60
61 #ifndef MAX_PRIMER_LENGTH
62 #error "Define MAX_PRIMER_LENGTH in Makefile..."
63 /* to ensure that MAX_PRIMER_LENGTH <= DPAL_MAX_ALIGN. */
64 #endif
65 #if (MAX_PRIMER_LENGTH > DPAL_MAX_ALIGN)
66 #error "MAX_PRIMER_LENGTH must be <= DPAL_MAX_ALIGN"
67 #endif
68
69 #define MAX_NN_TM_LENGTH 36 /* The maxium length for which to use the
70 nearest neighbor model when calculating
71 oligo Tms. */
72
73 #define ALIGN_SCORE_UNDEF SHRT_MIN
74
75 #define MACRO_CAT_2(A,B) A##B
76
77 #define INITIAL_LIST_LEN 2000 /* Initial size of oligo lists. */
78 #define INITIAL_NUM_RETURN 5 /* Initial space to allocate for pairs to
79 return. */
80
81 #define PAIR_OK 1
82 #define PAIR_FAILED 0
83
84 #define OK_OR_MUST_USE(H) ((H)->ok == OV_OK || (H)->must_use)
85
86 /* Function declarations. */
87 static void add_must_use_warnings(seq_args *, const char *,
88 const oligo_stats *);
89 static void add_pair(const primer_pair *, pair_array_t *);
90 static short align(const char *, const char*, const dpal_args *a);
91 static int choose_pair(const primer_args *,
92 seq_args *,
93 int, int, int, int,
94 const dpal_args*, const dpal_args*,
95 const dpal_args*,
96 pair_array_t *,
97 int * cancelFlag, int * progress, Primer3Context * ctx);
98 static void check_sequence_quality(const primer_args *, primer_rec *,
99 oligo_type, const seq_args *, int, int,
100 int *, int *);
101 static int choose_internal_oligo(const primer_rec *, const primer_rec *,
102 int, int *, seq_args *,
103 const primer_args *,
104 const dpal_args*,
105 const dpal_args *, const dpal_args *, Primer3Context * ctx);
106 void compute_position_penalty(const primer_args *, const seq_args *,
107 primer_rec *, oligo_type);
108
109 static void create_and_print_file(const seq_args *, int, const primer_rec[],
110 const oligo_type, const int, const int,
111 const char *);
112 static int find_stop_codon(const char *, int, int);
113 static void gc_and_n_content(const int, const int, const char *, primer_rec *);
114 static int make_internal_oligos_list(const primer_args *,
115 seq_args *,
116 int *, const dpal_args *,
117 const dpal_args*,
118 const dpal_args *, Primer3Context * ctx);
119 static int make_lists(const primer_args *,
120 seq_args *,
121 int *, int *, const dpal_args *,
122 const dpal_args*, const dpal_args *, int * cancelFlag, int * progress, Primer3Context * ctx);
123 static double obj_fn(const primer_args *, primer_pair *);
124 static short oligo_max_template_mispriming(const primer_rec *);
125 static int oligo_overlaps_interval(const int, const int,
126 interval_array_t, const int);
127 static int oligo_pair_seen(const primer_pair *, const pair_array_t *);
128 static void oligo_param(const primer_args *pa,
129 primer_rec *, oligo_type,
130 const dpal_args*,
131 const dpal_args *, const dpal_args *,
132 seq_args *, oligo_stats *, Primer3Context * ctx);
133 static int pair_param(const primer_args *, seq_args *,
134 int, int, int, primer_pair *,
135 const dpal_args*,
136 const dpal_args *, const dpal_args*, Primer3Context * ctx);
137 static int pair_spans_target(const primer_pair *, const seq_args *);
138 static int primer_pair_comp(const void *, const void*);
139 static int primer_rec_comp(const void *, const void *);
140 static void print_list(const seq_args *, const primer_args *,
141 int, int, int);
142 static void print_list_header(FILE *, oligo_type, int, int);
143 static void print_oligo(FILE *, const seq_args *, int, const primer_rec *,
144 oligo_type, int, int);
145 //static void print_usage();
146 // static void boulder_print_pairs(const program_args *, const primer_args *,
147 // const seq_args *, const pair_array_t *);
148 static void print_explain(const oligo_stats *, oligo_type);
149 static void print_all_explain(const primer_args *, const seq_args *);
150 static FILE *safe_fopen(const char*, const char*);
151 static void set_dpal_args(dpal_args *);
152 static void sig_handler(int);
153 static double p_obj_fn(const primer_args *, primer_rec *, int );
154 static void oligo_compl(primer_rec *, const primer_args *, seq_args *,
155 oligo_type, const dpal_args *,
156 const dpal_args*, const dpal_args *);
157 static void oligo_mispriming(primer_rec *,
158 const primer_args *,
159 seq_args *,
160 oligo_type,
161 const dpal_args *, Primer3Context * ctx);
162 static int pair_repeat_sim(primer_pair *, const primer_args *);
163 static void boulder_print_oligos(const primer_args *,
164 const seq_args *, int, oligo_type);
165 static void free_repeat_sim_score(int, int, int, Primer3Context * ctx);
166
167 /* edited by T. Koressaar for lowercase masking: */
168 static void check_if_lowercase_masked(const int position,
169 const char *sequence,
170 primer_rec *h);
171
172 // /* Other global variables. */
173 // const char *pr_program_name;
174 // int pr_program_name_len;
175
176 /* Primer lists and their lengths. */
177 // static primer_rec *f;
178 // static primer_rec *r;
179 // static primer_rec *mid;
180 // static int f_len, r_len, mid_len;
181
182 /* Argument structures for dpal (too cumbersome to pass around
183 individually.) */
184 //static dpal_args *lib_local_dpal_args, *lib_local_end_dpal_args;
185
186 typedef struct oligo_array {
187 int len;
188 primer_rec *data;
189 } oligo_array;
190
191 #define FREE_STUFF { \
192 free(sa); free(local_args); free(end_args); \
193 free(local_end_args); \
194 free((ctx->lib_local_dpal_args)); free((ctx->lib_local_end_dpal_args)); \
195 if (0 != best_pairs.storage_size) free(best_pairs.pairs); \
196 free_seq_lib(&pa->repeat_lib); \
197 free_seq_lib(&pa->io_mishyb_library); free(pa); \
198 }
199
200 //int
201 //primer3_main(argc,argv)
202 // int argc;
203 // char *argv[];
204 //{
205 // program_args prog_args;
206 // primer_args *pa;
207 // seq_args *sa;
208 // dpal_args *local_args;
209 // dpal_args *end_args;
210 // dpal_args *local_end_args;
211 // int input_found=0;
212 // pair_array_t best_pairs;
213 // int n_f, n_r, n_m; /* Will be initialized in pr_choice. */
214 //
215 // pr_program_name = argv[0];
216 // pr_program_name_len = strlen(argv[0]);
217 // /*
218 // * We set up some signal handlers in case someone starts up the program
219 // * from the command line, wonders why nothing is happening, and then kills
220 // * the program.
221 // */
222 // signal(SIGINT, sig_handler);
223 // signal(SIGTERM, sig_handler);
224 //
225 // /*
226 // * We allocate the following structures on the heap rather than on the
227 // * stack in order to take advantage of testcenter memory access checking
228 // * during testing.
229 // */
230 // pa = pr_safe_malloc(sizeof(*pa));
231 // sa = pr_safe_malloc(sizeof(*sa));
232 // local_args = pr_safe_malloc(sizeof(*local_args));
233 // end_args = pr_safe_malloc(sizeof(*end_args));
234 // local_end_args = pr_safe_malloc(sizeof(*local_end_args));
235 // lib_local_end_dpal_args = pr_safe_malloc(sizeof(*lib_local_end_dpal_args));
236 // lib_local_dpal_args = pr_safe_malloc(sizeof(*lib_local_dpal_args));
237 //
238 // best_pairs.storage_size = best_pairs.num_pairs = 0;
239 // pr_set_default_global_args(pa);
240 //
241 // memset(&prog_args, 0, sizeof(prog_args));
242 // while (--argc > 0) {
243 // argv++;
244 // if (!strcmp(*argv, "-format_output"))
245 // prog_args.format_output = 1;
246 // else if (!strcmp(*argv, "-2x_compat"))
247 // prog_args.twox_compat = 1;
248 // else if (!strcmp(*argv, "-strict_tags"))
249 // prog_args.strict_tags = 1;
250 // else {
251 // print_usage();
252 // FREE_STUFF;
253 // exit(-1);
254 // }
255 // }
256 //
257 // set_dpal_args(local_args);
258 // local_args->flag = DPAL_LOCAL;
259 //
260 // set_dpal_args(end_args);
261 // end_args->flag = DPAL_GLOBAL_END;
262 //
263 // set_dpal_args(local_end_args);
264 // local_end_args->flag = DPAL_LOCAL_END;
265 //
266 // /*
267 // * Read the data from input stream record by record and process it if
268 // * there are no errors.
269 // */
270 // while ((read_record(&prog_args, pa, sa)) > 0) {
271 // input_found = 1;
272 //
273 // *lib_local_dpal_args = *local_args;
274 // if (pa->lib_ambiguity_codes_consensus) {
275 // PR_ASSERT(dpal_set_ambiguity_code_matrix(lib_local_dpal_args));
276 // }
277 //
278 // *lib_local_end_dpal_args = *local_end_args;
279 // if (pa->lib_ambiguity_codes_consensus) {
280 // PR_ASSERT(dpal_set_ambiguity_code_matrix(lib_local_end_dpal_args));
281 // }
282 //
283 // n_f = n_m = n_r = 0;
284 // if (NULL == sa->error.data && NULL == pa->glob_err.data) {
285 // pr_choice(pa, sa, local_args, end_args, local_end_args,
286 // &best_pairs, &n_f, &n_r, &n_m);
287 // }
288 // if (NULL != pa->glob_err.data)
289 // pr_append_new_chunk(&sa->error, pa->glob_err.data);
290 //
291 // if (pick_pcr_primers == pa->primer_task
292 // || pick_pcr_primers_and_hyb_probe == pa->primer_task) {
293 // if (prog_args.format_output) {
294 // format_pairs(stdout, pa, sa, &best_pairs);
295 // }
296 // else {
297 // boulder_print_pairs(&prog_args, pa, sa, &best_pairs);
298 // }
299 // }
300 // else if(pa->primer_task == pick_left_only) {
301 // if (prog_args.format_output)
302 // format_oligos(stdout, pa, sa, f, n_f, OT_LEFT);
303 // else boulder_print_oligos(pa, sa, n_f, OT_LEFT);
304 // }
305 // else if(pa->primer_task == pick_right_only) {
306 // if (prog_args.format_output)
307 // format_oligos(stdout, pa, sa, r, n_r, OT_RIGHT);
308 // else boulder_print_oligos(pa, sa, n_r, OT_RIGHT);
309 // }
310 // else if(pa->primer_task == pick_hyb_probe_only) {
311 // if(prog_args.format_output)
312 // format_oligos(stdout, pa, sa, mid, n_m, OT_INTL);
313 // else boulder_print_oligos(pa, sa, n_m, OT_INTL);
314 // }
315 //
316 // best_pairs.num_pairs = 0;
317 //
318 // if(pa->repeat_lib.seq_num > 0 || pa->io_mishyb_library.seq_num > 0)
319 // free_repeat_sim_score(n_f, n_r, n_m);
320 // if (NULL != sa->internal_input) free(sa->internal_input);
321 // if (NULL != sa->left_input) free(sa->left_input);
322 // if (NULL != sa->right_input) free(sa->right_input);
323 // if (NULL != sa->sequence) free(sa->sequence);
324 // if (NULL != sa->quality) free(sa->quality);
325 // if (NULL != sa->trimmed_seq) free(sa->trimmed_seq);
326 //
327 // /* edited by T. Koressaar for lowercase masking */
328 // if (NULL != sa->trimmed_orig_seq) free(sa->trimmed_orig_seq);
329 //
330 // if (NULL != sa->upcased_seq) free(sa->upcased_seq);
331 // if (NULL != sa->upcased_seq_r) free(sa->upcased_seq_r);
332 // if (NULL != sa->sequence_name) free(sa->sequence_name);
333 // if (NULL != sa->error.data) free(sa->error.data);
334 // if (NULL != sa->warning.data) free(sa->warning.data);
335 //
336 // if (NULL != pa->glob_err.data) {
337 // fprintf(stderr, "%s: %s\n", pr_program_name, pa->glob_err.data);
338 // free(pa->glob_err.data);
339 // FREE_STUFF;
340 // exit(-4);
341 // }
342 // }
343 // FREE_STUFF; /* Free memory */
344 // if (0 == input_found) {
345 // print_usage();
346 // exit(-3);
347 // }
348 // return 0;
349 //}
350 //#undef FREE_STUFF
351
352 /*
353 * Find up to pa->num_return primer pairs for the sequence seq with t targets.
354 * Set sa->error and return 1 on error; otherwise return 0.
355 */
356 void
pr_choice(pa,sa,local_args,end_args,local_end_args,best_pairs,n_f,n_r,n_m,cancelFlag,progress,ctx)357 pr_choice(pa, sa, local_args, end_args, local_end_args, best_pairs,
358 n_f, n_r, n_m, cancelFlag, progress, ctx)
359 primer_args *pa;
360 seq_args *sa;
361 const dpal_args *local_args, *end_args, *local_end_args;
362 pair_array_t *best_pairs;
363 int *n_f; /* Number of acceptable left primers found. */
364 int *n_r; /* Number of acceptable right primers found. */
365 int *n_m; /* Number of acceptable internal oligos found. */
366 int * cancelFlag;
367 int * progress;
368 Primer3Context * ctx;
369 /* *n_f, *n_r, *n_m must initialized before calling this function. */
370 {
371 int i; /* Loop index. */
372 int int_num; /* Product size range counter. */
373 pair_array_t p;
374
375 PR_ASSERT(NULL != pa);
376 PR_ASSERT(NULL != sa);
377 PR_ASSERT(NULL != local_args);
378 PR_ASSERT(NULL != end_args);
379 PR_ASSERT(NULL != local_end_args);
380 PR_ASSERT(NULL != n_f);
381 PR_ASSERT(NULL != n_r);
382 PR_ASSERT(NULL != n_m);
383 PR_ASSERT(0 == *n_f);
384 PR_ASSERT(0 == *n_r);
385 PR_ASSERT(0 == *n_m);
386
387 if (_pr_data_control(pa, sa) !=0 ) return;
388
389 if (NULL == (ctx->f)) {
390 (ctx->f) = pr_safe_malloc(sizeof(*(ctx->f)) * INITIAL_LIST_LEN);
391 (ctx->r) = pr_safe_malloc(sizeof(*(ctx->r)) * INITIAL_LIST_LEN);
392 (ctx->f_len) = (ctx->r_len) = INITIAL_LIST_LEN;
393 }
394
395 if (make_lists(pa, sa, n_f, n_r,local_args,
396 end_args, local_end_args, cancelFlag, progress, ctx)!=0 ) {
397 return;
398 }
399
400 if ((pa->primer_task == pick_hyb_probe_only
401 || pa->primer_task == pick_pcr_primers_and_hyb_probe)
402 && make_internal_oligos_list(pa, sa, n_m, local_args,
403 end_args, local_end_args, ctx) != 0)
404 return;
405 /* Creates files with left, right, and internal oligos. */
406 // if (pa->file_flag) print_list(sa, pa, *n_f, *n_r, *n_m);
407 /* We sort _after_ printing lists to maintain the order of test output. */
408 if(pa->primer_task != pick_left_only
409 && pa->primer_task != pick_hyb_probe_only)
410 qsort(&(ctx->r)[0], *n_r, sizeof(*(ctx->r)), primer_rec_comp);
411 if(pa->primer_task != pick_right_only
412 && pa->primer_task != pick_hyb_probe_only)
413 qsort(&(ctx->f)[0], *n_f, sizeof(*(ctx->f)), primer_rec_comp);
414 if(pa->primer_task == pick_hyb_probe_only)
415 qsort(&(ctx->mid)[0], *n_m, sizeof(*(ctx->mid)), primer_rec_comp);
416 p.storage_size = p.num_pairs = 0;
417 if (pa->primer_task == pick_pcr_primers
418 || pa->primer_task == pick_pcr_primers_and_hyb_probe){
419
420 /* Look for pa->num_return best primer pairs. */
421 for(int_num=0; int_num < pa->num_intervals; int_num++) {
422 if( *cancelFlag ) {
423 break;
424 }
425 if(choose_pair(pa, sa, *n_f, *n_r, *n_m, int_num,
426 local_args, end_args, local_end_args, &p, cancelFlag, progress, ctx)!=0)
427 continue;
428
429 for (i = 0; i < p.num_pairs && best_pairs->num_pairs < pa->num_return;
430 i++) {
431 if (!oligo_pair_seen(&p.pairs[i], best_pairs))
432 add_pair(&p.pairs[i], best_pairs);
433 }
434
435 if (pa->num_return == best_pairs->num_pairs) break;
436 p.num_pairs = 0;
437 }
438 }
439
440 /* If it was necessary to use a left_input, right_input,
441 or internal_oligo_input primer that was
442 unacceptable, then add warnings. */
443 if (pa->pick_anyway) {
444 if (sa->left_input) {
445 add_must_use_warnings(sa, "Left primer", &sa->left_expl);
446 }
447 if (sa->right_input) {
448 add_must_use_warnings(sa, "Right primer", &sa->right_expl);
449 }
450 if (sa->internal_input) {
451 add_must_use_warnings(sa, "Hybridization probe", &sa->intl_expl);
452 }
453 }
454
455 if (0 != p.storage_size) free(p.pairs);
456 return;
457 }
458
459 /* Call this function only if the 'stat's contains
460 the _errors_ associated with a given primer
461 i.e. that primer was supplied by the caller
462 and pick_anyway is set. */
463 static void
add_must_use_warnings(sa,text,stats)464 add_must_use_warnings(sa, text, stats)
465 seq_args *sa;
466 const char* text;
467 const oligo_stats *stats;
468 {
469 const char *sep = "/";
470 pr_append_str s;
471
472 s.data = NULL;
473 s.storage_size = 0;
474
475 if (stats->ns) pr_append_w_sep(&s, sep, "Too many Ns");
476 if (stats->target) pr_append_w_sep(&s, sep, "Overlaps Target");
477 if (stats->excluded) pr_append_w_sep(&s, sep, "Overlaps Excluded Region");
478 if (stats->gc) pr_append_w_sep(&s, sep, "Unacceptable GC content");
479 if (stats->gc_clamp) pr_append_w_sep(&s, sep, "No GC clamp");
480 if (stats->temp_min) pr_append_w_sep(&s, sep, "Tm too low");
481 if (stats->temp_max) pr_append_w_sep(&s, sep, "Tm too high");
482 if (stats->compl_any) pr_append_w_sep(&s, sep, "High self complementarity");
483 if (stats->compl_end)
484 pr_append_w_sep(&s, sep, "High end self complementarity");
485 if (stats->repeat_score)
486 pr_append_w_sep(&s, sep, "High similarity to mispriming or mishyb library");
487 if (stats->poly_x) pr_append_w_sep(&s, sep, "Long poly-X");
488 if (stats->seq_quality) pr_append_w_sep(&s, sep, "Low sequence quality");
489 if (stats->stability) pr_append_w_sep(&s, sep, "High 3' stability");
490 if (stats->no_orf) pr_append_w_sep(&s, sep, "Would not amplify any ORF");
491
492 /* edited by T. Koressaar for lowercase masking: */
493 if (stats->gmasked)
494 pr_append_w_sep(&s, sep, "Masked with lowercase letter");
495
496 if (s.data) {
497 pr_append_new_chunk(&sa->warning, text);
498 pr_append(&sa->warning, " is unacceptable: ");
499 pr_append(&sa->warning, s.data);
500 free(s.data);
501 }
502
503 }
504
505 /* Return 1 iff pair is already in the first num_pairs elements of retpair. */
506 static int
oligo_pair_seen(pair,retpair)507 oligo_pair_seen(pair, retpair)
508 const primer_pair *pair;
509 const pair_array_t *retpair;
510 {
511 const primer_pair *q, *stop;
512
513 /* OLD CODE: const primer_pair *q = &retpair->pairs[0],
514 *stop = &retpair->pairs[retpair->num_pairs]; */
515
516 /* retpair might not have any pairs in it yet (add_pair
517 allocates memory for retpair->pairs. */
518 if (retpair->num_pairs == 0)
519 return 0;
520
521 q = &retpair->pairs[0];
522 stop = &retpair->pairs[retpair->num_pairs];
523
524 for (; q < stop; q++) {
525 if (q->left->start == pair->left->start
526 && q->left->length == pair->left->length
527 && q->right->start == pair->right->start
528 && q->right->length == pair->right->length) return 1;
529 }
530 return 0;
531 }
532
533 /* Add 'pair' to 'retpair'. */
534 static void
add_pair(pair,retpair)535 add_pair(pair, retpair)
536 const primer_pair *pair;
537 pair_array_t *retpair;
538 {
539 if (0 == retpair->storage_size) {
540 retpair->storage_size = INITIAL_NUM_RETURN;
541 retpair->pairs
542 = pr_safe_malloc(retpair->storage_size * sizeof(*retpair->pairs));
543 } else if (retpair->storage_size == retpair->num_pairs) {
544 retpair->storage_size *= 2;
545 retpair->pairs
546 = pr_safe_realloc(retpair->pairs,
547 retpair->storage_size * sizeof(*retpair->pairs));
548 }
549 retpair->pairs[retpair->num_pairs] = *pair;
550 retpair->num_pairs++;
551 }
552
553
554 /* Make lists of acceptable left and right primers. */
555 static int
make_lists(pa,sa,n_f,n_r,local_args,end_args,local_end_args,cancelFlag,progress,ctx)556 make_lists(pa, sa, n_f, n_r, local_args, end_args, local_end_args, cancelFlag, progress, ctx)
557 const primer_args *pa;
558 seq_args *sa;
559 int *n_f, *n_r;
560 const dpal_args *local_args, *end_args, *local_end_args;
561 int * cancelFlag;
562 int * progress;
563 Primer3Context * ctx;
564 {
565 int left, right;
566 int i,j,n,k,pr_min;
567 int tar_l, tar_r, f_b, r_b;
568
569 /*
570 * The position of the intial base of the rightmost stop codon that is
571 * to the left of sa->start_codon_pos; valid only if sa->start_codon_pos
572 * is "not null". We will not want to include a stop codon to the right
573 * of the the start codon in the amplicon.
574 */
575 int stop_codon1 = -1;
576
577 char s[MAX_PRIMER_LENGTH+1],s1[MAX_PRIMER_LENGTH+1];
578 primer_rec h;
579
580 left = right = 0;
581 if (!PR_START_CODON_POS_IS_NULL(sa)) {
582 stop_codon1 = find_stop_codon(sa->trimmed_seq,
583 sa->start_codon_pos, -1);
584
585 sa->stop_codon_pos = find_stop_codon(sa->trimmed_seq,
586 sa->start_codon_pos, 1);
587 sa->stop_codon_pos += sa->incl_s;
588 }
589
590 pr_min = INT_MAX;
591 for(i=0;i<pa->num_intervals;i++)if(pa->pr_min[i]<pr_min)
592 pr_min= pa->pr_min[i];
593
594 PR_ASSERT(INT_MAX > (n=strlen(sa->trimmed_seq)));
595 tar_r = 0;
596 tar_l = n;
597 for(i=0;i<sa->num_targets;i++) {
598 if(sa->tar[i][0]>tar_r)tar_r = sa->tar[i][0];
599 if(sa->tar[i][0]+sa->tar[i][1]-1<tar_l)tar_l=
600 sa->tar[i][0]+sa->tar[i][1]-1;
601 }
602
603 if (_PR_DEFAULT_POSITION_PENALTIES(pa)) {
604 if (0 == tar_r) tar_r = n;
605 if (tar_l == n) tar_l = 0;
606 } else {
607 tar_r = n;
608 tar_l = 0;
609 }
610
611 if (pa->primer_task == pick_left_only)
612 f_b = n - 1;
613 else if (tar_r - 1 < n - pr_min + pa->primer_max_size - 1
614 && !(pa->pick_anyway && sa->left_input))
615 f_b=tar_r - 1;
616 else
617 f_b = n - pr_min + pa->primer_max_size-1;
618 k = 0;
619 if(pa->primer_task != pick_right_only && pa->primer_task != pick_hyb_probe_only){
620 left=n; right=0;
621 for (i = f_b; i >= pa->primer_min_size - 1; i--) {
622 if(*cancelFlag) {
623 return -1;
624 }
625 *progress = 25*(f_b - i)/(f_b - (pa->primer_min_size - 2));
626 s[0]='\0';
627 for (j = pa->primer_min_size; j <= pa->primer_max_size; j++) {
628 if (i-j > n-pr_min-1 && pick_left_only != pa->primer_task) continue;
629 if (i-j+1>=0) {
630 if (k >= (ctx->f_len)) {
631 (ctx->f_len) += ((ctx->f_len) >> 1);
632 (ctx->f) = pr_safe_realloc((ctx->f), (ctx->f_len) * sizeof(*(ctx->f)));
633 }
634 h.start=i-j+1;
635 h.length=j;
636 h.repeat_sim.score = NULL;
637 _pr_substr(sa->trimmed_seq,h.start,h.length,s);
638
639 if (sa->left_input && strcmp_nocase(sa->left_input, s))
640 continue;
641
642 h.must_use = (sa->left_input && pa->pick_anyway);
643
644 if (pa->explain_flag) sa->left_expl.considered++;
645
646 if (!PR_START_CODON_POS_IS_NULL(sa)
647 /* Make sure the primer would amplify at least part of
648 the ORF. */
649 && (0 != (h.start - sa->start_codon_pos) % 3
650 || h.start <= stop_codon1
651 || (sa->stop_codon_pos != -1
652 && h.start >= sa->stop_codon_pos))) {
653 if (pa->explain_flag) sa->left_expl.no_orf++;
654 if (!pa->pick_anyway) continue;
655 }
656
657 h.repeat_sim.score = NULL;
658 oligo_param(pa, &h, OT_LEFT, local_args, end_args,
659 local_end_args, sa, &sa->left_expl, ctx);
660 if (OK_OR_MUST_USE(&h)) {
661 h.quality = p_obj_fn(pa, &h, 0);
662 (ctx->f)[k] = h;
663 if((ctx->f)[k].start < left) left=(ctx->f)[k].start;
664 k++;
665 } else if (h.ok==OV_TOO_MANY_NS || h.ok==OV_INTERSECT_TARGET
666 || h.ok==OV_SELF_ANY || h.ok==OV_END_STAB
667 || h.ok==OV_POLY_X || h.ok==OV_EXCL_REGION || h.ok==OV_GC_CLAMP
668 || h.ok == OV_SEQ_QUALITY || h.ok == OV_LIB_SIM ) {
669 /* Break from the inner for loop, because there is no
670 legal longer oligo with the same 3' sequence. */
671 break;
672 }
673 }
674 else break;
675 }
676 }
677 }
678 *n_f = k;
679
680
681 if (pa->primer_task == pick_right_only)
682 r_b = 0;
683 else if (tar_l+1>pr_min - pa->primer_max_size
684 && !(pa->pick_anyway && sa->right_input))
685 r_b = tar_l+1;
686 else
687 r_b = pr_min - pa->primer_max_size;
688 k = 0;
689 if(pa->primer_task != pick_left_only && pa->primer_task != pick_hyb_probe_only) {
690 for(i=r_b; i<=n-pa->primer_min_size; i++) {
691 if(*cancelFlag) {
692 return -1;
693 }
694 *progress = 25 + 25*(i - r_b)/(n - pa->primer_min_size - r_b + 1);
695 s[0]='\0';
696 for(j = pa->primer_min_size; j <= pa->primer_max_size; j++) {
697 if (i+j<pr_min && pa->primer_task != pick_right_only) continue;
698 if(i+j-1<n) {
699 if (k >= (ctx->r_len)) {
700 (ctx->r_len) += ((ctx->r_len) >> 1);
701 (ctx->r) = pr_safe_realloc((ctx->r), (ctx->r_len) * sizeof(*(ctx->r)));
702 }
703 h.start=i+j-1;
704 h.length=j;
705 h.repeat_sim.score = NULL;
706 _pr_substr(sa->trimmed_seq,i,j,s);
707 _pr_reverse_complement(s,s1);
708
709 if (sa->right_input && strcmp_nocase(sa->right_input, s1))
710 continue;
711 h.must_use = (sa->right_input && pa->pick_anyway);
712
713 h.repeat_sim.score = NULL;
714 oligo_param(pa, &h, OT_RIGHT, local_args, end_args,
715 local_end_args, sa, &sa->right_expl, ctx);
716 sa->right_expl.considered++;
717 if (OK_OR_MUST_USE(&h)) {
718 h.quality = p_obj_fn(pa, &h, 0);
719 (ctx->r)[k] = h;
720 if ((ctx->r)[k].start > right) right=(ctx->r)[k].start;
721 k++;
722 } else if (h.ok==OV_TOO_MANY_NS || h.ok==OV_INTERSECT_TARGET
723 || h.ok==OV_SELF_ANY || h.ok==OV_END_STAB
724 || h.ok==OV_POLY_X || h.ok==OV_EXCL_REGION || h.ok==OV_GC_CLAMP
725 || h.ok == OV_SEQ_QUALITY || h.ok == OV_LIB_SIM ) {
726 /* Break from the inner for loop, because there is no
727 legal longer oligo with the same 3' sequence. */
728 break;
729 }
730 }
731 else break;
732 }
733 }
734 }
735 *n_r=k;
736
737 /*
738 * Return 1 if one of lists is empty or if leftmost left primer and
739 * rightmost right primer do not provide sufficient product size.
740 */
741 sa->left_expl.ok = *n_f;
742 sa->right_expl.ok = *n_r;
743 if ((pa->primer_task != pick_right_only &&
744 pa->primer_task != pick_hyb_probe_only && 0 == *n_f) ||
745 ((pa->primer_task != pick_left_only &&
746 pa->primer_task != pick_hyb_probe_only) && 0 == *n_r)) return 1;
747 else if ((pa->primer_task == pick_pcr_primers ||
748 pa->primer_task == pick_pcr_primers_and_hyb_probe)
749 && right - left < pr_min - 1) {
750 sa->pair_expl.product = 1;
751 sa->pair_expl.considered = 1;
752 return 1;
753 } else return 0;
754 }
755
756 /*
757 * Make complete list of acceptable internal oligos in mid. Place the number
758 * of valid elements in mid in *n_m. Return 1 if there are no acceptable
759 * internal oligos; otherwise return 0.
760 */
761 static int
make_internal_oligos_list(pa,sa,n_m,local_args,end_args,local_end_args,ctx)762 make_internal_oligos_list(pa, sa, n_m, local_args, end_args, local_end_args, ctx)
763 const primer_args *pa;
764 seq_args *sa;
765 int *n_m;
766 const dpal_args *local_args, *end_args, *local_end_args;
767 Primer3Context * ctx;
768 {
769 int i, j, n, k;
770
771 char s[MAX_PRIMER_LENGTH+1];
772 primer_rec h;
773
774 if (NULL == (ctx->mid)) {
775 (ctx->mid_len) = INITIAL_LIST_LEN;
776 (ctx->mid) = pr_safe_malloc(sizeof(*(ctx->mid)) * (ctx->mid_len));
777 }
778
779 n = strlen(sa->trimmed_seq);
780 k = 0;
781 for(i = n - 1; i >= pa->io_primer_min_size-1; i--) {
782 s[0] = '\0';
783 for(j = pa->io_primer_min_size; j <=pa->io_primer_max_size; j++) {
784 if(i-j < -1) break;
785 if (k >= (ctx->mid_len)) {
786 (ctx->mid_len) += ((ctx->mid_len) >> 1);
787 (ctx->mid) = pr_safe_realloc((ctx->mid), (ctx->mid_len) * sizeof(*(ctx->mid)));
788 }
789 h.start = i - j +1;
790 h.length = j;
791 h.repeat_sim.score = NULL;
792 _pr_substr(sa->trimmed_seq, h.start, h.length,s);
793
794 if (sa->internal_input && strcmp_nocase(sa->internal_input, s))
795 continue;
796 h.must_use = (sa->internal_input && pa->pick_anyway);
797
798 h.repeat_sim.score = NULL;
799 oligo_param(pa, &h, OT_INTL, local_args, end_args,
800 local_end_args, sa, &sa->intl_expl, ctx);
801 sa->intl_expl.considered++;
802 if (OK_OR_MUST_USE(&h)) {
803 h.quality = p_obj_fn(pa, &h, 2);
804 (ctx->mid)[k] = h;
805 k++;
806 } else if (h.ok==OV_TOO_MANY_NS || h.ok==OV_INTERSECT_TARGET
807 || h.ok==OV_SELF_ANY || h.ok==OV_POLY_X
808 || h.ok==OV_EXCL_REGION || h.ok==OV_GC_CLAMP
809 || h.ok==OV_SEQ_QUALITY || h.ok==OV_LIB_SIM ) {
810 /* Break from the inner for loop, because there is no
811 legal longer oligo with the same 3' sequence. */
812 break;
813 }
814 }
815 }
816 *n_m = k;
817 sa->intl_expl.ok = *n_m;
818 if(*n_m==0)return 1;
819 else return 0;
820 }
821
822 /*
823 * Compute various characteristics of the oligo, and determine
824 * if it is acceptable.
825 */
826 #define OUTSIDE_START_WT 30.0
827 #define INSIDE_START_WT 20.0
828 #define INSIDE_STOP_WT 100.0
829 #define OUTSIDE_STOP_WT 0.5
830 static void
oligo_param(pa,h,l,local_args,end_args,local_end_args,sa,stats,ctx)831 oligo_param(pa, h, l, local_args, end_args, local_end_args, sa, stats, ctx)
832 const primer_args *pa;
833 primer_rec *h;
834 oligo_type l;
835 const dpal_args *local_args, *end_args, *local_end_args;
836 seq_args *sa;
837 oligo_stats *stats;
838 Primer3Context * ctx;
839 {
840 int i,j,k, min_q, min_end_q;
841 int poly_x, max_poly_x;
842 int must_use = h->must_use;
843 const char *seq = sa->trimmed_seq;
844 char s1[MAX_PRIMER_LENGTH+1], s1_rev[MAX_PRIMER_LENGTH+1];
845
846
847
848 h->ok = OV_UNINITIALIZED;
849 h->target = h->gc_content = h->num_ns=h->excl=0;
850
851 h->template_mispriming = h->template_mispriming_r = ALIGN_SCORE_UNDEF;
852
853 PR_ASSERT(OT_LEFT == l || OT_RIGHT == l || OT_INTL == l);
854
855 if (OT_LEFT == l || OT_INTL == l) {j = h->start; k=j+h->length-1;}
856 else {j = h->start-h->length+1; k=h->start;}
857
858 PR_ASSERT(k >= 0);
859 PR_ASSERT(k < TRIMMED_SEQ_LEN(sa));
860
861 /* edited by T. Koressaar for lowercase masking */
862 if(pa->lowercase_masking==1) {
863 if(l==OT_LEFT) {
864 check_if_lowercase_masked(k, sa->trimmed_orig_seq,h);
865 }
866 if(l==OT_RIGHT) {
867 check_if_lowercase_masked(j, sa->trimmed_orig_seq,h);
868 }
869 if(h->ok==OV_GMASKED) {
870 stats->gmasked++;
871 if (!must_use) return;
872 }
873 }
874 /* end T. Koressar's changes */
875
876 gc_and_n_content(j, k-j+1, sa->trimmed_seq, h);
877
878 if (((OT_LEFT == l || OT_RIGHT == l) &&
879 h->num_ns > pa->num_ns_accepted) ||
880 (OT_INTL == l && h->num_ns > pa->io_num_ns_accepted) ) {
881 h->ok = OV_TOO_MANY_NS;
882 stats->ns++;
883 if (!must_use) return;
884 }
885
886 /* Upstream error checking has ensured that we use non-default position
887 penalties only when there is 0 or 1 target. */
888 PR_ASSERT(sa->num_targets <= 1 || _PR_DEFAULT_POSITION_PENALTIES(pa));
889 if (l < 2
890 && _PR_DEFAULT_POSITION_PENALTIES(pa)
891 && oligo_overlaps_interval(j, k-j+1, sa->tar, sa->num_targets)) {
892 h->position_penalty = 0.0;
893 h->position_penalty_infinite = '\1';
894 h->target = 1;
895 } else if (l < 2 && !_PR_DEFAULT_POSITION_PENALTIES(pa)
896 && 1 == sa->num_targets) {
897 compute_position_penalty(pa, sa, h, l);
898 if (h->position_penalty_infinite) h->target = 1;
899 } else {
900 h->position_penalty = 0.0;
901 h->position_penalty_infinite = '\0';
902 }
903
904 if (!PR_START_CODON_POS_IS_NULL(sa)) {
905 if (OT_LEFT == l) {
906 if (sa->start_codon_pos > h->start)
907 h->position_penalty
908 = (sa->start_codon_pos - h->start) * OUTSIDE_START_WT;
909 else
910 h->position_penalty
911 = (h->start - sa->start_codon_pos) * INSIDE_START_WT;
912 }
913 else if (OT_RIGHT == l) {
914 if (-1 == sa->stop_codon_pos) {
915 h->position_penalty = (TRIMMED_SEQ_LEN(sa) - h->start - 1) * INSIDE_STOP_WT;
916 } else if (sa->stop_codon_pos < h->start)
917 h->position_penalty
918 = (h->start - sa->stop_codon_pos) * OUTSIDE_STOP_WT;
919 else
920 h->position_penalty
921 = (sa->stop_codon_pos - h->start) * INSIDE_STOP_WT;
922 }
923 }
924
925 if (l < 2 && oligo_overlaps_interval(j, k-j+1, sa->excl, sa->num_excl))
926 h->excl = 1;
927
928 if (l == 2 && oligo_overlaps_interval(j, k-j+1, sa->excl_internal,
929 sa->num_internal_excl))
930 h->excl = 1;
931
932 if(l < 2 && h->target==1) {
933 h->ok = OV_INTERSECT_TARGET;
934 stats->target++;
935 if (!must_use) return;
936 }
937
938 if(h->excl==1){
939 h->ok = OV_EXCL_REGION;
940 stats->excluded++;
941 if (!must_use) return;
942 }
943
944 if((l<2 && (h->gc_content< pa->min_gc || h->gc_content > pa->max_gc)) ||
945 (l==2 && (h->gc_content< pa->io_min_gc || h->gc_content > pa->io_max_gc))){
946 h->ok = OV_GC_CONTENT;
947 stats->gc++;
948 if (!must_use) return;
949 }
950 if(pa->gc_clamp != 0){
951 if(OT_LEFT == l){
952 for(i=k-pa->gc_clamp+1; i<= k; i++)if(seq[i] !='G'&&seq[i] !='C'){
953 h->ok = OV_GC_CLAMP;
954 stats->gc_clamp++;
955 if (!must_use) return; else break;
956 }
957 }
958 if(OT_RIGHT == l){
959 for(i=j; i<j+pa->gc_clamp; i++)if(seq[i] != 'G' && seq[i] != 'C'){
960 h->ok = OV_GC_CLAMP;
961 stats->gc_clamp++;
962 if (!must_use) return; else break;
963 }
964 }
965 }
966
967 check_sequence_quality(pa, h, l, sa, j, k, &min_q, &min_end_q);
968 if(OT_LEFT == l || OT_RIGHT == l) {
969 if (min_q < pa->min_quality) {
970 h->ok = OV_SEQ_QUALITY;
971 stats->seq_quality++;
972 if (!must_use) return;
973 } else if (min_end_q < pa->min_end_quality) {
974 h->ok = OV_SEQ_QUALITY;
975 stats->seq_quality++;
976 if (!must_use) return;
977 }
978 } else if (OT_INTL == l) {
979 if (min_q < pa->io_min_quality) {
980 h->ok = OV_SEQ_QUALITY;
981 stats->seq_quality++;
982 if (!must_use) return;
983 }
984 } else {
985 PR_ASSERT(0); /* Programming error. */
986 }
987
988 if(OT_LEFT == l || OT_RIGHT == l) max_poly_x = pa->max_poly_x;
989 else max_poly_x = pa->io_max_poly_x;
990 if(max_poly_x > 0) {
991 poly_x = 1;
992 for(i=j+1;i<=k;i++){
993 if(seq[i] == seq[i-1]||seq[i] == 'N'){
994 poly_x++;
995 if(poly_x > max_poly_x){
996 h->ok = OV_POLY_X;
997 stats->poly_x++;
998 if (!must_use) return; else break;
999 }
1000 }
1001 else poly_x = 1;
1002 }
1003 }
1004
1005 _pr_substr(seq,j,k-j+1,s1);
1006
1007 if(OT_LEFT == l || OT_RIGHT == l)
1008 h->temp
1009 = seqtm(s1, pa->dna_conc, pa->salt_conc, pa->divalent_conc, pa->dntp_conc,
1010 MAX_NN_TM_LENGTH,
1011 pa->tm_santalucia,
1012 pa->salt_corrections);
1013 else
1014 h->temp
1015 = seqtm(s1, pa->io_dna_conc, pa->io_salt_conc, pa->io_divalent_conc, pa->io_dntp_conc,
1016 MAX_NN_TM_LENGTH,
1017 pa->tm_santalucia,
1018 pa->salt_corrections);
1019
1020 if (((l == OT_LEFT || l == OT_RIGHT) && h->temp < pa->min_tm)
1021 || (l==OT_INTL && h->temp<pa->io_min_tm)) {
1022 h->ok = OV_TM_LOW;
1023 stats->temp_min++;
1024 if (!must_use) return;
1025 }
1026 if (((OT_LEFT == l || OT_RIGHT == l) && h->temp>pa->max_tm)
1027 || (OT_INTL == l && h->temp>pa->io_max_tm)) {
1028 h->ok = OV_TM_HIGH;
1029 stats->temp_max++;
1030 if (!must_use) return;
1031 }
1032 if (OT_LEFT == l) {
1033 if ((h->end_stability = end_oligodg(s1, 5,
1034 pa->tm_santalucia))
1035 > pa->max_end_stability) {
1036 h->ok = OV_END_STAB;
1037 stats->stability++;
1038 if (!must_use) return;
1039 }
1040 } else if (OT_RIGHT == l) {
1041 _pr_reverse_complement(s1, s1_rev);
1042 if ((h->end_stability = end_oligodg(s1_rev, 5,
1043 pa->tm_santalucia))
1044 > pa->max_end_stability) {
1045 h->ok = OV_END_STAB;
1046 stats->stability++;
1047 if (!must_use) return;
1048 }
1049 }
1050
1051 if (must_use
1052 || pa->file_flag
1053 || (pa->primer_task != pick_pcr_primers &&
1054 pa->primer_task != pick_pcr_primers_and_hyb_probe)
1055 || ((OT_RIGHT == l || OT_LEFT == l) &&
1056 (pa->primer_weights.compl_any || pa->primer_weights.compl_end))
1057 || (OT_INTL == l
1058 && (pa->io_weights.compl_any || pa->io_weights.compl_end))) {
1059
1060 oligo_compl(h, pa, sa, l, local_args, end_args, local_end_args);
1061
1062 if (h->ok != OV_UNINITIALIZED && !must_use) {
1063 PR_ASSERT(h->ok != OV_OK);
1064 return;
1065 }
1066
1067 } else
1068 h->self_any = h->self_end = ALIGN_SCORE_UNDEF;
1069
1070 if (must_use
1071 || pa->file_flag
1072 ||(pa->primer_task != pick_pcr_primers &&
1073 pa->primer_task != pick_pcr_primers_and_hyb_probe)
1074 || ((OT_RIGHT == l || OT_LEFT == l) && pa->primer_weights.repeat_sim)
1075 || ((OT_RIGHT == l || OT_LEFT == l)
1076 && pa->primer_weights.template_mispriming)
1077 || (OT_INTL == l && pa->io_weights.repeat_sim)) {
1078
1079 oligo_mispriming(h, pa, sa, l, local_end_args, ctx);
1080
1081 }
1082
1083 if (OV_UNINITIALIZED == h->ok) h->ok = OV_OK;
1084 }
1085 #undef OUTSIDE_START_WT
1086 #undef INSIDE_START_WT
1087 #undef INSIDE_STOP_WT
1088 #undef OUTSIDE_STOP_WT
1089
1090 static void
check_sequence_quality(pa,h,l,sa,j,k,r_min_q,r_min_q_end)1091 check_sequence_quality(pa, h, l, sa, j, k, r_min_q, r_min_q_end)
1092 const primer_args *pa;
1093 primer_rec *h;
1094 oligo_type l;
1095 const seq_args *sa;
1096 int j, k;
1097 int *r_min_q, *r_min_q_end;
1098 {
1099 int i,min_q, min_q_end, m, q;
1100
1101 q = min_q = min_q_end = pa->quality_range_max;
1102
1103 if (NULL != sa->quality) {
1104
1105 if(OT_LEFT == l || OT_RIGHT == l){
1106 min_q = pa->min_quality;
1107 min_q_end = pa->min_end_quality;
1108 }
1109 else {min_q = min_q_end = pa->io_min_quality;}
1110
1111 if (OT_LEFT == l || OT_INTL == l) {
1112
1113 for(i = k-4; i <= k; i++) {
1114 if(i < j) continue;
1115 m = sa->quality[i + sa->incl_s];
1116 if (m < q) q = m;
1117 }
1118 min_q_end = q;
1119
1120 for(i = j; i<=k-5; i++) {
1121 m = sa->quality[i + sa->incl_s];
1122 if (m < q) q = m;
1123 }
1124 min_q = q;
1125
1126 } else if (OT_RIGHT == l) {
1127 for(i = j; i < j+5; i++) {
1128 if(i > k) break;
1129 m = sa->quality[i + sa->incl_s];
1130 if (m < q) q = m;
1131 }
1132 min_q_end = q;
1133
1134 for(i = j+5; i <= k; i++) {
1135 m = sa->quality[i + sa->incl_s];
1136 if (m < q) q = m;
1137 }
1138 min_q = q;
1139 } else {
1140 PR_ASSERT(0); /* Programming error. */
1141 }
1142 }
1143 h->seq_quality = *r_min_q = min_q;
1144 *r_min_q_end = min_q_end;
1145 }
1146
1147 /*
1148 * Set h->gc_content to the GC % of the 'len' bases starting at 'start' in
1149 * 'sequence' (excluding 'N's). Set h->num_ns to the number of 'N's in that
1150 * subsequence.
1151 */
1152 static void
gc_and_n_content(start,len,sequence,h)1153 gc_and_n_content(start, len, sequence, h)
1154 const int start, len;
1155 const char *sequence;
1156 primer_rec *h;
1157 {
1158 const char* p = &sequence[start];
1159 const char* stop = p + len;
1160 int num_gc = 0, num_gcat = 0, num_n = 0;
1161 while (p < stop) {
1162 if ('N' == *p)
1163 num_n++;
1164 else {
1165 num_gcat++;
1166 if ('C' == *p || 'G' == *p) num_gc++;
1167 }
1168 p++;
1169 }
1170 h->num_ns = num_n;
1171 if (0 == num_gcat) h->gc_content= 0.0;
1172 else h->gc_content = 100.0 * ((double)num_gc)/num_gcat;
1173 }
1174
1175 static int
oligo_overlaps_interval(start,len,intervals,num_intervals)1176 oligo_overlaps_interval(start, len, intervals, num_intervals)
1177 const int start, len;
1178 interval_array_t intervals;
1179 const int num_intervals;
1180 {
1181 int i;
1182 int last = start + len - 1;
1183 for (i = 0; i < num_intervals; i++)
1184 if (!(last < intervals[i][0]
1185 || start > (intervals[i][0] + intervals[i][1] - 1)))
1186 return 1;
1187 return 0;
1188 }
1189
1190 /* Calculate the part of the objective function due to one primer. */
1191 static double
p_obj_fn(pa,h,j)1192 p_obj_fn(pa, h, j)
1193 const primer_args *pa;
1194 primer_rec *h;
1195 int j;
1196 {
1197 double sum;
1198
1199 sum = 0;
1200 if (j == OT_LEFT || j == OT_RIGHT) {
1201 if(pa->primer_weights.temp_gt && h->temp > pa->opt_tm)
1202 sum += pa->primer_weights.temp_gt * (h->temp - pa->opt_tm);
1203 if(pa->primer_weights.temp_lt && h->temp < pa->opt_tm)
1204 sum += pa->primer_weights.temp_lt * (pa->opt_tm - h->temp);
1205
1206 if(pa->primer_weights.gc_content_gt && h->gc_content > pa->opt_gc_content)
1207 sum += pa->primer_weights.gc_content_gt
1208 * (h->gc_content - pa->opt_gc_content);
1209 if(pa->primer_weights.gc_content_lt && h->gc_content < pa->opt_gc_content)
1210 sum += pa->primer_weights.gc_content_lt
1211 * (pa->opt_gc_content - h->gc_content);
1212
1213 if(pa->primer_weights.length_lt && h->length < pa->primer_opt_size)
1214 sum += pa->primer_weights.length_lt * (pa->primer_opt_size - h->length);
1215 if(pa->primer_weights.length_gt && h->length > pa->primer_opt_size)
1216 sum += pa->primer_weights.length_gt * (h->length - pa->primer_opt_size);
1217 if(pa->primer_weights.compl_any)
1218 sum += pa->primer_weights.compl_any * h->self_any
1219 / PR_ALIGN_SCORE_PRECISION;
1220 if(pa->primer_weights.compl_end)
1221 sum += pa->primer_weights.compl_end * h->self_end
1222 / PR_ALIGN_SCORE_PRECISION;
1223 if(pa->primer_weights.num_ns)
1224 sum += pa->primer_weights.num_ns * h->num_ns;
1225 if(pa->primer_weights.repeat_sim)
1226 sum += pa->primer_weights.repeat_sim
1227 * h->repeat_sim.score[h->repeat_sim.max]
1228 / PR_ALIGN_SCORE_PRECISION;
1229 if (!h->target) {
1230 /* We might be evaluating p_obj_fn with h->target if
1231 the client supplied 'pick_anyway' and specified
1232 a primer or oligo. */
1233 PR_ASSERT(!h->position_penalty_infinite);
1234 if(pa->primer_weights.pos_penalty)
1235 sum += pa->primer_weights.pos_penalty * h->position_penalty;
1236 }
1237 if(pa->primer_weights.end_stability)
1238 sum += pa->primer_weights.end_stability * h->end_stability;
1239 if(pa->primer_weights.seq_quality)
1240 sum += pa->primer_weights.seq_quality *
1241 (pa->quality_range_max - h->seq_quality);
1242
1243 if (pa->primer_weights.template_mispriming) {
1244 PR_ASSERT(oligo_max_template_mispriming(h) != ALIGN_SCORE_UNDEF);
1245 sum += pa->primer_weights.template_mispriming *
1246 oligo_max_template_mispriming(h);
1247 }
1248
1249 return sum;
1250 } else if (j == OT_INTL) {
1251 if(pa->io_weights.temp_gt && h->temp > pa->io_opt_tm)
1252 sum += pa->io_weights.temp_gt * (h->temp - pa->io_opt_tm);
1253 if(pa->io_weights.temp_lt && h->temp < pa->io_opt_tm)
1254 sum += pa->io_weights.temp_lt * (pa->io_opt_tm - h->temp);
1255
1256 if(pa->io_weights.gc_content_gt && h->gc_content > pa->io_opt_gc_content)
1257 sum += pa->io_weights.gc_content_gt
1258 * (h->gc_content - pa->io_opt_gc_content);
1259 if(pa->io_weights.gc_content_lt && h->gc_content < pa->io_opt_gc_content)
1260 sum += pa->io_weights.gc_content_lt
1261 * (pa->io_opt_gc_content - h->gc_content);
1262
1263 if(pa->io_weights.length_lt && h->length < pa->io_primer_opt_size)
1264 sum += pa->io_weights.length_lt * (pa->io_primer_opt_size - h->length);
1265 if(pa->io_weights.length_gt && h->length > pa->io_primer_opt_size)
1266 sum += pa->io_weights.length_gt * (h->length - pa->io_primer_opt_size);
1267 if(pa->io_weights.compl_any)
1268 sum += pa->io_weights.compl_any * h->self_any / PR_ALIGN_SCORE_PRECISION;
1269 if(pa->io_weights.compl_end)
1270 sum += pa->io_weights.compl_end * h->self_end / PR_ALIGN_SCORE_PRECISION;
1271 if(pa->io_weights.num_ns)
1272 sum += pa->io_weights.num_ns * h->num_ns;
1273 if(pa->io_weights.repeat_sim)
1274 sum += pa->io_weights.repeat_sim
1275 * h->repeat_sim.score[h->repeat_sim.max]
1276 / PR_ALIGN_SCORE_PRECISION;
1277 if(pa->io_weights.seq_quality)
1278 sum += pa->io_weights.seq_quality *
1279 (pa->quality_range_max - h->seq_quality);
1280 return sum;
1281 } else {
1282 PR_ASSERT(0); /* Programmig error. */
1283 }
1284 }
1285
1286 /* Return max of h->template_mispriming and h->template_mispriming_r (max
1287 template mispriming on either strand). */
1288 static short
oligo_max_template_mispriming(h)1289 oligo_max_template_mispriming(h)
1290 const primer_rec *h;
1291 {
1292 return h->template_mispriming > h->template_mispriming_r ?
1293 h->template_mispriming : h->template_mispriming_r;
1294 }
1295
1296
1297 static int
choose_pair(pa,sa,n_f,n_r,n_m,int_num,local_args,end_args,local_end_args,p,cancelFlag,progress,ctx)1298 choose_pair(pa, sa, n_f, n_r, n_m, int_num, local_args, end_args,
1299 local_end_args, p, cancelFlag, progress, ctx)
1300 const primer_args *pa;
1301 seq_args *sa;
1302 int n_f;
1303 int n_r;
1304 int n_m;
1305 int int_num;
1306 const dpal_args *local_args, *end_args, *local_end_args;
1307 pair_array_t *p;
1308 int * cancelFlag;
1309 int * progress;
1310 Primer3Context * ctx;
1311 {
1312 int i,j;
1313 int k; /*
1314 * The number of acceptable primer pairs saved in p.
1315 * (k <= pa->num_return.)
1316 */
1317 int i0, n_last, n_int;
1318 primer_pair worst_pair; /* The worst pair among those being "remembered". */
1319 int i_worst; /* The index within p of worst_pair. */
1320
1321 primer_pair h;
1322
1323 /* Mou must always initialize your variables! */
1324 h.left = NULL;
1325 h.right = NULL;
1326 h.intl = NULL;
1327
1328 k=0;
1329
1330 i_worst = 0;
1331 n_last = n_f;
1332 for(i=0;i<n_r;i++) {
1333 {
1334 const double basePow = 1.0/2;
1335 const double subPow = 1.0/1.01;
1336 *progress = 50 + (int)(50*(1 - pow(basePow, int_num)*(basePow + (1 - basePow)*pow(subPow, i))));
1337 }
1338
1339 /*
1340 * Make a quick cut based on the the quality of the best left
1341 * primer.
1342 *
1343 * worst_pair must be defined if k >= pa->num_return.
1344 */
1345 if (!OK_OR_MUST_USE(&(ctx->r)[i])) continue;
1346 if (k >= pa->num_return
1347 && ((ctx->r)[i].quality+(ctx->f)[0].quality > worst_pair.pair_quality
1348 || worst_pair.pair_quality == 0))
1349 break;
1350
1351 for(j=0;j<n_last;j++) {
1352 if( *cancelFlag ) {
1353 break;
1354 }
1355 /*
1356 * Invariant: if 2 pairs in p have the same pair_quality, then the
1357 * pair discovered second will have a higher index within p.
1358 *
1359 * This invariant is needed to produce consistent results when 2
1360 * pairs have the same pair_quality.
1361 */
1362
1363 if (!OK_OR_MUST_USE(&(ctx->r)[i])) break;
1364 if (!OK_OR_MUST_USE(&(ctx->f)[j])) continue;
1365 if(k>= pa->num_return
1366 && ((ctx->f)[j].quality+(ctx->r)[i].quality > worst_pair.pair_quality
1367 || worst_pair.pair_quality == 0)) {
1368 /* worst_pair must be defined if k >= pa->num_return. */
1369 n_last=j;
1370 break;
1371 }
1372
1373 if (PAIR_OK ==
1374 pair_param(pa, sa, j, i, int_num, &h, local_args, end_args,
1375 local_end_args, ctx)) {
1376
1377 if (!pa->pr_pair_weights.io_quality) {
1378 h.pair_quality = obj_fn(pa, &h);
1379 PR_ASSERT(h.pair_quality >= 0.0);
1380 }
1381
1382
1383 if ( pa->primer_task == pick_pcr_primers_and_hyb_probe
1384 && choose_internal_oligo(h.left, h.right,
1385 n_m, &n_int, sa, pa, local_args,
1386 end_args, local_end_args, ctx)!=0) {
1387 sa->pair_expl.internal++;
1388 continue;
1389 }
1390 sa->pair_expl.ok++;
1391 if (k < pa->num_return) {
1392 if ( pa->primer_task == pick_pcr_primers_and_hyb_probe)
1393 h.intl = &(ctx->mid)[n_int];
1394 if(pa->pr_pair_weights.io_quality) {
1395 h.pair_quality = obj_fn(pa, &h);
1396 PR_ASSERT(h.pair_quality >= 0.0);
1397 }
1398
1399 add_pair(&h, p);
1400 /* p->pairs[k] = h; */
1401 if (k == 0 || primer_pair_comp(&h, &worst_pair) > 0){
1402 worst_pair = h;
1403 i_worst = k;
1404 }
1405 k++;
1406 } else {
1407 if ( pa->primer_task == pick_pcr_primers_and_hyb_probe)
1408 h.intl = &(ctx->mid)[n_int];
1409 if(pa->pr_pair_weights.io_quality) {
1410 h.pair_quality = obj_fn(pa, &h);
1411 PR_ASSERT(h.pair_quality >= 0.0);
1412 }
1413
1414 if (primer_pair_comp(&h, &worst_pair) < 0) {
1415 /*
1416 * There are already pa->num_return results, and vl is better than
1417 * the pa->num_return the quality found so far.
1418 */
1419 p->pairs[i_worst] = h;
1420 worst_pair = h; /* h is a lower bound on the worst pair. */
1421 for (i0 = 0; i0<pa->num_return; i0++)
1422 if(primer_pair_comp(&p->pairs[i0], &worst_pair) > 0) {
1423 i_worst = i0;
1424 worst_pair = p->pairs[i0];
1425
1426 }
1427 }
1428 }
1429 }
1430 }
1431 }
1432 if(k!=0) qsort(p->pairs, k, sizeof(primer_pair), primer_pair_comp);
1433 p->num_pairs = k;
1434 if (k==0) return 1;
1435 else return 0;
1436 }
1437
1438 /* Choose best internal oligo for given pair of left and right primers. */
1439 static int
choose_internal_oligo(left,right,n_m,nm,sa,pa,local_args,end_args,local_end_args,ctx)1440 choose_internal_oligo(left, right, n_m, nm, sa, pa, local_args,
1441 end_args, local_end_args, ctx)
1442 const primer_rec *left, *right;
1443 int n_m;
1444 int *nm;
1445 seq_args *sa;
1446 const primer_args *pa;
1447 const dpal_args *local_args, *end_args, *local_end_args;
1448 Primer3Context * ctx;
1449 {
1450 int i,k;
1451 double min;
1452
1453 min = 1000000.;
1454 i = -1;
1455 for(k=0; k < n_m; k++){
1456 if (((ctx->mid)[k].start > (left->start + (left->length-1)))
1457 && (((ctx->mid)[k].start+((ctx->mid)[k].length-1)) < (right->start-right->length+1))
1458 && ((ctx->mid)[k].quality < min)
1459 && (OK_OR_MUST_USE(&(ctx->mid)[k]))) {
1460 if((ctx->mid)[k].self_any == ALIGN_SCORE_UNDEF){
1461 oligo_compl(&(ctx->mid)[k], pa, sa, OT_INTL, local_args,
1462 end_args, local_end_args);
1463 if (!OK_OR_MUST_USE(&(ctx->mid)[k])) continue;
1464 }
1465 if((ctx->mid)[k].repeat_sim.score == NULL) {
1466 oligo_mispriming(&(ctx->mid)[k], pa, sa, OT_INTL, local_args, ctx);
1467 if (!OK_OR_MUST_USE(&(ctx->mid)[k])) continue;
1468 }
1469 min = (ctx->mid)[k].quality;
1470 i=k;
1471 }
1472 }
1473 *nm = i;
1474 if(*nm < 0) return 1;
1475 return 0;
1476 }
1477
1478 /* Compare function for sorting primer records. */
1479 static int
primer_rec_comp(x1,x2)1480 primer_rec_comp(x1, x2)
1481 const void *x1, *x2;
1482 {
1483 const primer_rec *a1 = x1, *a2 = x2;
1484
1485 if(a1->quality < a2->quality) return -1;
1486 if (a1->quality > a2->quality) return 1;
1487
1488 /*
1489 * We want primer_rec_comp to always return a non-0 result, because
1490 * different implementations of qsort differ in how they treat "equal"
1491 * elements, making it difficult to compare test output on different
1492 * systems.
1493 */
1494 if(a1->start > a2->start) return -1;
1495 if(a1->start < a2->start) return 1;
1496 if(a1->length < a2->length) return -1;
1497 return 1;
1498 }
1499
1500 /* Compare function for sorting primer records. */
1501 static int
primer_pair_comp(x1,x2)1502 primer_pair_comp(x1, x2)
1503 const void *x1, *x2;
1504 {
1505 const primer_pair *a1 = x1, *a2 = x2;
1506 int y1, y2;
1507
1508 if(a1->pair_quality < a2->pair_quality) return -1;
1509 if (a1->pair_quality > a2->pair_quality) return 1;
1510
1511 if (a1->compl_measure < a2->compl_measure) return -1;
1512 if (a1->compl_measure > a2->compl_measure) return 1;
1513
1514 /*
1515 * The following statements ensure that sorting
1516 * produces the same order on all systems regardless
1517 * of whether the sorting function is stable.
1518 */
1519
1520 y1 = a1->left->start;
1521 y2 = a2->left->start;
1522 if (y1 > y2) return -1; /* prefer left primers to the right. */
1523 if (y1 < y2) return 1;
1524
1525 y1 = a1->right->start;
1526 y2 = a2->right->start;
1527 if (y1 < y2) return -1; /* prefer right primers to the left. */
1528 if (y1 > y2) return 1;
1529
1530 y1 = a1->left->length;
1531 y2 = a2->left->length;
1532 if (y1 < y2) return -1; /* prefer shorter primers. */
1533 if (y1 > y2) return 1;
1534
1535 y1 = a1->right->length;
1536 y2 = a2->right->length;
1537 if (y1 < y2) return -1; /* prefer shorter primers. */
1538 if (y1 > y2) return 1;
1539
1540 return 0;
1541 }
1542
1543 /*
1544 * Defines parameter values for given primer pair. Returns PAIR_OK if the pair is
1545 * acceptable; PAIR_FAILED otherwise.
1546 */
1547 static int
pair_param(pa,sa,m,n,int_num,h,local_args,end_args,local_end_args,ctx)1548 pair_param(pa, sa, m, n, int_num, h, local_args, end_args, local_end_args, ctx)
1549 const primer_args *pa;
1550 seq_args *sa;
1551 int m, n, int_num;
1552 primer_pair *h;
1553 const dpal_args *local_args, *end_args, *local_end_args;
1554 Primer3Context * ctx;
1555 {
1556 char s1[MAX_PRIMER_LENGTH+1], s2[MAX_PRIMER_LENGTH+1],
1557 s1_rev[MAX_PRIMER_LENGTH+1], s2_rev[MAX_PRIMER_LENGTH+1];
1558 short compl_end;
1559
1560 /* FUTURE CODE: we must use the pair if the caller specifed
1561 both the left and the right primer. */
1562 int must_use = 0;
1563
1564 int pair_failed_flag = 0;
1565 double min_oligo_tm;
1566
1567 h->left = &(ctx->f)[m];
1568 h->right = &(ctx->r)[n];
1569 h->product_size = (ctx->r)[n].start - (ctx->f)[m].start+1;
1570 h->target = 0;
1571 h->compl_any = h->compl_end = 0;
1572
1573 sa->pair_expl.considered++;
1574
1575 if(h->product_size < pa->pr_min[int_num] ||
1576 h->product_size > pa->pr_max[int_num]) {
1577 sa->pair_expl.product++;
1578 h->product_size = -1;
1579 if (!must_use) return PAIR_FAILED;
1580 else pair_failed_flag = 1;
1581 }
1582
1583 if (sa->num_targets > 0) {
1584 if (pair_spans_target(h, sa))
1585 h->target = 1;
1586 else {
1587 h->target = -1;
1588 sa->pair_expl.target++;
1589 if (!must_use) return PAIR_FAILED;
1590 else pair_failed_flag = 1;
1591 }
1592 }
1593
1594 /* Compute product Tm and related parameters; check constraints. */
1595 h->product_tm
1596 = long_seq_tm(sa->trimmed_seq, h->left->start,
1597 h->right->start - h->left->start + 1, pa->salt_conc, pa->divalent_conc, pa->dntp_conc);
1598
1599 PR_ASSERT(h->product_tm != OLIGOTM_ERROR);
1600
1601 min_oligo_tm
1602 = h->left->temp > h->right->temp ? h->right->temp : h->left->temp;
1603 h->product_tm_oligo_tm_diff = h->product_tm - min_oligo_tm;
1604 h->t_opt_a = 0.3 * min_oligo_tm + 0.7 * h->product_tm - 14.9;
1605
1606 if (pa->product_min_tm != PR_DEFAULT_PRODUCT_MIN_TM
1607 && h->product_tm < pa->product_min_tm) {
1608 sa->pair_expl.low_tm++;
1609 if (!must_use) return PAIR_FAILED;
1610 else pair_failed_flag = 1;
1611 }
1612
1613 if (pa->product_max_tm != PR_DEFAULT_PRODUCT_MAX_TM
1614 && h->product_tm > pa->product_max_tm) {
1615 sa->pair_expl.high_tm++;
1616 if (!must_use) return PAIR_FAILED;
1617 else pair_failed_flag = 1;
1618 }
1619
1620 h->diff_tm = fabs((ctx->f)[m].temp - (ctx->r)[n].temp);
1621 if (h->diff_tm > pa->max_diff_tm) {
1622 sa->pair_expl.temp_diff++;
1623 if (!must_use) return PAIR_FAILED;
1624 else pair_failed_flag = 1;
1625 }
1626
1627 _pr_substr(sa->trimmed_seq,(ctx->f)[m].start,(ctx->f)[m].length,s1);
1628 _pr_substr(sa->trimmed_seq,(ctx->r)[n].start-(ctx->r)[n].length+1,(ctx->r)[n].length,s2);
1629 if((ctx->f)[m].self_any == ALIGN_SCORE_UNDEF){
1630 oligo_compl(&(ctx->f)[m], pa, sa, OT_LEFT, local_args, end_args, local_end_args);
1631 if (!OK_OR_MUST_USE(&(ctx->f)[m])) {
1632 sa->pair_expl.considered--;
1633 return PAIR_FAILED;
1634 }
1635 }
1636 if((ctx->f)[m].repeat_sim.score == NULL){
1637 oligo_mispriming(&(ctx->f)[m], pa, sa, OT_LEFT, local_end_args, ctx);
1638 if (!OK_OR_MUST_USE(&(ctx->f)[m])) {
1639 sa->pair_expl.considered--;
1640 return PAIR_FAILED;
1641 }
1642 }
1643 if((ctx->r)[n].self_any == ALIGN_SCORE_UNDEF){
1644 oligo_compl(&(ctx->r)[n], pa, sa, OT_RIGHT, local_args, end_args, local_end_args);
1645 if (!OK_OR_MUST_USE(&(ctx->r)[n])) {
1646 sa->pair_expl.considered--;
1647 return PAIR_FAILED;
1648 }
1649 }
1650 if((ctx->r)[n].repeat_sim.score == NULL){
1651 oligo_mispriming(&(ctx->r)[n], pa, sa, OT_RIGHT, local_end_args, ctx);
1652 if (!OK_OR_MUST_USE(&(ctx->r)[n])) {
1653 sa->pair_expl.considered--;
1654 return PAIR_FAILED;
1655 }
1656 }
1657
1658 /*
1659 * Similarity between s1 and s2 is equivalent to complementarity between
1660 * s2's complement and s1. (Both s1 and s2 are taken from the same strand.)
1661 */
1662 h->compl_any = align(s1,s2,local_args);
1663 if (h->compl_any > pa->self_any) {
1664 sa->pair_expl.compl_any++;
1665 return PAIR_FAILED;
1666 }
1667
1668 if ((h->compl_end = align(s1,s2,end_args)) > pa->self_end) {
1669 sa->pair_expl.compl_end++;
1670 return PAIR_FAILED;
1671 }
1672 /*
1673 * It is conceivable (though very unlikely in practice) that
1674 * align(s2_rev, s1_rev, end_args) > align(s1,s2,end_args).
1675 */
1676 _pr_reverse_complement(s1, s1_rev);
1677 _pr_reverse_complement(s2, s2_rev);
1678 if((compl_end = align(s2_rev, s1_rev, end_args)) > h->compl_end) {
1679 if (compl_end > pa->self_end) {
1680 sa->pair_expl.compl_end++;
1681 return PAIR_FAILED;
1682 }
1683 h->compl_end = compl_end;
1684 }
1685 h->compl_measure =
1686 (h->right->self_end + h->left->self_end + h->compl_end) * 1.1
1687 + h->right->self_any + h->left->self_any + h->compl_any;
1688 if((h->repeat_sim = pair_repeat_sim(h, pa)) > pa->pair_repeat_compl){
1689 sa->pair_expl.repeat_sim++;
1690 return PAIR_FAILED;
1691 }
1692
1693 if (!_pr_need_pair_template_mispriming(pa))
1694 h->template_mispriming = ALIGN_SCORE_UNDEF;
1695 else {
1696 PR_ASSERT(h->left->template_mispriming != ALIGN_SCORE_UNDEF);
1697 PR_ASSERT(h->left->template_mispriming_r != ALIGN_SCORE_UNDEF);
1698 PR_ASSERT(h->right->template_mispriming != ALIGN_SCORE_UNDEF);
1699 PR_ASSERT(h->right->template_mispriming_r != ALIGN_SCORE_UNDEF);
1700 h->template_mispriming =
1701 h->left->template_mispriming + h->right->template_mispriming_r;
1702 if ((h->left->template_mispriming_r + h->right->template_mispriming)
1703 > h->template_mispriming)
1704 h->template_mispriming
1705 = h->left->template_mispriming_r + h->right->template_mispriming;
1706
1707 if (pa->pair_max_template_mispriming >= 0.0
1708 && h->template_mispriming > pa->pair_max_template_mispriming) {
1709 sa->pair_expl.template_mispriming++;
1710 return PAIR_FAILED;
1711 }
1712
1713 }
1714 return PAIR_OK;
1715 }
1716
1717 /* compute_position_penalty is experimental code. */
1718 void
compute_position_penalty(pa,sa,h,o_type)1719 compute_position_penalty(pa, sa, h, o_type)
1720 const primer_args *pa;
1721 const seq_args *sa;
1722 primer_rec *h;
1723 oligo_type o_type;
1724 {
1725 int three_prime_base;
1726 int inside_flag = 0;
1727 int target_begin, target_end;
1728
1729 PR_ASSERT(OT_LEFT == o_type || OT_RIGHT == o_type);
1730 PR_ASSERT(1 == sa->num_targets);
1731 target_begin = sa->tar[0][0];
1732 target_end = target_begin + sa->tar[0][1] - 1;
1733
1734 three_prime_base = OT_LEFT == o_type
1735 ? h->start + h->length - 1 : h->start - h->length + 1;
1736 h->position_penalty_infinite = '\1';
1737 h->position_penalty = 0.0;
1738
1739 if (OT_LEFT == o_type) {
1740 if (three_prime_base <= target_end) {
1741 h->position_penalty_infinite = '\0';
1742 if (three_prime_base < target_begin)
1743 h->position_penalty = target_begin - three_prime_base - 1;
1744 else {
1745 h->position_penalty = three_prime_base - target_begin + 1;
1746 inside_flag = 1;
1747 }
1748 }
1749 } else { /* OT_RIGHT == o_type */
1750 if (three_prime_base >= target_begin) {
1751 h->position_penalty_infinite = '\0';
1752 if (three_prime_base > target_end) {
1753 h->position_penalty = three_prime_base - target_end - 1;
1754 } else {
1755 h->position_penalty = target_end - three_prime_base + 1;
1756 inside_flag = 1;
1757 }
1758 }
1759 }
1760 if (!inside_flag)
1761 h->position_penalty *= pa->outside_penalty;
1762 else {
1763 h->position_penalty *= pa->inside_penalty;
1764 }
1765 }
1766
1767 /*
1768 * Return 1 if 'pair' spans any target (in sa->tar); otherwise return 0; An
1769 * oligo pair spans a target, t, if the last base of the left primer is to
1770 * left of the last base of t and the first base of the right primer is to
1771 * the right of the first base of t. Of course the primers must
1772 * still be in a legal position with respect to each other.
1773 */
1774 static int
pair_spans_target(pair,sa)1775 pair_spans_target(pair, sa)
1776 const primer_pair *pair;
1777 const seq_args *sa;
1778 {
1779 int i;
1780 int last_of_left = pair->left->start + pair->left->length - 1;
1781 int first_of_right = pair->right->start - pair->right->length + 1;
1782 int target_first, target_last;
1783 for (i = 0; i < sa->num_targets; i++) {
1784 target_first = sa->tar[i][0];
1785 target_last = target_first + sa->tar[i][1] - 1;
1786 if (last_of_left <= target_last
1787 && first_of_right >= target_first
1788 && last_of_left < first_of_right)
1789 return 1;
1790 }
1791 return 0;
1792 }
1793
1794 /*
1795 * Return the value of the objective function value for given primer pair.
1796 * We must know that the pair is acceptable before calling obj_fn.
1797 */
1798 static double
obj_fn(pa,h)1799 obj_fn(pa, h)
1800 const primer_args *pa;
1801 primer_pair *h;
1802 {
1803 double sum;
1804
1805 sum = 0.0;
1806
1807 if(pa->pr_pair_weights.primer_quality)
1808 sum += pa->pr_pair_weights.primer_quality * (h->left->quality + h->right->quality);
1809
1810 if(pa->pr_pair_weights.io_quality &&
1811 pa->primer_task == pick_pcr_primers_and_hyb_probe)
1812 sum += pa->pr_pair_weights.io_quality * h->intl->quality;
1813
1814 if(pa->pr_pair_weights.diff_tm)
1815 sum += pa->pr_pair_weights.diff_tm * h->diff_tm;
1816
1817 if(pa->pr_pair_weights.compl_any)
1818 sum += pa->pr_pair_weights.compl_any * h->compl_any / PR_ALIGN_SCORE_PRECISION;
1819
1820 if(pa->pr_pair_weights.compl_end)
1821 sum += pa->pr_pair_weights.compl_end * h->compl_end / PR_ALIGN_SCORE_PRECISION;
1822
1823 if(pa->pr_pair_weights.product_tm_lt && h->product_tm < pa->product_opt_tm)
1824 sum += pa->pr_pair_weights.product_tm_lt *
1825 (pa->product_opt_tm - h->product_tm);
1826
1827 if(pa->pr_pair_weights.product_tm_gt && h->product_tm > pa->product_opt_tm)
1828 sum += pa->pr_pair_weights.product_tm_gt *
1829 (h->product_tm - pa->product_opt_tm);
1830
1831 if(pa->pr_pair_weights.product_size_lt &&
1832 h->product_size < pa->product_opt_size)
1833 sum += pa->pr_pair_weights.product_size_lt *
1834 (pa->product_opt_size - h->product_size);
1835
1836 if(pa->pr_pair_weights.product_size_gt &&
1837 h->product_size > pa->product_opt_size)
1838 sum += pa->pr_pair_weights.product_size_gt *
1839 (h->product_size - pa->product_opt_size);
1840
1841 if(pa->pr_pair_weights.repeat_sim)
1842 sum += pa->pr_pair_weights.repeat_sim * h->repeat_sim;
1843
1844 if (pa->pr_pair_weights.template_mispriming) {
1845 PR_ASSERT(pa->pr_pair_weights.template_mispriming >= 0.0);
1846 PR_ASSERT(h->template_mispriming >= 0);
1847 sum += pa->pr_pair_weights.template_mispriming * h->template_mispriming;
1848 }
1849
1850 PR_ASSERT(sum >= 0.0);
1851
1852 return sum;
1853 }
1854
1855 char *
pr_gather_warnings(sa,pa)1856 pr_gather_warnings(sa, pa)
1857 const seq_args *sa;
1858 const primer_args *pa;
1859 {
1860 pr_append_str warning;
1861
1862 PR_ASSERT(NULL != sa);
1863 PR_ASSERT(NULL != pa);
1864
1865 warning.data = NULL;
1866 warning.storage_size = 0;
1867
1868 if (pa->repeat_lib.warning.data)
1869 pr_append_new_chunk(&warning, pa->repeat_lib.warning.data);
1870
1871 if(pa->io_mishyb_library.warning.data != NULL) {
1872 pr_append_new_chunk(&warning, pa->io_mishyb_library.warning.data);
1873 pr_append(&warning, " (for internal oligo)");
1874 }
1875
1876 if (sa->warning.data) pr_append_new_chunk(&warning, sa->warning.data);
1877 return pr_is_empty(&warning) ? NULL : warning.data;
1878 }
1879
1880 void
pr_append(x,s)1881 pr_append(x, s)
1882 pr_append_str *x;
1883 const char *s;
1884 {
1885 int xlen, slen;
1886 if (NULL == x->data) {
1887 x->storage_size = 24;
1888 x->data = pr_safe_malloc(x->storage_size);
1889 *x->data = '\0';
1890 }
1891 xlen = strlen(x->data);
1892 slen = strlen(s);
1893 if (xlen + slen + 1 > x->storage_size) {
1894 x->storage_size += 2 * (slen + 1);
1895 x->data = pr_safe_realloc(x->data, x->storage_size);
1896 }
1897 strcpy(x->data + xlen, s);
1898 }
1899
1900 void
pr_append_new_chunk(x,s)1901 pr_append_new_chunk(x, s)
1902 pr_append_str *x;
1903 const char *s;
1904 {
1905 pr_append_w_sep(x, "; ", s);
1906 }
1907
1908 void
pr_append_w_sep(x,sep,s)1909 pr_append_w_sep(x, sep, s)
1910 pr_append_str *x;
1911 const char *sep;
1912 const char *s;
1913 {
1914 if (pr_is_empty(x))
1915 pr_append(x, s);
1916 else {
1917 pr_append(x, sep);
1918 pr_append(x, s);
1919 }
1920 }
1921
1922 void
pr_set_empty(x)1923 pr_set_empty(x)
1924 pr_append_str *x;
1925 {
1926 PR_ASSERT(NULL != x);
1927 if (NULL != x->data) *x->data = '\0';
1928 }
1929
1930 int
pr_is_empty(x)1931 pr_is_empty(x)
1932 const pr_append_str *x;
1933 {
1934 PR_ASSERT(NULL != x);
1935 return NULL == x->data || '\0' == *x->data;
1936 }
1937
1938 static short
align(s1,s2,a)1939 align(s1, s2, a)
1940 const char *s1, *s2;
1941 const dpal_args *a;
1942 {
1943 dpal_results r;
1944
1945 if(a->flag == DPAL_LOCAL || a->flag == DPAL_LOCAL_END) {
1946 if (strlen(s2) < 3) {
1947 /* For extremely short alignments we simply
1948 max out the score, because the dpal subroutines
1949 for these cannot handle this case.
1950 FIX: this can probably be corrected in dpal. */
1951 return (short) (100 * strlen(s2));
1952 }
1953 }
1954 dpal(s1, s2, a, &r);
1955 PR_ASSERT(r.score <= SHRT_MAX);
1956 if (r.score == DPAL_ERROR_SCORE) {
1957 /* There was an error. */
1958 if (errno == ENOMEM) {
1959 OOM_ERROR;
1960 } else {
1961 fprintf(stderr, "%s", r.msg);
1962 /* Fix this later, when error handling
1963 in "primer_choice" is updated to
1964 no longer exit. */
1965 PR_ASSERT(r.score != DPAL_ERROR_SCORE);
1966 }
1967 }
1968 return ((r.score<0) ? 0 : (short)r.score);
1969 }
1970
1971 /* Set dpal args to appropriate values for primer picking. */
1972 /* IMPORTANT -- if the scoring system changes, then
1973 the value returned by align() for short target sequences
1974 must be adjusted. */
1975 static void
set_dpal_args(a)1976 set_dpal_args(a)
1977 dpal_args *a;
1978 {
1979 unsigned int i, j;
1980
1981 memset(a, 0, sizeof(*a));
1982 for (i = 0; i <= UCHAR_MAX; i++)
1983 for (j = 0; j <= UCHAR_MAX; j++)
1984 if (('A' == i || 'C' == i || 'G' == i || 'T' == i || 'N' == i)
1985 && ('A' == j || 'C' == j || 'G' == j || 'T' == j
1986 || 'N' == j)) {
1987 if (i == 'N' || j == 'N')
1988 a->ssm[i][j] = -25;
1989 else if (i == j)
1990 a->ssm[i][j] = 100;
1991 else
1992 a->ssm[i][j] = -100;
1993 } else
1994 a->ssm[i][j] = INT_MIN;
1995
1996 a->gap = -200;
1997 a->gapl = -200;
1998 a->flag = DPAL_LOCAL;
1999 a->max_gap = 1;
2000 a->fail_stop = 1;
2001 a->check_chars = 0;
2002 a->debug = 0;
2003 a->score_only = 1;
2004 a->force_generic = 0;
2005 a->force_long_generic = 0;
2006 a->force_long_maxgap1 = 0;
2007 }
2008
2009 // char *
2010 // pr_oligo_sequence(sa, o)
2011 // const seq_args *sa;
2012 // const primer_rec *o;
2013 // {
2014 // static char s[MAX_PRIMER_LENGTH+1];
2015 // int seq_len;
2016 // PR_ASSERT(NULL != sa);
2017 // PR_ASSERT(NULL != o);
2018 // seq_len = strlen(sa->sequence);
2019 // PR_ASSERT(o->start + sa->incl_s >= 0);
2020 // PR_ASSERT(o->start + sa->incl_s + o->length <= seq_len);
2021 // _pr_substr(sa->sequence, sa->incl_s + o->start, o->length, s);
2022 // return &s[0];
2023 // }
2024
2025 // char *
2026 // pr_oligo_rev_c_sequence(sa, o)
2027 // const seq_args *sa;
2028 // const primer_rec *o;
2029 // {
2030 // static char s[MAX_PRIMER_LENGTH+1], s1[MAX_PRIMER_LENGTH+1];
2031 // int seq_len, start;
2032 // PR_ASSERT(NULL != sa);
2033 // PR_ASSERT(NULL != o);
2034 // seq_len = strlen(sa->sequence);
2035 // start = sa->incl_s + o->start - o->length + 1;
2036 // PR_ASSERT(start >= 0);
2037 // PR_ASSERT(start + o->length <= seq_len);
2038 // _pr_substr(sa->sequence, start, o->length, s);
2039 // _pr_reverse_complement(s,s1);
2040 // return &s1[0];
2041 // }
2042
2043 /* Calculate self complementarity. */
2044 static void
oligo_compl(h,ha,sa,l,local_args,end_args,local_end_args)2045 oligo_compl(h, ha, sa, l, local_args, end_args, local_end_args)
2046 primer_rec *h;
2047 const primer_args *ha;
2048 seq_args *sa;
2049 oligo_type l;
2050 const dpal_args *local_args, *end_args, *local_end_args;
2051 {
2052 char s[MAX_PRIMER_LENGTH+1], s1[MAX_PRIMER_LENGTH+1];
2053 int j;
2054 short self_any, self_end;
2055
2056 if (OT_INTL == l) {
2057 self_any = ha->io_self_any;
2058 self_end = ha->io_self_end;
2059 } else {
2060 self_any = ha->self_any;
2061 self_end = ha->self_end;
2062 }
2063
2064 j = (OT_LEFT == l || OT_INTL == l) ? h->start : h->start-h->length+1;
2065
2066 _pr_substr(sa->trimmed_seq, j, h->length, s1);
2067 _pr_reverse_complement(s1, s);
2068
2069 h->self_any = align(s1, s, local_args);
2070 if(h->self_any > self_any){
2071 h->ok = OV_SELF_ANY;
2072 if (OT_LEFT == l) {
2073 sa->left_expl.compl_any++;
2074 sa->left_expl.ok--;
2075 }
2076 else if (OT_RIGHT == l) {
2077 sa->right_expl.compl_any++;
2078 sa->right_expl.ok--;
2079 }
2080 else {
2081 sa->intl_expl.compl_any++;
2082 sa->intl_expl.ok--;
2083 }
2084 if (!h->must_use) return;
2085 }
2086
2087 h->self_end = (l != OT_RIGHT) ? align(s1,s,end_args) : align(s,s1,end_args);
2088 if(h->self_end > self_end) {
2089 h->ok = OV_SELF_END;
2090 if (OT_LEFT == l) {
2091 sa->left_expl.compl_end++;
2092 sa->left_expl.ok--;
2093 }
2094 else if (OT_RIGHT == l) {
2095 sa->right_expl.compl_end++;
2096 sa->right_expl.ok--;
2097 }
2098 else {
2099 sa->intl_expl.compl_end++;
2100 sa->intl_expl.ok--;
2101 }
2102 return;
2103 }
2104 }
2105
2106 static void
oligo_mispriming(h,pa,sa,l,align_args,ctx)2107 oligo_mispriming(h, pa, sa, l, align_args, ctx)
2108 primer_rec *h;
2109 const primer_args *pa;
2110 seq_args *sa;
2111 oligo_type l;
2112 const dpal_args *align_args;
2113 Primer3Context * ctx;
2114 {
2115 char
2116 s[MAX_PRIMER_LENGTH+1], /* Will contain the oligo sequence. */
2117 s_tmp[MAX_PRIMER_LENGTH+1], /* Scratch buffer. */
2118 s_r[MAX_PRIMER_LENGTH+1]; /* Will contain s reverse complemented. */
2119
2120 char *oseq, *target, *target_r;
2121 double w;
2122 const seq_lib *lib;
2123 int i;
2124 int first, last; /* Indexes of first and last bases of the oligo in sa->trimmed_seq,
2125 that is, WITHIN THE INCLUDED REGION. */
2126 int first_untrimmed, last_untrimmed;
2127 /* Indexes of first and last bases of the oligo in sa->seq,
2128 that is, WITHIN THE TOTAL SEQUENCE INPUT. */
2129 int min, max, tmp;
2130 int seqlen;
2131 int debug = 0;
2132 int match_length;
2133 short lib_compl;
2134 short tmp_score;
2135 char tmp_char;
2136
2137 if (OT_INTL == l) {
2138 lib = &(pa->io_mishyb_library);
2139 lib_compl = pa->io_repeat_compl;
2140 } else {
2141 lib = &(pa->repeat_lib);
2142 lib_compl = pa->repeat_compl;
2143 }
2144
2145 first = (OT_LEFT == l || OT_INTL == l)
2146 ? h->start
2147 : h->start - h->length + 1;
2148 last = (OT_LEFT == l || OT_INTL == l)
2149 ? h->start + h->length - 1
2150 : h->start;
2151
2152 match_length = h->length;
2153
2154 _pr_substr(sa->trimmed_seq, first, h->length, s_tmp);
2155 _pr_substr(s_tmp, 0, match_length, s);
2156 _pr_reverse_complement(s, s_tmp); /* FIX ME -- is s_tmp needed? */
2157 _pr_substr(s_tmp, 0, match_length, s_r);
2158
2159 /*
2160 * Calculate maximum similarity to sequences from user defined repeat
2161 * library. Compare it with maximum allowed repeat similarity.
2162 */
2163 if(lib->seq_num > 0) {
2164 h->repeat_sim.score =
2165 pr_safe_malloc(lib->seq_num * sizeof(short));
2166 h->repeat_sim.max = h->repeat_sim.min = 0;
2167 max = min = 0;
2168 h->repeat_sim.name = lib->names[0];
2169 for(i = 0; i < lib->seq_num; i++){
2170 if (OT_LEFT == l) w = lib->weight[i] *
2171 align(s, lib->seqs[i], (ctx->lib_local_end_dpal_args));
2172 else if (OT_INTL == l) w = lib->weight[i] *
2173 align(s, lib->seqs[i], (ctx->lib_local_dpal_args));
2174 else w = lib->weight[i] *
2175 align(s_r, lib->rev_compl_seqs[i], (ctx->lib_local_end_dpal_args));
2176
2177 h->repeat_sim.score[i] = w;
2178 if(w > max){
2179 max = w;
2180 h->repeat_sim.max = i;
2181 h->repeat_sim.name = lib->names[i];
2182 }
2183 if(w < min){
2184 min = w;
2185 h->repeat_sim.min = i;
2186 }
2187 if (w > lib_compl) {
2188 h->ok = OV_LIB_SIM;
2189 if (OT_LEFT == l) {
2190 sa->left_expl.repeat_score++;
2191 sa->left_expl.ok--;
2192 }
2193 else if (OT_RIGHT == l) {
2194 sa->right_expl.repeat_score++;
2195 sa->right_expl.ok--;
2196 }
2197 else {
2198 sa->intl_expl.repeat_score++;
2199 sa->intl_expl.ok--;
2200 }
2201 if (!h->must_use) return;
2202 }
2203 }
2204 }
2205
2206 if (_pr_need_template_mispriming(pa) && (l == OT_RIGHT || l == OT_LEFT)) {
2207 /* Calculate maximum similarity to ectopic sites in the template. */
2208
2209 seqlen = strlen(sa->upcased_seq);
2210 first_untrimmed = sa->incl_s + first;
2211 last_untrimmed = sa->incl_s + last;
2212
2213 if (l == OT_LEFT) {
2214 oseq = &s[0];
2215 target = &sa->upcased_seq[0];
2216 target_r = &sa->upcased_seq_r[0];
2217 } else { /* l == OT_RIGHT */
2218 if (debug)
2219 fprintf(stderr, "first_untrimmed = %d, last_untrimmed = %d\n",
2220 first_untrimmed, last_untrimmed);
2221 oseq = &s_r[0];
2222 target = &sa->upcased_seq_r[0];
2223 target_r = &sa->upcased_seq[0];
2224 /* We need to adjust first_untrimmed and last_untrimmed so that
2225 they are correct in the reverse-complemented
2226 sequence.
2227 */
2228 tmp = (seqlen - last_untrimmed) - 1;
2229 last_untrimmed = (seqlen - first_untrimmed) - 1;
2230 first_untrimmed = tmp;
2231 }
2232
2233 /* 1. Align to the template 5' of the oligo. */
2234 tmp_char = target[first_untrimmed];
2235 target[first_untrimmed] = '\0';
2236
2237 /* if (strlen(target) < 3) {
2238 tmp_score = 3;
2239 } else { 2007-07-06, fix was moved to the function align() in this file. */
2240 tmp_score = align(oseq, target, align_args);
2241 /* } */
2242
2243 if (debug) {
2244 if (l == OT_LEFT) fprintf(stderr, "\n************ OLIGO = LEFT\n");
2245 else fprintf(stderr, "\n************ OLIGO = RIGHT\n");
2246 fprintf(stderr, "first_untrimmed = %d, last_untrimmed = %d\n",
2247 first_untrimmed, last_untrimmed);
2248
2249 fprintf(stderr, "5' of oligo: Score %d aligning %s against %s\n\n", tmp_score,
2250 oseq, target);
2251 }
2252
2253 target[first_untrimmed] = tmp_char;
2254
2255 /* 2. Align to the template 3' of the oligo. */
2256 h->template_mispriming
2257 = align(oseq, &target[0] + last_untrimmed + 1, align_args);
2258
2259 if (debug)
2260 fprintf(stderr, "3' of oligo Score %d aligning %s against %s\n\n",
2261 h->template_mispriming, oseq, &target[0] + last_untrimmed + 1);
2262
2263 /* 3. Take the max of 1. and 2. */
2264 if (tmp_score > h->template_mispriming)
2265 h->template_mispriming = tmp_score;
2266
2267 /* 4. Align to the reverse strand of the template. */
2268 h->template_mispriming_r
2269 = align(oseq, target_r, align_args);
2270
2271 if (debug)
2272 fprintf(stderr, "other strand Score %d aligning %s against %s\n\n",
2273 h->template_mispriming_r, oseq, target_r);
2274
2275
2276 if (pa->max_template_mispriming >= 0
2277 && oligo_max_template_mispriming(h) > pa->max_template_mispriming) {
2278 h->ok = OV_TEMPLATE_MISPRIMING;
2279 if (OT_LEFT == l) {
2280 sa->left_expl.template_mispriming++;
2281 sa->left_expl.ok--;
2282 } else if (OT_RIGHT == l) {
2283 sa->right_expl.template_mispriming++;
2284 sa->right_expl.ok--;
2285 } else PR_ASSERT(0); /* Should not get here. */
2286 if (!h->must_use) return;
2287 }
2288 }
2289 }
2290
2291 static int
pair_repeat_sim(h,pa)2292 pair_repeat_sim(h, pa)
2293 primer_pair *h;
2294 const primer_args *pa;
2295 {
2296 int i, n, max, w;
2297 primer_rec *fw, *rev;
2298
2299 fw = h->left;
2300 rev = h->right;
2301
2302 max = 0;
2303 n = pa->repeat_lib.seq_num;
2304 if(n == 0) return 0;
2305 h->rep_name = pa->repeat_lib.names[0] ;
2306 for(i = 0; i < n; i++) {
2307 if((w=(fw->repeat_sim.score[i] +
2308 rev->repeat_sim.score[i])) > max) {
2309 max = w;
2310 h->rep_name = pa->repeat_lib.names[i] ;
2311 }
2312 }
2313 return max;
2314 }
2315
2316 /*
2317 * 's' is the sequence in which to find the stop codon.
2318 * 'start' is the position of a start codon.
2319 *
2320 * There are two modes depending on 'direction'.
2321 *
2322 * If direction is 1:
2323 *
2324 * Return the index of the first stop codon to the right of
2325 * 'start' in 's' (in same frame).
2326 * If there is no such stop codon return -1.
2327 *
2328 * If direction is -1:
2329 *
2330 * If 'start' is negative then return -1.
2331 * Otherwise return the index the first stop codon to left of 'start'.
2332 * If there is no such stop codon return -1.
2333 *
2334 * Note: we don't insist that the start codon be ATG, since in some
2335 * cases the caller will not have the full sequence in 's', nor even
2336 * know the postion of the start codon relative to s.
2337 */
2338 int
find_stop_codon(s,start,direction)2339 find_stop_codon(s, start, direction)
2340 const char* s;
2341 {
2342 const char *p, *q;
2343 int increment = 3 * direction;
2344 int len = strlen(s);
2345
2346 PR_ASSERT(s != NULL);
2347 PR_ASSERT(direction == 1 || direction == -1);
2348 PR_ASSERT(len >= 3);
2349 PR_ASSERT(start <= (len - 3));
2350
2351 if (start < 0) {
2352 if (direction == 1)
2353 while (start < 0) start += increment;
2354 else
2355 return -1;
2356 }
2357
2358 for (p = &s[start];
2359 p >= &s[0]
2360 && *p
2361 && *(p + 1)
2362 && *(p + 2);
2363 p += increment) {
2364 if ('T' != *p && 't' != *p) continue;
2365 q = p + 1;
2366 if ('A' == *q || 'a' == *q) {
2367 q++;
2368 if ('G' == *q || 'g' == *q || 'A' == *q || 'a' == *q)
2369 return p - &s[0];
2370 } else if ('G' == *q || 'g' == *q) {
2371 q++;
2372 if ('A' == *q || 'a' == *q)
2373 return p - &s[0];
2374 }
2375 }
2376 return -1;
2377 }
2378
2379 int
strcmp_nocase(s1,s2)2380 strcmp_nocase(s1, s2)
2381 const char *s1, *s2;
2382 {
2383 // static char M[UCHAR_MAX];
2384 // static int f = 0;
2385 //it is safe to remove static here. But we will get EPIC performance hit.
2386 char M[UCHAR_MAX];
2387 int f = 0;
2388
2389 int i;
2390 const char *p, *q;
2391
2392 if(f != 1){
2393 for(i = 0; i < UCHAR_MAX; i++) M[i] = i;
2394 i = 'a'; M[i] = 'A'; i = 'b'; M[i] = 'B'; i = 'c'; M[i] = 'C';
2395 i = 'A'; M[i] = 'a'; i = 'B'; M[i] = 'b'; i = 'C'; M[i] = 'c';
2396 i = 'd'; M[i] = 'D'; i = 'e'; M[i] = 'E'; i = 'f'; M[i] = 'F';
2397 i = 'D'; M[i] = 'd'; i = 'E'; M[i] = 'e'; i = 'F'; M[i] = 'f';
2398 i = 'g'; M[i] = 'G'; i = 'h'; M[i] = 'H'; i = 'i'; M[i] = 'I';
2399 i = 'G'; M[i] = 'g'; i = 'H'; M[i] = 'h'; i = 'I'; M[i] = 'i';
2400 i = 'k'; M[i] = 'K'; i = 'l'; M[i] = 'L'; i = 'm'; M[i] = 'M';
2401 i = 'K'; M[i] = 'k'; i = 'L'; M[i] = 'l'; i = 'M'; M[i] = 'm';
2402 i = 'n'; M[i] = 'N'; i = 'o'; M[i] = 'O'; i = 'p'; M[i] = 'P';
2403 i = 'N'; M[i] = 'n'; i = 'O'; M[i] = 'o'; i = 'P'; M[i] = 'p';
2404 i = 'q'; M[i] = 'Q'; i = 'r'; M[i] = 'R'; i = 's'; M[i] = 'S';
2405 i = 'Q'; M[i] = 'q'; i = 'R'; M[i] = 'r'; i = 'S'; M[i] = 's';
2406 i = 't'; M[i] = 'T'; i = 'u'; M[i] = 'U'; i = 'v'; M[i] = 'V';
2407 i = 'T'; M[i] = 't'; i = 'U'; M[i] = 'u'; i = 'V'; M[i] = 'v';
2408 i = 'w'; M[i] = 'W'; i = 'x'; M[i] = 'X'; i = 'y'; M[i] = 'Y';
2409 i = 'W'; M[i] = 'w'; i = 'X'; M[i] = 'x'; i = 'Y'; M[i] = 'y';
2410 i = 'z'; M[i] = 'Z'; i = 'Z'; M[i] = 'z'; i = 'j'; M[i] = 'J';
2411 i = 'J'; M[i] = 'j';
2412 f = 1;
2413 }
2414
2415 if(s1 == NULL || s2 == NULL) return 1;
2416 if(strlen(s1) != strlen(s2)) return 1;
2417 p = s1; q = s2;
2418 while(*p != '\0' && *p != '\n' && *q != '\0' && *q != '\n'){
2419 i = *p;
2420 if(*p == *q || M[i] == *q ) {p++; q++; continue;}
2421
2422 return 1;
2423 }
2424 return 0;
2425 }
2426
2427 static void
free_repeat_sim_score(n_f,n_r,n_m,ctx)2428 free_repeat_sim_score(n_f, n_r, n_m, ctx)
2429 int n_f, n_r, n_m;
2430 Primer3Context * ctx;
2431 {
2432 int i;
2433
2434 for(i = 0; i < n_f; i++) if((ctx->f)[i].repeat_sim.score != NULL)
2435 { free((ctx->f)[i].repeat_sim.score); (ctx->f)[i].repeat_sim.score = NULL; }
2436
2437 for(i = 0; i < n_r; i++) if((ctx->r)[i].repeat_sim.score != NULL)
2438 { free((ctx->r)[i].repeat_sim.score); (ctx->r)[i].repeat_sim.score = NULL; }
2439
2440 for(i = 0; i < n_m; i++) if((ctx->mid)[i].repeat_sim.score != NULL)
2441 { free((ctx->mid)[i].repeat_sim.score); (ctx->mid)[i].repeat_sim.score = NULL; }
2442 }
2443
2444 /* Edited by T. Koressaar for lowercase masking. This function checks
2445 if the 3' end of the primer has been masked by lowercase letter.
2446 Function created/Added by Eric Reppo, July 9, 2002
2447 */
2448 static void
check_if_lowercase_masked(position,sequence,h)2449 check_if_lowercase_masked(position, sequence, h)
2450 const int position;
2451 const char *sequence;
2452 primer_rec *h;
2453 {
2454 const char* p = &sequence[position];
2455 if ('a' == *p || 'c' == *p ||'g' == *p || 't' == *p) {
2456 h->ok=OV_GMASKED;
2457 }
2458 }
2459
2460 /* =========================================================== */
2461 /* Various fail-stop wrappers for standard library functions. */
2462
2463 void *
pr_safe_malloc(x)2464 pr_safe_malloc(x)
2465 size_t x;
2466 {
2467 void *r = malloc(x);
2468 if (NULL == r) OOM_ERROR;
2469 return r;
2470 }
2471
2472 void *
pr_safe_realloc(p,x)2473 pr_safe_realloc(p, x)
2474 void *p;
2475 size_t x;
2476 {
2477 void *r = realloc(p, x);
2478 if (NULL == r) OOM_ERROR;
2479 return r;
2480 }
2481
2482 static FILE *
safe_fopen(path,mode)2483 safe_fopen(path, mode)
2484 const char *path, *mode;
2485 {
2486 FILE *r = fopen(path, mode);
2487 if (NULL == r) {
2488 // fprintf(stderr, "%s: unable to open file %s:",
2489 // pr_program_name, path);
2490 perror("");
2491 //exit (-1);
2492 }
2493 return r;
2494 }
2495
2496 /* =========================================================== */
2497
2498 primers_t
runPrimer3(primer_args * pa,seq_args * sa,int * cancelFlag,int * progress)2499 runPrimer3(primer_args *pa, seq_args *sa, int * cancelFlag, int * progress)
2500 {
2501 Primer3Context ctx;
2502 primers_t primers;
2503 // program_args prog_args;
2504 dpal_args *local_args;
2505 dpal_args *end_args;
2506 dpal_args *local_end_args;
2507 int input_found=0;
2508 int n_f, n_r, n_m; /* Will be initialized in pr_choice. */
2509
2510 ctx.f = NULL; ctx.r = NULL; ctx.mid = NULL;
2511 ctx.f_len=0; ctx.r_len=0; ctx.mid_len=0;
2512
2513 /*
2514 * We set up some signal handlers in case someone starts up the program
2515 * from the command line, wonders why nothing is happening, and then kills
2516 * the program.
2517 */
2518 //signal(SIGINT, sig_handler);
2519 //signal(SIGTERM, sig_handler);
2520
2521 local_args = pr_safe_malloc(sizeof(*local_args));
2522 end_args = pr_safe_malloc(sizeof(*end_args));
2523 local_end_args = pr_safe_malloc(sizeof(*local_end_args));
2524 (ctx.lib_local_end_dpal_args) = pr_safe_malloc(sizeof(*(ctx.lib_local_end_dpal_args)));
2525 (ctx.lib_local_dpal_args) = pr_safe_malloc(sizeof(*(ctx.lib_local_dpal_args)));
2526
2527 // memset(&prog_args, 0, sizeof(prog_args));
2528
2529
2530 // prog_args.format_output = 1;
2531 primers.best_pairs.storage_size = primers.best_pairs.num_pairs = 0;
2532
2533 set_dpal_args(local_args);
2534 local_args->flag = DPAL_LOCAL;
2535
2536 set_dpal_args(end_args);
2537 end_args->flag = DPAL_GLOBAL_END;
2538
2539 set_dpal_args(local_end_args);
2540 local_end_args->flag = DPAL_LOCAL_END;
2541
2542 /*
2543 * Read the data from input stream record by record and process it if
2544 * there are no errors.
2545 */
2546
2547 *(ctx.lib_local_dpal_args) = *local_args;
2548 if (pa->lib_ambiguity_codes_consensus) {
2549 PR_ASSERT(dpal_set_ambiguity_code_matrix((ctx.lib_local_dpal_args)));
2550 }
2551
2552 *(ctx.lib_local_end_dpal_args) = *local_end_args;
2553 if (pa->lib_ambiguity_codes_consensus) {
2554 PR_ASSERT(dpal_set_ambiguity_code_matrix((ctx.lib_local_end_dpal_args)));
2555 }
2556
2557 n_f = n_m = n_r = 0;
2558 if (NULL == sa->error.data && NULL == pa->glob_err.data) {
2559 pr_choice(pa, sa, local_args, end_args, local_end_args,
2560 &primers.best_pairs, &n_f, &n_r, &n_m, cancelFlag, progress, &ctx);
2561 }
2562 if (NULL != pa->glob_err.data)
2563 pr_append_new_chunk(&sa->error, pa->glob_err.data);
2564
2565 /*if (pick_pcr_primers == pa->primer_task
2566 || pick_pcr_primers_and_hyb_probe == pa->primer_task) {
2567 if (prog_args.format_output) {
2568 format_pairs(stdout, pa, sa, &best_pairs);
2569 }
2570 else {
2571 boulder_print_pairs(&prog_args, pa, sa, &best_pairs);
2572 }
2573 }
2574 else if(pa->primer_task == pick_left_only) {
2575 if (prog_args.format_output)
2576 format_oligos(stdout, pa, sa, f, n_f, OT_LEFT);
2577 else boulder_print_oligos(pa, sa, n_f, OT_LEFT);
2578 }
2579 else if(pa->primer_task == pick_right_only) {
2580 if (prog_args.format_output)
2581 format_oligos(stdout, pa, sa, r, n_r, OT_RIGHT);
2582 else boulder_print_oligos(pa, sa, n_r, OT_RIGHT);
2583 }
2584 else if(pa->primer_task == pick_hyb_probe_only) {
2585 if(prog_args.format_output)
2586 format_oligos(stdout, pa, sa, mid, n_m, OT_INTL);
2587 else boulder_print_oligos(pa, sa, n_m, OT_INTL);
2588 } */
2589
2590 if(pa->repeat_lib.seq_num > 0 || pa->io_mishyb_library.seq_num > 0)
2591 free_repeat_sim_score(n_f, n_r, n_m, &ctx);
2592 free(local_args);
2593 free(end_args);
2594 free(local_end_args);
2595 free(ctx.lib_local_dpal_args);
2596 free(ctx.lib_local_end_dpal_args);
2597 // free_seq_lib(&pa->repeat_lib);
2598 // free_seq_lib(&pa->io_mishyb_library);
2599 primers.left = ctx.f;
2600 primers.right = ctx.r;
2601 primers.intl = ctx.mid;
2602 return primers;
2603 }
2604
2605 /* End of fail-stop wrappers. */
2606 /* =========================================================== */
2607